Cell aggregation on p4est meshes

Author: ericneiva

src/AgFEM/CellAggregation.jl

@@ -89,11 +89,18 @@ function aggregate(
   cut::EmbeddedDiscretization,
   geo::CSG.Geometry,
   lid_to_gid::AbstractVector,
-  in_or_out)
+  in_or_out;
+  grid_topology=nothing)
 
   facet_to_inoutcut = compute_bgfacet_to_inoutcut(cut.bgmodel,geo)
   cell_to_lcellin, cell_to_value, _ =
-    _aggregate_by_threshold(strategy.threshold,cut,geo,in_or_out,facet_to_inoutcut,lid_to_gid)
+    _aggregate_by_threshold(strategy.threshold,
+                            cut,
+                            geo,
+                            in_or_out,
+                            facet_to_inoutcut,
+                            lid_to_gid,
+                            grid_topology)
   cell_to_gcellin = fill(0,length(lid_to_gid))
   for i in eachindex(cell_to_gcellin)
     if !iszero(cell_to_lcellin[i])
@@ -132,7 +139,13 @@ function aggregate(
   cell_to_cellin
 end
 
-function _aggregate_by_threshold(threshold,cut,geo,loc,facet_to_inoutcut,lid_to_gid)
+function _aggregate_by_threshold(threshold,
+                                 cut,
+                                 geo,
+                                 loc,
+                                 facet_to_inoutcut,
+                                 lid_to_gid,
+                                 topo=nothing)
   @assert loc in (IN,OUT)
 
   cutinorout = loc == IN ? (CUT_IN,IN) : (CUT_OUT,OUT)
@@ -146,7 +159,8 @@ function _aggregate_by_threshold(threshold,cut,geo,loc,facet_to_inoutcut,lid_to_
   cell_to_inoutcut = compute_bgcell_to_inoutcut(cut,geo)
 
   cell_to_coords = get_cell_coordinates(bgtrian)
-  topo = get_grid_topology(model)
+  topo === nothing ? topo = get_grid_topology(model) : topo
+
   D = num_cell_dims(model)
   cell_to_faces = get_faces(topo,D,D-1)
   face_to_cells = get_faces(topo,D-1,D)

test/dev/eric_dev.jl

@@ -2,6 +2,7 @@ module DistributedAggregationP4estMeshes
 
   using Gridap
   using GridapEmbedded
+  using GridapEmbedded.Interfaces: CUT
   using GridapDistributed
   using PartitionedArrays
   using MPI
@@ -27,16 +28,72 @@ module DistributedAggregationP4estMeshes
                                             num_uniform_refinements;
                                             num_ghost_layers=num_ghost_layers)
 
+    geo1 = quadrilateral(;x0=Point(-0.9,-0.9),
+                         d1=VectorValue(1.8,0.0),
+                         d2=VectorValue(0.0,1.8))
+    geo2 = ! disk(0.4)
+    geo = intersect(geo1,geo2)
+
+    cutgeo = cut(dmodel, geo)
+    cell_to_inoutcut = compute_bgcell_to_inoutcut(cutgeo,geo)
     fmodel_refine_coarsen_flags = 
-      map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices
+      map(ranks,
+          partition(get_cell_gids(dmodel.dmodel)),
+          cell_to_inoutcut) do rank,indices,cell_to_inoutcut
         flags = zeros(Int,length(indices))
-        flags .= nothing_flag        
-        flags[1] = refine_flag
+        flags .= nothing_flag
+        toref = findall(c->c==CUT,cell_to_inoutcut)
+        flags[toref] .= refine_flag
         flags
     end
     fmodel,_ = Gridap.Adaptivity.adapt(dmodel,fmodel_refine_coarsen_flags);
 
+    for i in 1:4
+      cutgeo = cut(fmodel, geo)
+      cell_to_inoutcut = compute_bgcell_to_inoutcut(cutgeo,geo)
+      fmodel_refine_coarsen_flags = 
+        map(ranks,
+            partition(get_cell_gids(fmodel)),
+            cell_to_inoutcut) do rank,indices,cell_to_inoutcut
+          flags = zeros(Int,length(indices))
+          flags .= nothing_flag
+          toref = findall(c->c==CUT,cell_to_inoutcut)
+          flags[toref] .= refine_flag
+          flags
+      end
+      fmodel,_ = Gridap.Adaptivity.adapt(fmodel,fmodel_refine_coarsen_flags);
+    end
+
+    writevtk(EmbeddedBoundary(cutgeo),"data/quad_bnd");
+    # The next line fails unless I comment out the assertion 
+    # in the constructor of AppendedTriangulations.jl:
+    # @assert get_background_model(a) === get_background_model(b)
+    writevtk(Triangulation(cutgeo,PHYSICAL_IN),"data/quad_phys");
+
+    cell_gids = get_cell_gids(fmodel)
+    cell_indices = partition(cell_gids)
+
+    cutgeo = cut(fmodel, geo)
     ncgt = NonConformingGridTopology(fmodel)
+    strategy = AggregateCutCellsByThreshold(1.0)
+    lcell_to_lroot, lcell_to_root, lcell_to_value =
+      map(local_views(cutgeo),cell_indices,ncgt) do cutgeo,cell_indices,ncgt
+        lid_to_gid = local_to_global(cell_indices)
+        aggregate(strategy,cutgeo,geo,lid_to_gid,IN,grid_topology=ncgt)
+      end |> tuple_of_arrays
+
+    # Output for verification of lcell_to_root map
+    ocell_to_root = map(lcell_to_root,own_to_local(cell_gids)) do agg,o_to_l
+      map(Reindex(agg),o_to_l)
+    end
+
+    writevtk(EmbeddedBoundary(cutgeo),"data/quad_bnd");
+    writevtk(Triangulation(cutgeo,PHYSICAL),"data/quad_phys");
+    writevtk(
+      Triangulation(fmodel), "data/quad_aggregates", 
+      celldata = ["aggregate" => ocell_to_root],
+    );
+
     map(ncgt) do ncgt
       get_faces(ncgt,D,0)
       get_faces(ncgt,0,D)