review comments

Author: pbrubeck

firedrake/cython/dmcommon.pyx

@@ -2203,11 +2203,11 @@ def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray):
 
 def _get_periodicity(dm: PETSc.DM) -> tuple[tuple[bool, bool], ...]:
     """Return mesh periodicity information.
-    
+
     This function returns a 2-tuple of bools per dimension where the first entry indicates
     whether the mesh is periodic in that dimension, and the second indicates whether the
     mesh is single-cell periodic in that dimension.
-    
+
     """
     cdef:
         const PetscReal *maxCell, *L
@@ -4325,3 +4325,41 @@ def get_dm_cell_types(PETSc.DM dm):
     return tuple(
         polytope_type_enum for polytope_type_enum, found in enumerate(found_all) if found
     )
+
+
+def create_label_intersection(PETSc.DM dm, label_name, label_values):
+    """Return the intersection of the closure of a subdomains of a DMPlex.
+
+    Parameters
+    ----------
+    dm : PETSc.DM
+        The DMPlex.
+    label_name : str
+        The name of the label
+    label_values : Sequence[int]
+        The values of the subdomain label to intersect
+
+    Returns
+    -------
+    tuple
+        A PETSc.IS with the points in the intersection.
+
+    """
+    cdef:
+        PETSc.DMLabel label
+        PETSc.PetscIS is1, is2
+        PetscInt val = label_values[0]
+
+    label = dm.getLabel(label_name)
+    CHKERR(DMPlexLabelComplete(dm.dm, label.dmlabel))
+    CHKERR(DMLabelGetStratumIS(<DMLabel>label.dmlabel, val, &is1))
+
+    for i in range(1, len(label_values)):
+        iout = PETSc.IS()
+        val = label_values[i]
+        CHKERR(DMLabelGetStratumIS(<DMLabel>label.dmlabel, val, &is2))
+        CHKERR(ISIntersect(is1, is2, &(<PETSc.IS?>iout).iset))
+        CHKERR(ISDestroy(&is1))
+        CHKERR(ISDestroy(&is2))
+        is1 = (<PETSc.IS?>iout).iset
+    return iout

firedrake/cython/petschdr.pxi

@@ -142,6 +142,7 @@ cdef extern from "petscis.h" nogil:
     PetscErrorCode ISLocalToGlobalMappingGetBlockIndices(PETSc.PetscLGMap, const PetscInt**)
     PetscErrorCode ISLocalToGlobalMappingRestoreBlockIndices(PETSc.PetscLGMap, const PetscInt**)
     PetscErrorCode ISDestroy(PETSc.PetscIS*)
+    PetscErrorCode ISIntersect(PETSc.PetscIS, PETSc.PetscIS, PETSc.PetscIS*)
 
 cdef extern from "petscsf.h" nogil:
     struct PetscSFNode_:

firedrake/mesh.py

@@ -4863,6 +4863,8 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
         raise NotImplementedError("Can not create a submesh of a ``VertexOnlyMesh``")
     if subdim is None:
         subdim = mesh.topological_dimension
+        if subdomain_id == "on_boundary":
+            subdim = subdim - 1
     plex = mesh.topology_dm
     dim = plex.getDimension()
     if subdim not in {dim, dim - 1}:
@@ -4880,20 +4882,24 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
             label_name = dmcommon.FACE_SETS_LABEL
 
     # Parse non-integer subdomain_id
-    if subdomain_id == "on_boundary":
-        subdomain_id = tuple(mesh.exterior_facets.unique_markers)
+    if isinstance(subdomain_id, str):
+        if subdomain_id == "on_boundary":
+            subdomain_id = tuple(mesh.exterior_facets.unique_markers)
+        else:
+            raise ValueError(f"Submesh construction got invalid subdomain_id {subdomain_id}.")
 
     if isinstance(subdomain_id, Sequence):
         # Create a temporary DMLabel with the union of the labels in the list
         icomm = comm or mesh.comm
         iset = PETSc.IS().createGeneral([], comm=icomm)
         for sub in subdomain_id:
+            try:
+                sub, = sub
+            except ValueError:
+                pass
             if isinstance(sub, Sequence):
                 # Take the intersection of the (closure of the) labels from nested lists
-                ises = [plex.getStratumIS(label_name, subi) for subi in sub]
-                closure = [[plex.getTransitiveClosure(p)[0] for p in i.indices] for i in ises]
-                indices = reduce(np.intersect1d, closure)
-                cur = PETSc.IS().createGeneral(indices, comm=icomm)
+                cur = dmcommon.create_label_intersection(plex, label_name, sub)
             else:
                 cur = plex.getStratumIS(label_name, sub)
             iset = iset.union(cur)

tests/firedrake/submesh/test_submesh_interface.py

@@ -3,7 +3,7 @@
 from firedrake import *
 
 
-def test_submesh_subdomain_id_tuple():
+def test_submesh_subdomain_id_union():
     mesh = UnitSquareMesh(4, 4)
     x, y = SpatialCoordinate(mesh)
     M = FunctionSpace(mesh, "DG", 0)
@@ -25,7 +25,7 @@ def test_submesh_subdomain_id_tuple():
     assert np.allclose(submesh2.coordinates.dat.data, submesh1.coordinates.dat.data)
 
 
-def test_submesh_subdomain_id_nested_tuple():
+def test_submesh_subdomain_id_intersection():
     mesh = UnitSquareMesh(4, 4)
     x, y = SpatialCoordinate(mesh)
     M = FunctionSpace(mesh, "DG", 0)
@@ -48,7 +48,7 @@ def test_submesh_subdomain_id_nested_tuple():
 
 
 @pytest.mark.parametrize("subdomain_id", ["on_boundary", (1, 3, 6)])
-def test_submesh_facet_subdomain_id_tuple(subdomain_id):
+def test_submesh_facet_subdomain_id_union(subdomain_id):
     mesh = UnitCubeMesh(2, 2, 2)
     submesh1 = Submesh(mesh, mesh.topological_dimension - 1, subdomain_id=subdomain_id)
     if subdomain_id == "on_boundary":
@@ -67,7 +67,7 @@ def test_submesh_facet_subdomain_id_tuple(subdomain_id):
     assert np.allclose(submesh2.coordinates.dat.data, submesh1.coordinates.dat.data)
 
 
-def test_submesh_facet_subdomain_id_nested_tuple():
+def test_submesh_facet_subdomain_id_intersection():
     mesh = UnitSquareMesh(4, 4)
     x, y = SpatialCoordinate(mesh)
     M = FunctionSpace(mesh, "DG", 0)