Fixes

* Do not try and distribute when serial
* Use a slightly different approach to building a topology_dm (to
  preserve labels)

Author: connorjward

firedrake/cython/dmcommon.pyx

@@ -3972,6 +3972,7 @@ def dmplex_migrate(PETSc.DM dm, PETSc.SF sf) -> PETSc.DM:
     return migrated_dm
 
 
+# not used currently
 def densify_sf(PETSc.DM topology_dm, PETSc.SF sparse_sf) -> PETSc.SF:
     cdef:
         PETSc.SF dense_sf
@@ -4008,10 +4009,16 @@ def densify_sf(PETSc.DM topology_dm, PETSc.SF sparse_sf) -> PETSc.SF:
     return dense_sf
 
 
+# not used currently
 def dmplex_create_overlap_migration_sf(PETSc.DM topology_dm, PETSc.SF overlap_sf) -> PETSc.SF:
     cdef:
         PETSc.SF migration_sf
 
     migration_sf = PETSc.SF().create(comm=topology_dm.comm)
     CHKERR(DMPlexCreateOverlapMigrationSF(topology_dm.dm, overlap_sf.sf, &migration_sf.sf))
     return migration_sf
+
+
+# not used currently
+def set_label(PETSc.DM dm, PETSc.DMLabel label) -> None:
+    CHKERR(DMSetLabel(dm.dm, label.dmlabel))

firedrake/cython/petschdr.pxi

@@ -77,6 +77,7 @@ cdef extern from "petscdmlabel.h" nogil:
 cdef extern from "petscdm.h" nogil:
     int DMCreateLabel(PETSc.PetscDM,char[])
     int DMGetLabel(PETSc.PetscDM,char[],DMLabel*)
+    int DMSetLabel(PETSc.PetscDM,PETSc.PetscDMLabel)
     int DMGetPointSF(PETSc.PetscDM,PETSc.PetscSF*)
     int DMSetLabelValue(PETSc.PetscDM,char[],PetscInt,PetscInt)
     int DMGetLabelValue(PETSc.PetscDM,char[],PetscInt,PetscInt*)

firedrake/mesh.py

@@ -2318,16 +2318,11 @@ def callback(self):
             del self._callback
             # Finish the initialisation of mesh topology
             self.topology.init()
-            # geometry_dm = self._input_geometry_dm.migrate(self.topology.overlap_sf)
-
-            # make the overlap SF dense
-            # dense_overlap_sf = dmcommon.densify_sf(self.topology.topology_dm, self.topology.overlap_sf)
-            # breakpoint()
-            # migration_sf = dmcommon.dmplex_create_overlap_migration_sf(self.topology.topology_dm, self.topology.overlap_sf)
-            # geometry_dm = dmcommon.dmplex_migrate(self._input_geometry_dm, dense_overlap_sf)
-            # geometry_dm = dmcommon.dmplex_migrate(self._input_geometry_dm, self.topology.overlap_sf)
-            # geometry_dm = dmcommon.dmplex_migrate(self._input_geometry_dm, migration_sf)
-            geometry_dm = dmcommon.dmplex_migrate(self._input_geometry_dm, self.topology.sfBC)
+            if self.comm.size > 1:
+                assert self.topology.sfBC is not None
+                geometry_dm = dmcommon.dmplex_migrate(self._input_geometry_dm, self.topology.sfBC)
+            else:
+                geometry_dm = self._input_geometry_dm
 
             coordinates_fs = functionspace.FunctionSpace(self.topology, self.ufl_coordinate_element())
 
@@ -3128,22 +3123,32 @@ def Mesh(meshfile, **kwargs):
                                % (meshfile, ext[1:]))
         plex.setName(_generate_default_mesh_topology_name(name))
 
+    # TODO: Push other labels onto the coordinate DM, this should be done
+    # by DMPlex
+    # NOTE: OR we instead set the plex as the coordinate DM by dropping
+    # coordinates
+    # topology_dm = plex.getCoordinateDM()
+    #
+    # face_sets_label = plex.getLabel(dmcommon.FACE_SETS_LABEL)
+    # topology_dm.createLabel(dmcommon.FACE_SETS_LABEL)
+    # dmcommon.set_label(topology_dm, face_sets_label)
+    topology_dm = plex.clone()
+    topology_dm.setCoordinates(plex.getCoordinateDM().getCoordinates())
+
     # Create mesh topology
     # Pass the coordinate DM because we only want the topology here
-    topology = MeshTopology(plex.getCoordinateDM(), name=plex.getName(), reorder=reorder,
+    # topology = MeshTopology(plex.getCoordinateDM(), name=plex.getName(), reorder=reorder,
+    topology = MeshTopology(topology_dm, name=plex.getName(), reorder=reorder,
                             distribution_parameters=distribution_parameters,
                             distribution_name=kwargs.get("distribution_name"),
                             permutation_name=kwargs.get("permutation_name"),
                             comm=user_comm)
 
-    # distributed_plex = dmcommon.dmplex_migrate(plex, topology.sfBC)
-    distributed_plex = plex
-
     if netgen and isinstance(meshfile, netgen.libngpy._meshing.Mesh):
         netgen_firedrake_mesh.createFromTopology(topology, name=name, comm=user_comm)
         mesh = netgen_firedrake_mesh.firedrakeMesh
     else:
-        mesh = make_mesh_from_mesh_topology(topology, distributed_plex, name)
+        mesh = make_mesh_from_mesh_topology(topology, plex, name)
     mesh._tolerance = tolerance
     return mesh
 

firedrake/utility_meshes.py

@@ -736,13 +736,12 @@ def TensorRectangleMesh(
     )
 
     # mark boundary facets
-    topology_dm = plex.getCoordinateDM()
-    topology_dm.createLabel(dmcommon.FACE_SETS_LABEL)
-    topology_dm.markBoundaryFaces("boundary_faces")
+    plex.createLabel(dmcommon.FACE_SETS_LABEL)
+    plex.markBoundaryFaces("boundary_faces")
     coords = plex.getCoordinates()
     coord_sec = plex.getCoordinateSection()
-    if topology_dm.getStratumSize("boundary_faces", 1) > 0:
-        boundary_faces = topology_dm.getStratumIS("boundary_faces", 1).getIndices()
+    if plex.getStratumSize("boundary_faces", 1) > 0:
+        boundary_faces = plex.getStratumIS("boundary_faces", 1).getIndices()
         xtol = 0.5 * min(xcoords[1] - xcoords[0], xcoords[-1] - xcoords[-2])
         ytol = 0.5 * min(ycoords[1] - ycoords[0], ycoords[-1] - ycoords[-2])
         x0 = xcoords[0]
@@ -752,14 +751,14 @@ def TensorRectangleMesh(
         for face in boundary_faces:
             face_coords = plex.vecGetClosure(coord_sec, coords, face)
             if abs(face_coords[0] - x0) < xtol and abs(face_coords[2] - x0) < xtol:
-                topology_dm.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 1)
+                plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 1)
             if abs(face_coords[0] - x1) < xtol and abs(face_coords[2] - x1) < xtol:
-                topology_dm.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 2)
+                plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 2)
             if abs(face_coords[1] - y0) < ytol and abs(face_coords[3] - y0) < ytol:
-                topology_dm.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 3)
+                plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 3)
             if abs(face_coords[1] - y1) < ytol and abs(face_coords[3] - y1) < ytol:
-                topology_dm.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 4)
-    topology_dm.removeLabel("boundary_faces")
+                plex.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 4)
+    plex.removeLabel("boundary_faces")
     m = mesh.Mesh(
         plex,
         reorder=reorder,