From 9f7fd1aba8b389f7eb26e41b4048202501b96bc4 Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Tue, 19 Dec 2023 12:34:48 +1100 Subject: [PATCH 1/4] Misc changes required to support facet integrals on non-conforming mesh interfaces --- src/Geometry/BoundaryTriangulations.jl | 89 ++++++++++++++------------ src/Geometry/Triangulations.jl | 8 ++- 2 files changed, 54 insertions(+), 43 deletions(-) diff --git a/src/Geometry/BoundaryTriangulations.jl b/src/Geometry/BoundaryTriangulations.jl index d040d833b..53431876f 100644 --- a/src/Geometry/BoundaryTriangulations.jl +++ b/src/Geometry/BoundaryTriangulations.jl @@ -86,11 +86,10 @@ end struct BoundaryTriangulation{Dc,Dp,A,B} <: Triangulation{Dc,Dp} trian::A glue::B - function BoundaryTriangulation( trian::BodyFittedTriangulation, - glue::FaceToCellGlue) - + glue) + Dc = num_cell_dims(trian) Dp = num_point_dims(trian) A = typeof(trian) @@ -208,9 +207,7 @@ function get_glue(trian::BoundaryTriangulation,::Val{D},::Val{D}) where D FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface) end -function get_facet_normal(trian::BoundaryTriangulation) - - glue = trian.glue +function get_facet_normal(trian::BoundaryTriangulation, boundary_trian_glue::FaceToCellGlue) cell_grid = get_grid(get_background_model(trian.trian)) ## Reference normal @@ -223,24 +220,28 @@ function get_facet_normal(trian::BoundaryTriangulation) lface_pindex_to_n end ctype_lface_pindex_to_nref = map(f, get_reffes(cell_grid)) - face_to_nref = FaceCompressedVector(ctype_lface_pindex_to_nref,glue) + face_to_nref = FaceCompressedVector(ctype_lface_pindex_to_nref,boundary_trian_glue) face_s_nref = lazy_map(constant_field,face_to_nref) # Inverse of the Jacobian transpose cell_q_x = get_cell_map(cell_grid) cell_q_Jt = lazy_map(∇,cell_q_x) cell_q_invJt = lazy_map(Operation(pinvJt),cell_q_Jt) - face_q_invJt = lazy_map(Reindex(cell_q_invJt),glue.face_to_cell) + face_q_invJt = lazy_map(Reindex(cell_q_invJt),boundary_trian_glue.face_to_cell) # Change of domain D = num_cell_dims(cell_grid) - glue = get_glue(trian,Val(D)) - face_s_q = glue.tface_to_mface_map + boundary_trian_glue = get_glue(trian,Val(D)) + face_s_q = boundary_trian_glue.tface_to_mface_map face_s_invJt = lazy_map(∘,face_q_invJt,face_s_q) face_s_n = lazy_map(Broadcasting(Operation(push_normal)),face_s_invJt,face_s_nref) Fields.MemoArray(face_s_n) end +function get_facet_normal(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A} + glue = trian.glue + get_facet_normal(trian, glue) +end function push_normal(invJt,n) v = invJt⋅n m = sqrt(inner(v,v)) @@ -251,42 +252,46 @@ function push_normal(invJt,n) end end -function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation) - d = num_cell_dims(trian) - cell_grid = get_grid(get_background_model(trian.trian)) - polytopes = map(get_polytope, get_reffes(cell_grid)) - cell_to_ctype = trian.glue.cell_to_ctype - ctype_to_lvertex_to_qcoords = map(get_vertex_coordinates, polytopes) - ctype_to_lface_to_lvertices = map((p)->get_faces(p,d,0), polytopes) - ctype_to_lface_to_pindex_to_perm = map( (p)->get_face_vertex_permutations(p,d), polytopes) - - P = eltype(eltype(ctype_to_lvertex_to_qcoords)) - D = num_components(P) - T = eltype(P) - ctype_to_lface_to_pindex_to_qcoords = Vector{Vector{Vector{Point{D,T}}}}[] - - for (ctype, lface_to_pindex_to_perm) in enumerate(ctype_to_lface_to_pindex_to_perm) - lvertex_to_qcoods = ctype_to_lvertex_to_qcoords[ctype] - lface_to_pindex_to_qcoords = Vector{Vector{Point{D,T}}}[] - for (lface, pindex_to_perm) in enumerate(lface_to_pindex_to_perm) - cfvertex_to_lvertex = ctype_to_lface_to_lvertices[ctype][lface] - nfvertices = length(cfvertex_to_lvertex) - pindex_to_qcoords = Vector{Vector{Point{D,T}}}(undef,length(pindex_to_perm)) - for (pindex, cfvertex_to_ffvertex) in enumerate(pindex_to_perm) - ffvertex_to_qcoords = zeros(Point{D,T},nfvertices) - for (cfvertex, ffvertex) in enumerate(cfvertex_to_ffvertex) - lvertex = cfvertex_to_lvertex[cfvertex] - qcoords = lvertex_to_qcoods[lvertex] - ffvertex_to_qcoords[ffvertex] = qcoords - end - pindex_to_qcoords[pindex] = ffvertex_to_qcoords +function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation,glue) + d = num_cell_dims(trian) + cell_grid = get_grid(get_background_model(trian.trian)) + polytopes = map(get_polytope, get_reffes(cell_grid)) + cell_to_ctype = glue.cell_to_ctype + ctype_to_lvertex_to_qcoords = map(get_vertex_coordinates, polytopes) + ctype_to_lface_to_lvertices = map((p)->get_faces(p,d,0), polytopes) + ctype_to_lface_to_pindex_to_perm = map( (p)->get_face_vertex_permutations(p,d), polytopes) + + P = eltype(eltype(ctype_to_lvertex_to_qcoords)) + D = num_components(P) + T = eltype(P) + ctype_to_lface_to_pindex_to_qcoords = Vector{Vector{Vector{Point{D,T}}}}[] + + for (ctype, lface_to_pindex_to_perm) in enumerate(ctype_to_lface_to_pindex_to_perm) + lvertex_to_qcoods = ctype_to_lvertex_to_qcoords[ctype] + lface_to_pindex_to_qcoords = Vector{Vector{Point{D,T}}}[] + for (lface, pindex_to_perm) in enumerate(lface_to_pindex_to_perm) + cfvertex_to_lvertex = ctype_to_lface_to_lvertices[ctype][lface] + nfvertices = length(cfvertex_to_lvertex) + pindex_to_qcoords = Vector{Vector{Point{D,T}}}(undef,length(pindex_to_perm)) + for (pindex, cfvertex_to_ffvertex) in enumerate(pindex_to_perm) + ffvertex_to_qcoords = zeros(Point{D,T},nfvertices) + for (cfvertex, ffvertex) in enumerate(cfvertex_to_ffvertex) + lvertex = cfvertex_to_lvertex[cfvertex] + qcoords = lvertex_to_qcoods[lvertex] + ffvertex_to_qcoords[ffvertex] = qcoords end - push!(lface_to_pindex_to_qcoords,pindex_to_qcoords) + pindex_to_qcoords[pindex] = ffvertex_to_qcoords end - push!(ctype_to_lface_to_pindex_to_qcoords,lface_to_pindex_to_qcoords) + push!(lface_to_pindex_to_qcoords,pindex_to_qcoords) end + push!(ctype_to_lface_to_pindex_to_qcoords,lface_to_pindex_to_qcoords) + end + + FaceCompressedVector(ctype_to_lface_to_pindex_to_qcoords,glue) +end - FaceCompressedVector(ctype_to_lface_to_pindex_to_qcoords,trian.glue) +function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A} + _compute_face_to_q_vertex_coords(trian,trian.glue) end struct FaceCompressedVector{T,G<:FaceToCellGlue} <: AbstractVector{T} diff --git a/src/Geometry/Triangulations.jl b/src/Geometry/Triangulations.jl index ed82f0f70..dffc2eb31 100644 --- a/src/Geometry/Triangulations.jl +++ b/src/Geometry/Triangulations.jl @@ -147,7 +147,13 @@ get_background_model(trian::BodyFittedTriangulation) = trian.model get_grid(trian::BodyFittedTriangulation) = trian.grid function get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dt}) where Dt - tface_to_mface = trian.tface_to_mface + # union(...) here is required for those cases in which an mface in + # trian.tface_to_mface is listed more than once. Otherwise, the + # constructor of PosNegPartition fails because it does not admit + # the same mface to be the image of more than one tface. + # In turn, I have required this for the computation of facet integrals + # on non-conforming cell interfaces. + tface_to_mface = union(trian.tface_to_mface) tface_to_mface_map = Fill(GenericField(identity),num_cells(trian)) if isa(tface_to_mface,IdentityVector) && num_faces(trian.model,Dt) == num_cells(trian) mface_to_tface = tface_to_mface From 77348aef88b4b304f4e5a726eacb5bbfa7369f11 Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Tue, 19 Dec 2023 13:07:55 +1100 Subject: [PATCH 2/4] fixed get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dt}) where Dt so that it does not return a new instance whenever there are not repeated elements in trian.tface_to_mface --- src/Geometry/Triangulations.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/Geometry/Triangulations.jl b/src/Geometry/Triangulations.jl index dffc2eb31..7cc169a0e 100644 --- a/src/Geometry/Triangulations.jl +++ b/src/Geometry/Triangulations.jl @@ -147,13 +147,17 @@ get_background_model(trian::BodyFittedTriangulation) = trian.model get_grid(trian::BodyFittedTriangulation) = trian.grid function get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dt}) where Dt - # union(...) here is required for those cases in which an mface in + # unique(...) here is required for those cases in which an mface in # trian.tface_to_mface is listed more than once. Otherwise, the # constructor of PosNegPartition fails because it does not admit # the same mface to be the image of more than one tface. # In turn, I have required this for the computation of facet integrals # on non-conforming cell interfaces. - tface_to_mface = union(trian.tface_to_mface) + if !(allunique(trian.tface_to_mface)) + tface_to_mface = unique(trian.tface_to_mface) + else + tface_to_mface = trian.tface_to_mface + end tface_to_mface_map = Fill(GenericField(identity),num_cells(trian)) if isa(tface_to_mface,IdentityVector) && num_faces(trian.model,Dt) == num_cells(trian) mface_to_tface = tface_to_mface From 7505760e761161d2ce9a85ed7d6b533ae651f1cd Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" <38347633+amartinhuertas@users.noreply.github.com> Date: Wed, 20 Dec 2023 14:21:56 +1100 Subject: [PATCH 3/4] Update NEWS.md --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index 6e080ebfe..82f0d517b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - Changed how `allocate_vector` works. Now it only allocates, instead of allocating+initialising to zero. Since PR[#963](https://github.com/gridap/Gridap.jl/pull/963). +- Misc changes required to support facet integration on non-conforming meshes. These changes do not involve methods of the public API. Since PR[#967](https://github.com/gridap/Gridap.jl/pull/967) ## [0.17.21] - 2023-12-04 From dc1f4cf2ae3df2147c174185aa44b93842df7e9b Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Wed, 28 Aug 2024 22:28:01 +1000 Subject: [PATCH 4/4] Updated NEWS Bump version --- NEWS.md | 7 ++++++- Project.toml | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 3d4a8c76e..8321b6566 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.18.5] - 2024-08-28 + +### Changed + +- Misc changes required to support facet integration on non-conforming meshes. These changes do not involve methods of the public API. Since PR[#967](https://github.com/gridap/Gridap.jl/pull/967) + ## [0.18.4] - 2024-08-09 ### Changed @@ -77,7 +83,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - Changed how `allocate_vector` works. Now it only allocates, instead of allocating+initialising to zero. Since PR[#963](https://github.com/gridap/Gridap.jl/pull/963). -- Misc changes required to support facet integration on non-conforming meshes. These changes do not involve methods of the public API. Since PR[#967](https://github.com/gridap/Gridap.jl/pull/967) ### Fixed diff --git a/Project.toml b/Project.toml index 69a423369..310df8f9b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Gridap" uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e" authors = ["Santiago Badia ", "Francesc Verdugo ", "Alberto F. Martin "] -version = "0.18.4" +version = "0.18.5" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"