cleanup

Author: pbrubeck

firedrake/cython/dmcommon.pyx

@@ -4327,37 +4327,22 @@ def get_dm_cell_types(PETSc.DM dm):
     )
 
 
-def create_label_intersection(PETSc.DM dm, label_name, label_values):
-    """Return the intersection of the closure of subdomains of a DMPlex.
+def intersectIS(PETSc.IS i1, PETSc.IS i2):
+    """Return the intersection of two IS objects.
 
     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
+    i1 : PETSc.IS
+        The first IS.
+    i2 : PETSc.IS
+        The second IS.
 
     Returns
     -------
     PETSc.IS
-        A PETSc.IS with the points in the intersection.
+        A PETSc.IS with the intersection.
 
     """
-    cdef:
-        PETSc.IS iout, i1, i2
-        PETSc.DMLabel label
-
-    if len(label_values) == 0:
-        return PETSc.IS().createGeneral([], comm=dm.comm)
-
-    label = dm.getLabel(label_name)
-    CHKERR(DMPlexLabelComplete(dm.dm, label.dmlabel))
-    iout = label.getStratumIS(label_values[0])
-    for val in label_values[1:]:
-        i1 = iout
-        i2 = label.getStratumIS(val)
-        iout = PETSc.IS()
-        CHKERR(ISIntersect(i1.iset, i2.iset, &iout.iset))
+    cdef PETSc.IS iout = PETSc.IS()
+    CHKERR(ISIntersect(i1.iset, i2.iset, &iout.iset))
     return iout

firedrake/mesh.py

@@ -25,7 +25,7 @@
 from pyop2.mpi import (
     MPI, COMM_WORLD, temp_internal_comm
 )
-from functools import cached_property
+from functools import cached_property, reduce
 from pyop2.utils import as_tuple
 import petsctools
 from petsctools import OptionsManager, get_external_packages
@@ -4954,14 +4954,19 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
     if isinstance(subdomain_id, str):
         raise ValueError(f"Submesh construction got invalid subdomain_id {subdomain_id}.")
     elif isinstance(subdomain_id, Sequence):
+        label = plex.getLabel(label_name)
+        if subdim != dim:
+            plex.labelComplete(label)
         # Take the union of the labels in the list
         iset = PETSc.IS().createGeneral([], comm=mesh.comm)
         for sub in subdomain_id:
             if isinstance(sub, Sequence):
                 # Take the intersection of the labels from nested lists
-                cur = dmcommon.create_label_intersection(plex, label_name, sub)
+                if len(sub) == 0:
+                    continue
+                cur = reduce(dmcommon.intersectIS, map(label.getStratumIS, sub))
             else:
-                cur = plex.getStratumIS(label_name, sub)
+                cur = label.getStratumIS(sub)
             iset = iset.union(cur)
         # Create a temporary label
         label_name = "temp_label"