Merge pull request #91 from pmartorell/distributed_agfem_failing_case

Fix distributed AgFEM corner case

NEWS.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 - Unused `name` keywork argument in `square` and `quadrilateral` analytical geometries. Since PR [#86](https://github.com/gridap/GridapEmbedded.jl/pull/86).
+- Gluing of the remote root cells to the local+ghost mesh. Since PR [#91](https://github.com/gridap/GridapEmbedded.jl/pull/91).
 
 ## [0.9.3] - 2024-05-20
 

src/AgFEM/AgFEMSpaces.jl

@@ -76,9 +76,9 @@ function _setup_agfem_constraints(
       if dof > 0
         fdof = dof
         acell_dof = fdof_to_acell[fdof]
+        fdof_to_isagg[fdof] &= iscut
         if acell_dof == 0 || gcell > acell_to_gcell[acell_dof]
           fdof_to_acell[fdof] = acell
-          fdof_to_isagg[fdof] = iscut && fdof_to_isagg[fdof]
           fdof_to_ldof[fdof] = ldof
          end
       end

src/Distributed/Distributed.jl

@@ -11,6 +11,8 @@ using Gridap.Geometry
 using Gridap.Helpers
 using Gridap.ReferenceFEs
 
+using PartitionedArrays: VectorFromDict
+
 using GridapEmbedded.CSG
 using GridapEmbedded.LevelSetCutters
 using GridapEmbedded.Interfaces
@@ -25,6 +27,8 @@ using GridapEmbedded.AgFEM: AggregateCutCellsByThreshold
 using GridapEmbedded.MomentFittedQuadratures: MomentFitted
 using Gridap.Geometry: AppendedTriangulation
 using Gridap.Geometry: get_face_to_parent_face
+using Gridap.FESpaces: FESpaceWithLinearConstraints
+using Gridap.FESpaces: _dof_to_DOF, _DOF_to_dof
 using GridapDistributed: DistributedDiscreteModel
 using GridapDistributed: DistributedTriangulation
 using GridapDistributed: DistributedFESpace

src/Distributed/DistributedAgFEM.jl

@@ -18,9 +18,7 @@ function AgFEMSpace(
   trian = add_ghost_cells(trian)
   trian_gids = generate_cell_gids(trian)
   cell_to_cellin = _active_aggregates(bgcell_to_bgcellin)
-  cell_to_ldofs = map(get_cell_dof_ids,spaces)
-  cell_to_ldofs = map(i->map(sort,i),cell_to_ldofs)
-  _remove_improper_cell_ldofs!(cell_to_ldofs,cell_to_cellin)
+  cell_to_ldofs = cell_ldof_to_mdof(spaces,cell_to_cellin)
   nldofs = map(num_free_dofs,spaces)
   gids = generate_gids(trian_gids,cell_to_ldofs,nldofs)
   vector_type = _find_vector_type(spaces,gids)
@@ -296,6 +294,39 @@ function _get_cell_measure(trian1::Triangulation,trian2::Triangulation)
   end
 end
 
+function merge_nodes(model::DistributedDiscreteModel)
+  node_gids = get_face_gids(model,0)
+  cell_gids = get_cell_gids(model)
+  models = map(local_views(model),local_views(node_gids)) do model,node_ids
+    merge_nodes(model,node_ids)
+  end
+  DistributedDiscreteModel(models,cell_gids)
+end
+
+function merge_nodes(model::DiscreteModel,ids)
+  l_to_g = local_to_global(ids)
+  n_global = length(global_to_local(ids))
+  n_local = length(l_to_g)
+  g_to_l = VectorFromDict(reverse(l_to_g),reverse(1:length(l_to_g)),n_global)
+  l_to_lparent = g_to_l[l_to_g]
+  lnew_to_l = findall(map(==,l_to_lparent,1:n_local))
+  l_to_lnew = zeros(Int,n_local)
+  l_to_lnew[lnew_to_l] = 1:length(lnew_to_l)
+  g_to_lnew = map(l->iszero(l) ? l : l_to_lnew[l], g_to_l)
+  l_to_lnew = g_to_lnew[l_to_g]
+
+  grid = get_grid(model)
+  coords = get_node_coordinates(grid)
+  conn = get_cell_node_ids(grid)
+  reffes = get_reffes(grid)
+  ctypes = get_cell_type(grid)
+
+  coords = coords[lnew_to_l]
+  conn = map(Broadcasting(Reindex(l_to_lnew)),conn) |> Table
+  grid = UnstructuredGrid(coords,conn,reffes,ctypes)
+  UnstructuredDiscreteModel(grid)
+end
+
 function add_remote_cells(model::DistributedDiscreteModel,remote_cells,remote_parts)
   # Send remote gids to owners
   snd_ids = remote_parts
@@ -343,7 +374,7 @@ function add_remote_cells(model::DistributedDiscreteModel,remote_cells,remote_pa
   grids = map(lazy_append,lgrids,rgrids)
   models = map(UnstructuredDiscreteModel,grids)
   agids = add_remote_ids(gids,remote_cells,remote_parts)
-  DistributedDiscreteModel(models,agids)
+  DistributedDiscreteModel(models,agids) |> merge_nodes
 end
 
 function add_remote_aggregates(model::DistributedDiscreteModel,aggregates,aggregate_owner)
@@ -436,29 +467,44 @@ function _active_aggregates(bgcell_to_bgcellin)
   bgcell_to_acell[ acell_to_bgcellin ]
 end
 
-function _remove_improper_cell_ldofs!(
-  cell_to_ldofs::AbstractVector{<:AbstractVector{<:AbstractVector}},
-  bgcell_to_bgcellin::AbstractVector{<:AbstractVector})
+function cell_ldof_to_mdof(
+  spaces::AbstractArray{<:FESpace},
+  cell_to_cellin::AbstractArray{<:AbstractVector})
 
-  map(_remove_improper_cell_ldofs!,cell_to_ldofs,bgcell_to_bgcellin)
+  map(cell_ldof_to_mdof,spaces,cell_to_cellin)
 end
 
-
-function _remove_improper_cell_ldofs!(cell_to_ldofs,cell_to_cellin)
-  for cell in 1:length(cell_to_ldofs)
-    cell_to_cellin[cell] != cell || continue
-    cell_to_ldofs[cell] = empty!(cell_to_ldofs[cell])
+function cell_ldof_to_mdof(
+  space::FESpaceWithLinearConstraints,
+  cell_to_cellin::AbstractVector)
+
+  DOF_to_mDOFs = space.DOF_to_mDOFs
+  cell_ldof_to_dof = space.cell_to_ldof_to_dof
+  n_fdofs = space.n_fdofs
+  n_fmdofs = space.n_fmdofs
+  cell_ldof_to_mdof = map(cell_ldof_to_dof) do ldof_to_dof
+    map(ldof_to_dof) do dof
+      DOF = _dof_to_DOF(dof,n_fdofs)
+      mDOFs = DOF_to_mDOFs[DOF]
+      length(mDOFs) == 1 ? _DOF_to_dof(mDOFs[1],n_fmdofs) : zero(eltype(mDOFs))
+    end
+  end
+  for (cell,ldof_to_mdof) in enumerate(cell_ldof_to_mdof)
+    if cell_to_cellin[cell] != cell
+     empty!(ldof_to_mdof)
+    end
   end
-  cell_to_ldofs
+  cell_ldof_to_mdof
 end
 
 function _local_aggregates(cell_to_gcellin,gids::PRange)
   map(_local_aggregates,cell_to_gcellin,global_to_local(gids))
 end
 
 function _local_aggregates(cell_to_gcellin,gcell_to_cell)
+  T = eltype(cell_to_gcellin)
   map(cell_to_gcellin) do gcin
-    iszero(gcin) ? gcin : gcell_to_cell[ gcin ]
+    iszero(gcin) ? gcin : T(gcell_to_cell[ gcin ])
   end
 end
 

test/DistributedTests/PoissonTests.jl

@@ -103,7 +103,7 @@ function main(distribute,parts;
       "gid"=>own_to_global(gids)])#,
   #  cellfields=["uh"=>uh])
 
-  writevtk(Ω,"trian_O",cellfields=["uh"=>uh])
+  writevtk(Ω,"trian_O",cellfields=["uh"=>uh,"eh"=>e])
   writevtk(Γ,"trian_G")
   @test el2/ul2 < 1.e-8
   @test eh1/uh1 < 1.e-7

test/DistributedTests/sequential/PoissonTests.jl

@@ -3,6 +3,8 @@ using PartitionedArrays
 include("../PoissonTests.jl")
 with_debug() do distribute
     PoissonTests.main(distribute,(2,2))
+    PoissonTests.main(distribute,(1,4),cells=(4,8))
+    PoissonTests.main(distribute,(2,4),cells=(8,8))
     PoissonTests.main(distribute,(4,1),cells=(12,12),geometry=:remotes)
 end
 end