Submesh: support tuple subdomain_id

Author: pbrubeck

firedrake/mesh.py

@@ -4806,10 +4806,11 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
     subdim : int | None
         Topological dimension of the submesh.
         Defaults to ``mesh.topological_dimension``.
-    subdomain_id : int | None
+    subdomain_id : int | Sequence | None
         Subdomain ID representing the submesh.
-        If `None` the submesh will cover the entire domain.
-        This is useful to obtain a codim-1 submesh over all facets or
+        If multiple subdomain IDs are provided, their union is taken.
+        If `None` the submesh will cover the entire domain,
+        this is useful to obtain a codim-1 submesh over all facets or
         a submesh over a different communicator.
     label_name : str | None
         Name of the label to search ``subdomain_id`` in.
@@ -4876,8 +4877,23 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
             label_name = dmcommon.CELL_SETS_LABEL
         elif subdim == dim - 1:
             label_name = dmcommon.FACE_SETS_LABEL
+
+    if isinstance(subdomain_id, (tuple, list)):
+        # A list of subdomain ids requires us to build an internal DM label with the union
+        iset = PETSc.IS().createGeneral([], comm=comm or mesh.comm)
+        for sub in subdomain_id:
+            iset = iset.union(plex.getStratumIS(label_name, sub))
+        label_name = "temp_union"
+        subdomain_id = 1
+        plex.createLabel(label_name)
+        label = plex.getLabel(label_name)
+        label.setStratumIS(subdomain_id, iset)
+
     subplex = dmcommon.submesh_create(plex, subdim, label_name, subdomain_id, ignore_halo, comm=comm)
 
+    if label_name == "temp_union":
+        plex.removeLabel(label_name)
+
     comm = comm or mesh.comm
     name = name or _generate_default_submesh_name(mesh.name)
     subplex.setName(_generate_default_mesh_topology_name(name))

tests/firedrake/submesh/test_submesh_facet.py

@@ -134,3 +134,18 @@ def test_submesh_facet_all_facets():
     rmesh = RelabeledMesh(mesh, [facet_function], [facet_value])
     submesh2 = Submesh(rmesh, mesh.topological_dimension - 1, facet_value)
     assert submesh2.cell_set.size == submesh1.cell_set.size
+
+
+def test_submesh_facet_subdomain_id_tuple():
+    mesh = UnitCubeMesh(2, 2, 2)
+    subdomain_id = (1, 3, 6)
+    submesh1 = Submesh(mesh, mesh.topological_dimension - 1, subdomain_id=subdomain_id)
+    assert abs(assemble(1*dx(domain=submesh1)) - len(subdomain_id)) < 1E-12
+
+    V = FunctionSpace(mesh, "HDiv Trace", 0)
+    facet_function = Function(V)
+    DirichletBC(V, 1, subdomain_id).apply(facet_function)
+    facet_value = 999
+    rmesh = RelabeledMesh(mesh, [facet_function], [facet_value])
+    submesh2 = Submesh(rmesh, mesh.topological_dimension - 1, facet_value)
+    assert submesh2.cell_set.size == submesh1.cell_set.size