From 06a47192ce4ba242a5e97dc88f1a101dbdf4bd1e Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 5 Oct 2024 12:31:22 +1000 Subject: [PATCH] SubFacetSkeletonTriangulations --- scripts/Embedded/SubFacetSkeletons.jl | 193 ++++++++++++++++++++-- scripts/Embedded/differentiable_trians.jl | 26 +-- 2 files changed, 186 insertions(+), 33 deletions(-) diff --git a/scripts/Embedded/SubFacetSkeletons.jl b/scripts/Embedded/SubFacetSkeletons.jl index e3986531..f3a8e222 100644 --- a/scripts/Embedded/SubFacetSkeletons.jl +++ b/scripts/Embedded/SubFacetSkeletons.jl @@ -1,26 +1,189 @@ using GridapEmbedded + using GridapEmbedded.Interfaces +using GridapEmbedded.Interfaces: SubFacetData, SubFacetTriangulation -using Interfaces: SubFacetData, SubFacetTriangulation +using Gridap.Geometry +using Gridap.Geometry: CompositeTriangulation -struct SubFacetSkeletonData{Dp,T,Tn} - subfacets::SubFacetData{Dp,T,Tn} -end +Base.round(a::VectorValue{D,T};kwargs...) where {D,T} = VectorValue{D,T}(round.(a.data;kwargs...)) + +function consistent_facet_to_points( + facet_to_points::Table, point_to_coords::Vector +) + f(pt::Number) = round(pt;sigdigits=12) + f(id::Integer) = f(point_to_coords[id]) -struct SubFacetSkeletonTriangulation{Dc,Dp} <: Triangulation{Dc,Dp} - subfacets :: SubFacetData{Dp} - subgrid :: UnstructuredGrid - function SubFacetSkeletonTriangulation( - subfacets::SubFacetData{Dp} - ) where {Dp} - Dc = Dp-2 - new{Dc,Dp}(subfacets,subgrid) + # Create a list of the unique points composing the facets + npts = length(point_to_coords) + nfaces = length(facet_to_points) + touched = zeros(Bool,npts) + for face in 1:nfaces + pts = view(facet_to_points,face) + touched[pts] .= true end + touched_ids = findall(touched) + unique_ids = unique(f,touched_ids) + + # Create a mapping from the old point ids to the new ones + touched_to_uid = collect(Int32,indexin(f.(touched_ids),f.(unique_ids))) + point_to_uid = extend(touched_to_uid,PosNegPartition(touched_ids,npts)) + + facet_to_uids = Table( + collect(Int32,lazy_map(Reindex(point_to_uid),facet_to_points.data)), + facet_to_points.ptrs + ) + return facet_to_uids, unique_ids end -function SkeletonTriangulation(trian::SubFacetTriangulation{Dc,Dp}) - subfacets = trian.subfacets +function Geometry.get_active_model( + trian::SubFacetTriangulation{Dc,Dp} +) where {Dc,Dp} subgrid = trian.subgrid - return SubFacetSkeletonTriangulation(subfacets,subgrid) + subfacets = trian.subfacets + facet_to_uids, uid_to_point = consistent_facet_to_points( + subfacets.facet_to_points,subfacets.point_to_coords + ) + topo = UnstructuredGridTopology( + subgrid,facet_to_uids,uid_to_point + ) + return UnstructuredDiscreteModel(subgrid,topo,FaceLabeling(topo)) +end + +function generate_ghost_trian( + trian::CompositeTriangulation, + bgmodel::UnstructuredDiscreteModel{Dc} +) where {Dc} + cell_glue = get_glue(trian,Val(Dc)) + + topo = get_grid_topology(bgmodel) + face_to_cell = get_faces(topo,Dc-1,Dc) + cell_to_face = get_faces(topo,Dc,Dc-1) + + n_bgfaces = num_faces(bgmodel,1) + n_faces = num_cells(cell_skeleton) + ghost_faces = zeros(Int32,n_faces) + p_lcell = ones(Int8,n_bgfaces) + m_lcell = ones(Int8,n_bgfaces) + for (i,(p_cell, m_cell)) in enumerate(zip(cell_glue.plus.tface_to_mface,cell_glue.minus.tface_to_mface)) + inter = intersect(cell_to_face[p_cell],cell_to_face[m_cell]) + @assert length(inter) == 1 + face = first(inter) + ghost_faces[i] = face + + nbors = face_to_cell[ghost_faces[i]] + p_lcell[face] = findfirst(==(p_cell),nbors) + m_lcell[face] = findfirst(==(m_cell),nbors) + end + + plus = BoundaryTriangulation(bgmodel,ghost_faces,p_lcell) + minus = BoundaryTriangulation(bgmodel,ghost_faces,m_lcell) + return SkeletonTriangulation(plus,minus) +end +struct SubFacetSkeletonTriangulation{Di,Df,Dp} <: Triangulation{Di,Dp} + cell_skeleton::CompositeTriangulation{Di,Dp} # Interface -> BG Cell pair + face_skeleton::SkeletonTriangulation{Df,Dp} # Interface -> Cut Facet pair + ghost_skeleton::SkeletonTriangulation{Df,Dp} # Ghost Facet -> BG Cell pair + face_trian::SubFacetTriangulation{Df,Dp} # Cut Facet -> BG Cell + face_model::UnstructuredDiscreteModel{Df,Dp} +end + +function SkeletonTriangulation(face_trian::SubFacetTriangulation) + bgmodel = get_background_model(face_trian) + face_model = get_active_model(face_trian) + face_skeleton = SkeletonTriangulation(face_model) + cell_skeleton = CompositeTriangulation(face_trian,face_skeleton) + ghost_skeleton = generate_ghost_trian(cell_skeleton,bgmodel) + return SubFacetSkeletonTriangulation(cell_skeleton,face_skeleton,ghost_skeleton,face_trian,face_model) +end + +function Geometry.get_background_model(t::SubFacetSkeletonTriangulation) + get_background_model(t.cell_skeleton) +end + +function Geometry.get_active_model(t::SubFacetSkeletonTriangulation) + get_active_model(t.cell_skeleton) +end + +function Geometry.get_grid(t::SubFacetSkeletonTriangulation) + get_grid(t.cell_skeleton) +end + +function Geometry.get_glue(ttrian::SubFacetSkeletonTriangulation{Di,Df,Dp},::Val{D}) where {D,Di,Df,Dp} + if D == Dp + get_glue(ttrian.cell_skeleton,Val(D)) # Maps to bgmodel + elseif D == Df + get_glue(ttrian.face_skeleton,Val(D)) # Maps to SubFacetTriangulation + else + @notimplemented + end +end + +# Normal vector to the cut facets , n_∂Ω +function get_subfacet_normal_vector(trian::SubFacetSkeletonTriangulation) + return get_normal_vector(trian.face_skeleton) +end + +# Normal vector to the ghost facets, n_k +function get_ghost_normal_vector(trian::SubFacetSkeletonTriangulation) + n_out = get_normal_vector(trian.ghost_skeleton) + n_in = (-1) * n_out + plus = GenericCellField(CellData.get_data(n_in.plus),trian.cell_skeleton.plus,ReferenceDomain()) + minus = GenericCellField(CellData.get_data(n_in.minus),trian.cell_skeleton.minus,ReferenceDomain()) + return SkeletonPair(plus,minus) +end + +# Tangent vector to a skeleton triangulation +function get_tangent_vector(trian::SkeletonTriangulation{1}) + # TODO: Can we create this in the ReferenceDomain? + # If so, we need to be careful to orient both sides so that their physical + # representations are consistent. + function t(c) + @assert length(c) == 2 + t = c[2] - c[1] + return t/norm(t) + end + edge_coords = get_cell_coordinates(trian) + data = lazy_map(t,CellData.get_data(edge_coords)) + plus = GenericCellField(data,trian.plus,PhysicalDomain()) + minus = GenericCellField(data,trian.minus,PhysicalDomain()) + return SkeletonPair(plus,minus) +end + +# Normal vector to the cut interface, n_S +function get_normal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di} + if Di == 0 + n_S = get_tangent_vector(trian.ghost_skeleton) + elseif Di == 1 + n_k = get_ghost_normal_vector(trian) + t_S = get_tangent_vector(trian.cell_skeleton) + return Operation(cross)(n_k,t_S) # nk = tS x nS -> nS = nk x tS (eq 6.25) + else + @notimplemented + end + return n_S +end + +# Tangent vector to the cut interface, t_S +function get_tangent_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di} + @notimplementedif Di != 1 + n_S = get_normal_vector(trian) + n_k = get_ghost_normal_vector(trian) + return Operation(cross)(n_S,n_k) +end + +# Conormal vectors, m_k = t_S x n_∂Ω +function get_conormal_vector(trian::SubFacetSkeletonTriangulation{Di}) where {Di} + op = Operation(GridapEmbedded._normal_vector) + n_∂Ω = get_subfacet_normal_vector(trian) + if Di == 0 # 2D + m_k = op(n_∂Ω) + elseif Di == 1 # 3D + t_S = get_tangent_vector(trian) + m_k = op(t_S,n_∂Ω) # m_k = t_S x n_∂Ω (eq 6.26) + else + @notimplemented + end + return m_k end diff --git a/scripts/Embedded/differentiable_trians.jl b/scripts/Embedded/differentiable_trians.jl index c01d7d97..80219d5a 100644 --- a/scripts/Embedded/differentiable_trians.jl +++ b/scripts/Embedded/differentiable_trians.jl @@ -114,17 +114,12 @@ function precompute_autodiff_caches( bgmodel = get_background_model(trian) subcells = trian.subcells - cell_to_bgcell = subcells.cell_to_bgcell - cell_to_points = subcells.cell_to_points - point_to_rcoords = subcells.point_to_rcoords - point_to_coords = subcells.point_to_coords - precompute_autodiff_caches( bgmodel, - cell_to_bgcell, - cell_to_points, - point_to_rcoords, - point_to_coords, + subcells.cell_to_bgcell, + subcells.cell_to_points, + subcells.point_to_rcoords, + subcells.point_to_coords, ) end @@ -134,17 +129,12 @@ function precompute_autodiff_caches( bgmodel = get_background_model(trian) subfacets = trian.subfacets - facet_to_bgcell = subfacets.facet_to_bgcell - facet_to_points = subfacets.facet_to_points - point_to_rcoords = subfacets.point_to_rcoords - point_to_coords = subfacets.point_to_coords - precompute_autodiff_caches( bgmodel, - facet_to_bgcell, - facet_to_points, - point_to_rcoords, - point_to_coords, + subfacets.facet_to_bgcell, + subfacets.facet_to_points, + subfacets.point_to_rcoords, + subfacets.point_to_coords, ) end