Added support for overlapping boundary triangulations

src/Arrays/Tables.jl

@@ -365,10 +365,10 @@ Base.IndexStyle(::Type{<:LocalItemFromTable}) = IndexLinear()
 end
 
 """
-    find_local_index(a_to_b, b_to_la_to_a)
+    find_local_index(a_to_b, b_to_la_to_a) -> a_to_la
 """
 function find_local_index(a_to_b, b_to_la_to_a)
-  @notimplemented "find_local_index only implemented for table"
+  @notimplemented "find_local_index only implemented for Table"
 end
 
 function find_local_index(a_to_b, b_to_la_to_a::Table)
@@ -398,6 +398,27 @@ Base.IndexStyle(::Type{<:LocalIndexFromTable}) = IndexStyle(Table)
   return T(UNSET)
 end
 
+"""
+    find_local_index(c_to_a, c_to_b, b_to_la_to_a) -> c_to_lc
+"""
+function find_local_index(
+  c_to_a, c_to_b, b_to_lc_to_a :: Table
+)
+  c_to_lc = fill(Int8(-1),length(c_to_a))
+  for (c,a) in enumerate(c_to_a)
+    b = c_to_b[c]
+    pini = b_to_lc_to_a.ptrs[b]
+    pend = b_to_lc_to_a.ptrs[b+1]-1
+    for (lc,p) in enumerate(pini:pend)
+      if a == b_to_lc_to_a.data[p]
+        c_to_lc[c] = Int8(lc)
+        break
+      end
+    end
+  end
+  return c_to_lc
+end
+
 """
     flatten_partition(a_to_bs::Table,nb::Integer)
     flatten_partition(a_to_bs::Table)

src/Geometry/BoundaryTriangulations.jl

@@ -16,8 +16,8 @@ function FaceToCellGlue(
   cell_grid::Grid,
   face_grid::Grid,
   face_to_bgface::AbstractVector,
-  bgface_to_lcell::AbstractVector)
-
+  bgface_to_lcell::AbstractVector
+)
   D = num_cell_dims(cell_grid)
   cD = num_cell_dims(face_grid)
   bgface_to_cells = get_faces(topo,cD,D)
@@ -55,6 +55,49 @@ function FaceToCellGlue(
     ctype_to_lface_to_ftype)
 end
 
+function OverlappingFaceToCellGlue(
+  topo::GridTopology,
+  cell_grid::Grid,
+  face_grid::Grid,
+  face_to_bgface::AbstractVector,
+  face_to_lcell::AbstractVector
+)
+  Dc = num_cell_dims(cell_grid)
+  Df = num_cell_dims(face_grid)
+  bgface_to_cells = get_faces(topo,Df,Dc)
+  cell_to_bgfaces = get_faces(topo,Dc,Df)
+  cell_to_lface_to_pindex = Table(get_cell_permutations(topo,Df))
+
+  face_to_cells = lazy_map(Reindex(bgface_to_cells), face_to_bgface)
+  face_to_cell  = collect(Int32,lazy_map(getindex,face_to_cells,face_to_lcell))
+  face_to_lface = find_local_index(face_to_bgface,face_to_cell,cell_to_bgfaces)
+
+  f = (p) -> fill(Int8(UNSET), num_faces(p,Df))
+  ctype_to_lface_to_ftype = map(f, get_polytopes(cell_grid))
+  face_to_ftype = get_cell_type(face_grid)
+  cell_to_ctype = get_cell_type(cell_grid)
+
+  _fill_ctype_to_lface_to_ftype!(
+    ctype_to_lface_to_ftype,
+    face_to_cell,
+    face_to_lface,
+    face_to_ftype,
+    cell_to_ctype
+  )
+
+  FaceToCellGlue(
+    face_to_bgface,
+    face_to_lcell,
+    face_to_cell,
+    face_to_lface,
+    face_to_lcell,
+    face_to_ftype,
+    cell_to_ctype,
+    cell_to_lface_to_pindex,
+    ctype_to_lface_to_ftype
+  )
+end
+
 function _fill_ctype_to_lface_to_ftype!(
   ctype_to_lface_to_ftype,
   face_to_cell,

src/Geometry/Geometry.jl

@@ -234,4 +234,6 @@ include("CompressedCellArrays.jl")
 
 include("Polytopal.jl")
 
+include("HybridTriangulations.jl")
+
 end # module

src/Geometry/HybridTriangulations.jl

@@ -0,0 +1,54 @@
+
+struct HybridTriangulation{Dc,Dp,A,B} <: Triangulation{Dc,Dp}
+  cell_trian::A
+  face_trian::B
+  function HybridTriangulation(
+    cell_trian :: BodyFittedTriangulation{Dc,Dp},
+    face_trian :: BoundaryTriangulation{Df,Dp}
+  ) where {Dc,Df,Dp}
+    @assert Dc == Df+1 # TODO: Can we relax this?
+    A, B = typeof(cell_trian), typeof(face_trian)
+    new{Df,Dp,A,B}(cell_trian,face_trian)
+  end
+end
+
+get_background_model(t::HybridTriangulation) = get_background_model(t.face_trian)
+get_grid(t::HybridTriangulation) = get_grid(t.face_trian)
+get_glue(t::HybridTriangulation,::Val{d}) where d = get_glue(t.face_trian,Val(d))
+get_facet_normal(t::HybridTriangulation) = get_facet_normal(t.face_trian)
+
+function HybridTriangulation(model::DiscreteModel{Dc},cell_to_bgcell) where Dc
+  Df = Dc-1 # TODO: Expand to arbitrary face dimension
+
+  # Cell triangulation
+  cell_trian = Triangulation(model,cell_to_bgcell)
+
+  # Get faces around each cell and their local index
+  topo = get_grid_topology(model)
+  bgcell_to_bgfaces = get_faces(topo,Dc,Df)
+  bgface_to_bgcells = get_faces(topo,Df,Dc)
+
+  nfaces = sum(bgcell -> length(view(bgcell_to_bgfaces,bgcell)), cell_to_bgcell)
+  face_to_bgface = Vector{Int}(undef,nfaces)
+  face_to_lcell = Vector{Int8}(undef,nfaces)
+
+  face = 1
+  for (cell,bgcell) in enumerate(cell_to_bgcell)
+    bgfaces = view(bgcell_to_bgfaces,bgcell)
+    for (lface,bgface) in enumerate(bgfaces)
+      lcell = findfirst(x -> x == bgcell, view(bgface_to_bgcells,bgface))
+      face_to_bgface[face] = bgface
+      face_to_lcell[face] = lcell
+      face += 1
+    end
+  end
+
+  # Face triangulation around the selected cells
+  cell_grid = get_grid(model)
+  face_grid = view(Grid(ReferenceFE{Df},model),face_to_bgface)
+  glue = OverlappingFaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,face_to_lcell)
+  trian = BodyFittedTriangulation(model,face_grid,face_to_bgface)
+  face_trian = BoundaryTriangulation(trian,glue)
+
+  return HybridTriangulation(cell_trian,face_trian)
+end

src/Geometry/Triangulations.jl

@@ -138,69 +138,46 @@ function _restrict_model(model,tface_to_mface::IdentityVector)
   model
 end
 
-abstract type TrianFaceModelFaceMapInjectivity end;
-struct Injective    <: TrianFaceModelFaceMapInjectivity end;
-struct NonInjective <: TrianFaceModelFaceMapInjectivity end;
-
 # This is the most basic Triangulation
 # It represents a physical domain built using the faces of a DiscreteModel
-struct BodyFittedTriangulation{Dt,Dp,A,B,C,D<:TrianFaceModelFaceMapInjectivity} <: Triangulation{Dt,Dp}
+struct BodyFittedTriangulation{Dt,Dp,A,B,C} <: Triangulation{Dt,Dp}
   model::A
   grid::B
   tface_to_mface::C
+  injective::Bool
   function BodyFittedTriangulation(model::DiscreteModel,grid::Grid,tface_to_mface)
     Dp = num_point_dims(model)
     @assert Dp == num_point_dims(grid)
     Dt = num_cell_dims(grid)
     A = typeof(model)
     B = typeof(grid)
     C = typeof(tface_to_mface)
-
-    # While we do not have a more definitive solution, we need to distinguish 
-    # between injective and non-injective tface_to_mface maps.
-    # The inverse map, mface_to_tface, relies on PosNegPartition, which fails 
-    # whenever the same mface is the image of more than one tface.
-    # In turn, I have required non-injective mappings for the computation of facet 
-    # integrals on non-conforming cell interfaces.
-    if !(allunique(tface_to_mface))
-      tface_to_mface_injectivity = NonInjective()
-      D = typeof(tface_to_mface_injectivity)
-      new{Dt,Dp,A,B,C,D}(model,grid,tface_to_mface)
-    else 
-      tface_to_mface_injectivity = Injective()
-      D = typeof(tface_to_mface_injectivity)
-      new{Dt,Dp,A,B,C,D}(model,grid,tface_to_mface)
-    end
+    injective = isa(tface_to_mface,IdentityVector) || allunique(tface_to_mface)
+    new{Dt,Dp,A,B,C}(model,grid,tface_to_mface,injective)
   end
 end
 
 get_background_model(trian::BodyFittedTriangulation) = trian.model
 get_grid(trian::BodyFittedTriangulation) = trian.grid
 
-function get_glue(trian::BodyFittedTriangulation{Dt,Dp,A,B,C,Injective},::Val{Dt}) where {Dt,Dp,A,B,C}
+function get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dt}) where Dt
+  n_mfaces = num_faces(trian.model,Dt)
+  n_faces = num_cells(trian)
   tface_to_mface_map = Fill(GenericField(identity),num_cells(trian))
-  if isa(trian.tface_to_mface,IdentityVector) && num_faces(trian.model,Dt) == num_cells(trian)
+  if isa(trian.tface_to_mface,IdentityVector) && (n_faces == n_mfaces)
+    # Case 1: The triangulation spans the whole model injectively
     mface_to_tface = trian.tface_to_mface
+  elseif trian.injective
+    # Case 2: The triangulation spans part of the model injectively
+    mface_to_tface = PosNegPartition(trian.tface_to_mface,Int32(n_mfaces))
   else
-    nmfaces = num_faces(trian.model,Dt)
-    mface_to_tface = PosNegPartition(trian.tface_to_mface,Int32(nmfaces))
+    # Case 3: The triangulation is non-injective, so we cannot map information 
+    # back to the model (1-way glue model -> triangulation)
+    mface_to_tface = nothing
   end
   FaceToFaceGlue(trian.tface_to_mface,tface_to_mface_map,mface_to_tface)
 end
 
-function get_glue(trian::BodyFittedTriangulation{Dt,Dp,A,B,C,NonInjective},::Val{Dt}) where {Dt,Dp,A,B,C}
-  tface_to_mface_map = Fill(GenericField(identity),num_cells(trian))
-  mface_to_tface = nothing
-  # Whenever tface_to_mface is non-injective, we currently avoid the computation of 
-  # mface_to_tface, which relies on PosNegPartition. This is a limitation that we should 
-  # face in the future on those scenarios on which we need mface_to_tface.
-  FaceToFaceGlue(trian.tface_to_mface,tface_to_mface_map,mface_to_tface)
-end
-
-#function get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dm}) where {Dt,Dm}
-#  @notimplemented
-#end
-
 function Base.view(t::BodyFittedTriangulation,ids::AbstractArray)
   model = t.model
   grid = view(t.grid,ids)

test/GeometryTests/HybridTriangulationsTests.jl

@@ -0,0 +1,31 @@
+
+using Gridap
+using Gridap.Geometry
+
+using Gridap.Geometry: HybridTriangulation, get_faces
+using Gridap.CellData: change_domain
+
+model = CartesianDiscreteModel((0,1,0,1),(4,4))
+
+Ω = Triangulation(model)
+Γ = Triangulation(ReferenceFE{1},model)
+Λ = HybridTriangulation(model,[6,7,10,11,14,15])
+
+reffe = ReferenceFE(lagrangian,Float64,1)
+VΩ = FESpace(Ω,reffe)
+VΓ = FESpace(Γ,reffe)
+X = MultiFieldFESpace([VΩ,VΓ])
+
+n = get_normal_vector(Λ)
+
+ΠΛ(u) = change_domain(u,Λ,DomainStyle(u))
+dΛ = Measure(Λ,2)
+mass(u,v) = ∫(u*v)*dΛ
+function mf_mass((_u1,_u2),(_v1,_v2))
+  u1, u2, v1, v2 = ΠΛ(_u1), ΠΛ(_u2), ΠΛ(_v1), ΠΛ(_v2)
+  ∫(u1*v1 + u1*v2 + u2*v1 + u2*v2)*dΛ
+end
+
+assemble_matrix(mass,VΩ,VΩ)
+assemble_matrix(mass,VΓ,VΓ)
+assemble_matrix(mf_mass,X,X)