Bugfix: fix topological consistency for polytopal models

Author: JordiManyer

src/Adaptivity.jl

@@ -909,6 +909,7 @@ function Adaptivity.coarsen(fmodel::DistributedDiscreteModel, ptopo::Distributed
   ctopo = GridapDistributed.DistributedGridTopology(
     map(Geometry.PolytopalGridTopology, cnode_coords, connectivity, polys), face_gids
   )
+  _setup_consistent_faces!(ctopo)
   labels = Geometry.FaceLabeling(ctopo)
   cmodels = map(local_views(ctopo), local_views(labels)) do ctopo, labels
     cgrid = Geometry.PolytopalGrid(ctopo)

src/Geometry.jl

@@ -99,13 +99,72 @@ function _setup_face_gids!(topo::DistributedGridTopology{Dc},dim) where {Dc}
   end
 end
 
+# In some cases, the orientation of locally computed faces is NOT consistent. 
+# The following functions can be used to check for consistent orientation and fix it.
+function _setup_consistent_faces!(topo::DistributedGridTopology)
+  # Setting up consistent face-to-vertex maps should be enough
+  # to guarantee consistent face orientation if it is done before 
+  # any other face-to-face map is setup. So we should call this function 
+  # just after creating the new models.
+  D = num_cell_dims(topo)
+  for dimfrom in 1:D-1
+    _setup_consistent_faces!(topo, dimfrom, 0)
+  end
+end
+
+function _setup_consistent_faces!(topo::DistributedGridTopology, dimfrom::Integer, dimto::Integer)
+  @check 0 <= dimto <= dimfrom <= num_cell_dims(topo)
+  gids_from = partition(get_face_gids(topo, dimfrom))
+  gids_to   = partition(get_face_gids(topo, dimto))
+  lfrom_to_gto = map(local_views(topo), gids_to) do topo, gids_to
+    lfrom_to_lto = get_faces(topo, dimfrom, dimto)
+    to_global!(lfrom_to_lto.data, gids_to)
+    JaggedArray(lfrom_to_lto.data, lfrom_to_lto.ptrs)
+  end
+  wait(consistent!(PVector(lfrom_to_gto, gids_from)))
+  map(lfrom_to_gto, gids_to) do lfrom_to_gto, gids_to
+    to_local!(lfrom_to_gto.data, gids_to)
+  end
+  return nothing
+end
+
+function isconsistent_faces(topo::DistributedGridTopology)
+  D = num_cell_dims(topo)
+  for dimfrom in 1:D-1
+    for dimto in 0:dimfrom-1
+      !isconsistent_faces(topo, dimfrom, dimto) && return false
+    end
+  end
+  return true
+end
+
+function isconsistent_faces(topo::DistributedGridTopology, dimfrom::Integer, dimto::Integer)
+  @check 0 <= dimto <= dimfrom <= num_cell_dims(topo)
+  gids_from = partition(get_face_gids(topo, dimfrom))
+  gids_to   = partition(get_face_gids(topo, dimto))
+
+  lfrom_to_lto = map(local_views(topo)) do topo
+    get_faces(topo, dimfrom, dimto)
+  end
+  lfrom_to_gto = map(lfrom_to_lto, gids_to) do lfrom_to_lto, gids_to
+    lto_gto = local_to_global(gids_to)
+    JaggedArray(lto_gto[lfrom_to_lto.data],lfrom_to_lto.ptrs)
+  end
+  wait(consistent!(PVector(lfrom_to_gto, gids_from)))
+  isconsistent = map(lfrom_to_lto, lfrom_to_gto, gids_to) do lfrom_to_lto, lfrom_to_gto, gids_to
+    gto_to_lto = global_to_local(gids_to)
+    lfrom_to_lto.data == gto_to_lto[lfrom_to_gto.data]
+  end
+  return reduce(&, isconsistent)
+end
+
 function Geometry.get_isboundary_face(topo::DistributedGridTopology, d::Integer)
   face_gids = get_face_gids(topo, d)
   is_local_boundary = map(local_views(topo)) do topo
     get_isboundary_face(topo,d)
   end
-  t = assemble!(&,PVector(is_local_boundary, partition(face_gids))) |> fetch |> consistent!
-  is_global_boundary = partition(fetch(t))
+  t = assemble!(&,PVector(is_local_boundary, partition(face_gids)))
+  is_global_boundary = partition(fetch(consistent!(fetch(t))))
   return is_global_boundary
 end
 
@@ -515,10 +574,12 @@ end
 # PolytopalDiscreteModel
 
 function Geometry.PolytopalDiscreteModel(model::GenericDistributedDiscreteModel)
-  return GenericDistributedDiscreteModel(
+  pmodel = GenericDistributedDiscreteModel(
     map(Geometry.PolytopalDiscreteModel,local_views(model)),
     get_cell_gids(model)
   )
+  _setup_consistent_faces!(get_grid_topology(pmodel))
+  return pmodel
 end
 
 # Simplexify

test/GeometryTests.jl

@@ -6,8 +6,75 @@ using PartitionedArrays
 using LinearAlgebra
 using Test
 
-function main(distribute,parts)
+using Gridap.Geometry
+
+function test_local_part_face_labelings_consistency(lmodel::CartesianDiscreteModel{D},gids,gmodel) where {D}
+  local_topology         = lmodel.grid_topology
+  global_topology        = gmodel.grid_topology
+  local_labelings        = lmodel.face_labeling
+  global_labelings       = gmodel.face_labeling
+  l_d_to_dface_to_entity = local_labelings.d_to_dface_to_entity
+  g_d_to_dface_to_entity = global_labelings.d_to_dface_to_entity
+  loc_to_glo             = local_to_global(gids)
+  #traverse local cells
+  for cell_lid=1:num_cells(lmodel)
+    cell_gid=loc_to_glo[cell_lid]
+    for d = 0:D-1
+      local_cell_to_faces = local_topology.n_m_to_nface_to_mfaces[D+1,d+1]
+      global_cell_to_faces = global_topology.n_m_to_nface_to_mfaces[D+1,d+1]
+      la = local_cell_to_faces.ptrs[cell_lid]
+      lb = local_cell_to_faces.ptrs[cell_lid+1]
+      ga = global_cell_to_faces.ptrs[cell_gid]
+      gb = global_cell_to_faces.ptrs[cell_gid+1]
+      @assert (lb-la)==(gb-ga)
+      for i = 0:lb-la-1
+        face_lid = local_cell_to_faces.data[la+i]
+        face_gid = global_cell_to_faces.data[ga+i]
+        local_entity = l_d_to_dface_to_entity[d+1][face_lid]
+        global_entity = g_d_to_dface_to_entity[d+1][face_gid]
+        if (local_entity != global_entity)
+          return false
+        end
+      end
+    end
+  end
+  return true
+end
+
+function test_model(model)
+  D = num_cell_dims(model)
+
+  grid = get_grid(model)
+  labels = get_grid(model)
+
+  topo = get_grid_topology(model)
+  GridapDistributed.isconsistent_faces(topo)
+
+  for d in 0:D
+    fgids = get_face_gids(model,d)
+    trian = Triangulation(no_ghost,ReferenceFE{d},model)
+    @test num_cells(trian) == length(fgids)
+    trian = Triangulation(with_ghost,ReferenceFE{d},model)
+    @test num_cells(trian) == length(fgids)
+  end
+
+  Ω = Triangulation(with_ghost,model)
+  @test num_cells(Ω) == num_cells(model)
 
+  Ω = Triangulation(no_ghost,model)
+  @test num_cells(Ω) == num_cells(model)
+
+  Γ = Boundary(with_ghost,model,tags="boundary")
+  nbfacets = num_cells(Γ)
+
+  Γ = Boundary(no_ghost,model,tags="boundary")
+  @test num_cells(Γ) == nbfacets
+
+  return true
+end
+
+function main_cartesian(distribute,parts)
+  ranks = distribute(LinearIndices((prod(parts),)))
   output = mkpath(joinpath(@__DIR__,"output"))
 
   if length(parts) == 2
@@ -18,14 +85,12 @@ function main(distribute,parts)
     cells = (4,4,4)
   end
 
-  ranks = distribute(LinearIndices((prod(parts),)))
-
-
   model = CartesianDiscreteModel(ranks,parts,domain,cells)
   writevtk(model,joinpath(output,"model"))
 
-  @test num_cells(model)==prod(cells)
-  @test num_vertices(model)==prod(cells .+ 1)
+  @test test_model(model)
+  @test num_cells(model) == prod(cells)
+  @test num_vertices(model) == prod(cells .+ 1)
   @test num_cell_dims(model) == length(cells)
   @test num_point_dims(model) == length(cells)
 
@@ -54,25 +119,6 @@ function main(distribute,parts)
     @test test_local_part_face_labelings_consistency(lmodel,gids,gmodel)
   end
 
-  grid = get_grid(model)
-  labels = get_grid(model)
-
-  Ω = Triangulation(with_ghost,model)
-  writevtk(Ω,joinpath(output,"Ω"))
-  @test num_cells(Ω) == num_cells(model)
-
-  Ω = Triangulation(no_ghost,model)
-  writevtk(Ω,joinpath(output,"Ω"))
-  @test num_cells(Ω) == num_cells(model)
-
-  Γ = Boundary(with_ghost,model,tags="boundary")
-  writevtk(Γ,joinpath(output,"Γ"))
-  nbfacets = num_cells(Γ)
-
-  Γ = Boundary(no_ghost,model,tags="boundary")
-  writevtk(Γ,joinpath(output,"Γ"))
-  @test num_cells(Γ) == nbfacets
-
   function is_in(coords)
     R = 1.6
     n = length(coords)
@@ -111,40 +157,48 @@ function main(distribute,parts)
 
   # Multiple ghost layers
   model = CartesianDiscreteModel(ranks,parts,domain,cells;ghost=map(i->2,parts))
+  @test test_model(model)
+
+  # Unstructured conversion
+  model = Geometry.UnstructuredDiscreteModel(
+    CartesianDiscreteModel(ranks,parts,domain,cells)
+  )
+  @test test_model(model)
+
+  # Simplexify
+  model = simplexify(
+    CartesianDiscreteModel(ranks,parts,domain,cells)
+  )
+  @test test_model(model)
 
 end
 
-function test_local_part_face_labelings_consistency(lmodel::CartesianDiscreteModel{D},gids,gmodel) where {D}
-   local_topology         = lmodel.grid_topology
-   global_topology        = gmodel.grid_topology
-   local_labelings        = lmodel.face_labeling
-   global_labelings       = gmodel.face_labeling
-   l_d_to_dface_to_entity = local_labelings.d_to_dface_to_entity
-   g_d_to_dface_to_entity = global_labelings.d_to_dface_to_entity
-   loc_to_glo             = local_to_global(gids)
-   #traverse local cells
-   for cell_lid=1:num_cells(lmodel)
-        cell_gid=loc_to_glo[cell_lid]
-        for d=0:D-1
-             local_cell_to_faces = local_topology.n_m_to_nface_to_mfaces[D+1,d+1]
-             global_cell_to_faces = global_topology.n_m_to_nface_to_mfaces[D+1,d+1]
-             la = local_cell_to_faces.ptrs[cell_lid]
-             lb = local_cell_to_faces.ptrs[cell_lid+1]
-             ga = global_cell_to_faces.ptrs[cell_gid]
-             gb = global_cell_to_faces.ptrs[cell_gid+1]
-             @assert (lb-la)==(gb-ga)
-             for i=0:lb-la-1
-                 face_lid = local_cell_to_faces.data[la+i]
-                 face_gid = global_cell_to_faces.data[ga+i]
-                 local_entity = l_d_to_dface_to_entity[d+1][face_lid]
-                 global_entity = g_d_to_dface_to_entity[d+1][face_gid]
-                 if (local_entity != global_entity)
-                     return false
-                 end
-             end
-        end
-   end
-   return true
+function main_polytopal(distribute,parts)
+  ranks = distribute(LinearIndices((prod(parts),)))
+
+  if length(parts) == 2
+    domain = (0,4,0,4)
+    cells = (4,4)
+  elseif length(parts) == 3
+    domain = (0,4,0,4,0,4)
+    cells = (4,4,4)
+  end
+
+  model = Geometry.PolytopalDiscreteModel(
+    CartesianDiscreteModel(ranks,parts,domain,cells)
+  )
+  @test test_model(model)
+
+  model = Geometry.PolytopalDiscreteModel(
+    simplexify(CartesianDiscreteModel(ranks,parts,domain,cells))
+  )
+  @test test_model(model)
+
+end
+
+function main(distribute,parts)
+  main_cartesian(distribute,parts)
+  main_polytopal(distribute,parts)
 end
 
 end # module

test/PolytopalCoarseningTests.jl

@@ -1,5 +1,6 @@
 module PolytopalCoarseningTests
 
+using Test
 using Gridap
 using GridapDistributed, PartitionedArrays
 
@@ -28,7 +29,30 @@ function main(distribute, np)
   ptopo = Geometry.PatchTopology(get_grid_topology(fmodel), patch_cells)
 
   cmodel, glues = Adaptivity.coarsen(fmodel,ptopo; return_glue=true)
+  @test GridapDistributed.isconsistent_faces(get_grid_topology(cmodel))
   #writevtk(cmodel, "tmp/cmodel")
 end
 
+# ranks = collect(1:2)
+# 
+# fmodel_serial = Geometry.PolytopalDiscreteModel(CartesianDiscreteModel((0,1,0,1),(4,4)))
+# cmodel_serial = Adaptivity.coarsen(
+#   fmodel_serial, Geometry.PatchTopology(
+#     get_grid_topology(fmodel_serial), Table([[1,2,5,6],[3,4,7,8],[9,10,13,14],[11,12,15,16]])
+#   )
+# )
+# fmodel_good = DiscreteModel(ranks,fmodel_serial,[1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,2])
+# cmodel_good = DiscreteModel(ranks,cmodel_serial,[1,2,1,2])
+# 
+# fmodel = Geometry.PolytopalDiscreteModel(CartesianDiscreteModel(ranks,(2,1),(0,1,0,1),(4,4)))
+# cmodel = Adaptivity.coarsen(
+#   fmodel, Geometry.PatchTopology(
+#     get_grid_topology(fmodel),[Table([[1,2,4,5],[7,8,10,11]]),Table([[2,3,5,6],[8,9,11,12]])]
+#   )
+# )
+# 
+# GridapDistributed.isconsistent_faces(get_grid_topology(fmodel))
+# GridapDistributed.isconsistent_faces(get_grid_topology(cmodel))
+
+
 end # module