@@ -733,3 +733,201 @@ function refine_cell_gids(
733733 end
734734 return refine_cell_gids (cmodel,fmodels,f_own_to_local)
735735end
736+
737+ # Coarsening for polytopal meshes
738+ #
739+ # We assume some properties of the coarsening:
740+ # - The patch topology must the a partition of the owned fine cells, i.e
741+ # + All patches are disjoint
742+ # + All owned fine cells are part of a patch
743+ # - All patches are fully owned by a single processor. I.e we do NOT allow
744+ # patches that are split between multiple processors.
745+ # The second property is slightly restrictive, but I think it is very reasonable... allowing
746+ # for split patches would add an insane amount of complexity to the code. Basically,
747+ # this property means that we can always coarsen the owned fine cells locally, then
748+ # communicate to find out ghosts.
749+ # If we ever need to have split patches, one should instead re-distribute the model accordingly
750+ # and then coarsen.
751+ #
752+ # The main idea of the algorithm is as follows:
753+ # 1. First, we coarsen the owned part of the models
754+ # 2. We then create a global numering of the coarse cells
755+ # 3. We can use the above to build the cell-to-node connectivity
756+ # 4. We also need a global numbering of the coarse nodes, which we need to
757+ # communicate the model coordinates
758+ # 5. We can then create our coarse model
759+ function Adaptivity. coarsen (fmodel:: DistributedDiscreteModel , ptopo:: DistributedPatchTopology ; return_glue= false )
760+
761+ # First, some preliminary checks
762+ fgids = partition (get_cell_gids (fmodel))
763+ map (local_views (ptopo), fgids) do ptopo, fids
764+ patch_cells = Geometry. get_patch_cells (ptopo)
765+ lcell_to_ocell = local_to_own (fids)
766+ ocell_to_lcell = own_to_local (fids)
767+ @check allunique (patch_cells. data) " Patches must not overlap"
768+ @check all (c -> ! iszero (lcell_to_ocell[c]), patch_cells. data) " All local patches must be fully owned"
769+ @check issubset (ocell_to_lcell, patch_cells. data) " All owned cells must be part of a patch"
770+ end
771+
772+ # 1. Local coarsening on the owned part of the model
773+ own_polys, own_connectivity = map (local_views (fmodel), local_views (ptopo)) do fmodel, ptopo
774+ Adaptivity. generate_patch_polytopes (fmodel,ptopo)
775+ end |> tuple_of_arrays
776+
777+ # 2. Coarse cell gids
778+ # - Each processor can easily create the global numbering of their owned coarse cells.
779+ # - Through the fine cell ghosts, we have to communicate the ghost coarse cell gids, and
780+ # from that get the local-to-global mapping.
781+ Dc = num_cell_dims (fmodel)
782+ fcell_gids = partition (get_face_gids (fmodel, Dc))
783+ n_own_ccells = map (Geometry. num_patches, local_views (ptopo))
784+ n_cgids = sum (n_own_ccells)
785+ first_cgid = scan (+ ,n_own_ccells,type= :exclusive ,init= 1 )
786+ fcell_to_cgid = map (local_views (ptopo), fcell_gids, first_cgid) do ptopo, fcell_gids, first_cgid
787+ patch_cells = Geometry. get_patch_cells (ptopo)
788+ fcell_to_cgid = zeros (Int, local_length (fcell_gids))
789+ Arrays. flatten_partition! (fcell_to_cgid,patch_cells)
790+ fcell_to_cgid .+ = first_cgid - 1
791+ return fcell_to_cgid
792+ end
793+ consistent! (PVector (fcell_to_cgid, fcell_gids)) |> wait
794+
795+ ccell_gids, fcell_to_ccell = map (fcell_gids, fcell_to_cgid, first_cgid, n_own_ccells) do fids, fcell_to_cgid, first_cgid, n_own_ccells
796+ owner = part_id (fids)
797+ own_range = first_cgid: (first_cgid+ n_own_ccells- 1 )
798+ c_own_to_global = collect (Int, own_range)
799+ c_own_to_flid = collect (Int32, indexin (c_own_to_global, fcell_to_cgid))
800+
801+ is_ghost (gid) = ! iszero (gid) && (gid ∉ own_range)
802+ c_ghost_to_global = filter! (is_ghost, unique (fcell_to_cgid))
803+ c_ghost_to_flid = collect (Int32, indexin (c_ghost_to_global, fcell_to_cgid))
804+ c_ghost_to_owner = local_to_owner (fids)[c_ghost_to_flid]
805+
806+ own_cgids = OwnIndices (n_cgids, owner, c_own_to_global)
807+ ghost_cgids = GhostIndices (n_cgids, c_ghost_to_global, c_ghost_to_owner)
808+ cgids = OwnAndGhostIndices (own_cgids, ghost_cgids)
809+
810+ cgid_to_clid = global_to_local (cgids)
811+ fcell_to_ccell = zeros (Int32, local_length (fids))
812+ for (flid, cgid) in enumerate (fcell_to_cgid)
813+ fcell_to_ccell[flid] = cgid_to_clid[cgid]
814+ end
815+ return cgids, fcell_to_ccell
816+ end |> tuple_of_arrays
817+
818+ # 3 & 4. Coarse node gids:
819+ # - Each coarse node corresponds to a fine node. However: not all local coarse nodes
820+ # are also local fine nodes. Therefore we will have to work with the
821+ # global fine node ids to have a unique numbering.
822+ # - We will communicate the coarse-cell-to-fine-node-gid connectivity, then
823+ # build a coarse node numbering from that.
824+
825+ # First: We communicate the number of nodes per coarse cell
826+ ccell_to_nnodes = map (ccell_gids, own_connectivity) do ccids, own_connectivity
827+ nnodes = zeros (Int32, local_length (ccids))
828+ nnodes[own_to_local (ccids)] .= map (length, own_connectivity)
829+ return nnodes
830+ end
831+ consistent! (PVector (ccell_to_nnodes, ccell_gids)) |> wait
832+
833+ # We communicate the connectivity:
834+ # For each coarse cell, the fine node gids in that cell
835+ fnode_gids = partition (get_face_gids (fmodel, 0 ))
836+ connectivity = map (
837+ own_connectivity, ccell_to_nnodes, ccell_gids, fnode_gids
838+ ) do own_connectivity, ccell_to_nnodes, ccids, fnids
839+ fnode_local_to_global = local_to_global (fnids)
840+ ptrs = pushfirst! (ccell_to_nnodes, 0 )
841+ Arrays. length_to_ptrs! (ptrs)
842+ data = zeros (Int32, ptrs[end ]- 1 )
843+ for (oid, lid) in enumerate (own_to_local (ccids))
844+ node_lids = view (own_connectivity, oid)
845+ node_gids = view (fnode_local_to_global,node_lids)
846+ data[ptrs[lid]: (ptrs[lid+ 1 ]- 1 )] .= node_gids
847+ end
848+ return JaggedArray (data, ptrs)
849+ end
850+ consistent! (PVector (connectivity, ccell_gids)) |> wait
851+ connectivity = map (c -> Table (c. data, c. ptrs), connectivity)
852+
853+ # We create a local numbering of the coarse nodes and renumber the connectivity.
854+ # For later, we also return the local-to-local mapping of the fine nodes to coarse nodes.
855+ n_cnodes, fnode_to_cnode = map (fnode_gids, connectivity) do fnids, conn
856+ n_lid = 0
857+ fgid_to_clid = Dict {Int,Int32} ()
858+ for k in eachindex (conn. data)
859+ fgid = conn. data[k]
860+ clid = get! (fgid_to_clid, fgid, n_lid + 1 )
861+ n_lid += (clid == n_lid + 1 ) # Increment if it's new
862+ conn. data[k] = clid # Renumber the connectivity
863+ end
864+ fnode_to_cnode = zeros (Int32, local_length (fnids))
865+ for (flid, fgid) in enumerate (local_to_global (fnids))
866+ fnode_to_cnode[flid] = get! (fgid_to_clid,fgid,0 )
867+ end
868+ return n_lid, fnode_to_cnode
869+ end |> tuple_of_arrays
870+
871+ # Finally, we use pre-existing routines to generate the coarse node gids
872+ cnode_gids = generate_gids (PRange (ccell_gids), connectivity, n_cnodes) |> partition
873+
874+ # Coarse node coordinates:
875+ cnode_coords = map (local_views (fmodel), cnode_gids, fnode_to_cnode) do fmodel, cngids, fnode_to_cnode
876+ @check issubset (own_to_local (cngids), fnode_to_cnode)
877+ fnode_coords = Geometry. get_vertex_coordinates (get_grid_topology (fmodel))
878+ cnode_coords = Vector {eltype(fnode_coords)} (undef, local_length (cngids))
879+ for (flid, clid) in enumerate (fnode_to_cnode)
880+ if clid > 0
881+ cnode_coords[clid] = fnode_coords[flid] # Fill owned info
882+ end
883+ end
884+ return cnode_coords
885+ end
886+ consistent! (PVector (cnode_coords, cnode_gids)) |> wait # Communicate coarse coords
887+
888+ # Coarse polytopes:
889+ # - We already have the owned polytopes from the local coarsening
890+ # - We create the ghost polytopes from the connectivity and the coarse node coordinates
891+ # TODO : The creation of 3D polyhedra requires a bit more information! This is the only thing missing for 3D.
892+ @notimplementedif Dc != 2 " Coarsening is only implemented for 2D polytopal meshes"
893+ polys = map (ccell_gids, connectivity, own_polys, cnode_coords) do ccids, conn, own_polys, cnode_coords
894+ polys = Vector {eltype(own_polys)} (undef, local_length (ccids))
895+ for (lid, oid) in enumerate (local_to_own (ccids))
896+ if oid > 0
897+ polys[lid] = own_polys[oid] # Owned polytope
898+ else
899+ polys[lid] = Polygon (cnode_coords[view (conn, lid)]) # Ghost polytope
900+ end
901+ end
902+ return polys
903+ end
904+
905+ # 5. We can finally create our coarse model
906+ face_gids = Vector {PRange} (undef, Dc+ 1 )
907+ face_gids[end ] = PRange (ccell_gids)
908+ face_gids[1 ] = PRange (cnode_gids)
909+ ctopo = GridapDistributed. DistributedGridTopology (
910+ map (Geometry. PolytopalGridTopology, cnode_coords, connectivity, polys), face_gids
911+ )
912+ _setup_consistent_faces! (ctopo)
913+ labels = Geometry. FaceLabeling (ctopo)
914+ cmodels = map (local_views (ctopo), local_views (labels)) do ctopo, labels
915+ cgrid = Geometry. PolytopalGrid (ctopo)
916+ Geometry. PolytopalDiscreteModel (cgrid, ctopo, labels)
917+ end
918+ cmodel = GridapDistributed. DistributedDiscreteModel (cmodels, face_gids)
919+ (! return_glue) && (return cmodel)
920+
921+ # If required, create the adaptivity glues
922+ ftopo = get_grid_topology (fmodel)
923+ glues = map (
924+ ccell_gids, cnode_gids, fcell_to_ccell, fnode_to_cnode, local_views (ftopo), local_views (ctopo)
925+ ) do ccids, cnids, fcell_to_ccell, fnode_to_cnode, ftopo, ctopo
926+ ccell_to_fcell = Arrays. inverse_table (fcell_to_ccell, local_length (ccids))
927+ cnode_to_fnode = find_inverse_index_map (fnode_to_cnode, local_length (cnids))
928+ Adaptivity. generate_patch_adaptivity_glue (
929+ ftopo, ctopo, fcell_to_ccell, ccell_to_fcell, fnode_to_cnode, cnode_to_fnode,
930+ )
931+ end
932+ return cmodel, glues
933+ end
0 commit comments