add a bunch of examples

Author: pbrubeck

firedrake/mesh.py

@@ -25,7 +25,7 @@
 from pyop2.mpi import (
     MPI, COMM_WORLD, temp_internal_comm
 )
-from functools import cached_property, reduce
+from functools import cached_property
 from pyop2.utils import as_tuple
 import petsctools
 from petsctools import OptionsManager, get_external_packages
@@ -4854,17 +4854,86 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
     ridges to be contained in the quad mesh are shared by at most two
     facets to make the quad mesh orientation algorithm work.
 
+    Examples
+    --------
+
+    Mark a cell subdomain and construct a codim-0 submesh from all cells in the subdomain
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> x, y = SpatialCoordinate(mesh)
+    >>> DG = FunctionSpace(mesh, "DG", 0)
+    >>> cell_marker = assemble(interpolate(conditional(lt(x, 0.5), 1, 0), DG))
+    >>> mesh.mark_entities(cell_marker, 111)
+    >>> submesh = Submesh(mesh, subdomain_id=111)
+
+    Mark a facet subdomain and construct a codim-1 submesh from all facets in the subdomain
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> x, y = SpatialCoordinate(mesh)
+    >>> DGT = FunctionSpace(mesh, "DGT", 0)
+    >>> facet_marker = assemble(interpolate(conditional(lt(abs(x-0.5), 1E-12), 1, 0), DGT))
+    >>> mesh.mark_entities(facet_marker, 222)
+    >>> submesh = Submesh(mesh, dim=mesh.topological_dimension-1, subdomain_id=222)
+
+    Construct a codim-0 submesh of the union of multiple subdomains by passing a list
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> x, y = SpatialCoordinate(mesh)
+    >>> DG = FunctionSpace(mesh, "DG", 0)
+    >>> mesh.mark_entities(assemble(interpolate(conditional(lt(x, 0.5), 1, 0), DG)), 1)
+    >>> mesh.mark_entities(assemble(interpolate(conditional(lt(y, 0.5), 1, 0), DG)), 2)
+    >>> submesh = Submesh(mesh, subdomain_id=[1, 2])
+
+    Construct a codim-0 submesh of the intersection of multiple subdomains by passing a nested list
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> x, y = SpatialCoordinate(mesh)
+    >>> DG = FunctionSpace(mesh, "DG", 0)
+    >>> mesh.mark_entities(assemble(interpolate(conditional(lt(x, 0.5), 1, 0), DG)), 1)
+    >>> mesh.mark_entities(assemble(interpolate(conditional(lt(y, 0.5), 1, 0), DG)), 2)
+    >>> submesh = Submesh(mesh, subdomain_id=[(1, 2)])
+
+    Construct a codim-1 submesh of all the facets (the skeleton mesh)
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> submesh = Submesh(mesh, subdim=1)
+
+    Construct a codim-1 submesh of the entire boundary
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> submesh = Submesh(mesh, subdomain_id="on_boundary")
+
+    Construct a codim-1 submesh of the union of multiple boundaries
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> submesh = Submesh(mesh, subdim=mesh.topological_dimension-1, subdomain_id=[1, 2, 3])
+
+    Construct a codim-0 submesh of the part of the mesh owned by each MPI rank
+
+    >>> mesh = UnitSquareMesh(4, 4)
+    >>> submesh = Submesh(mesh, comm=COMM_SELF)
+
     """
     if not isinstance(mesh, MeshGeometry):
         raise TypeError("Parent mesh must be a `MeshGeometry`")
     if isinstance(mesh.topology, ExtrudedMeshTopology):
         raise NotImplementedError("Can not create a submesh of an ``ExtrudedMesh``")
     elif isinstance(mesh.topology, VertexOnlyMeshTopology):
         raise NotImplementedError("Can not create a submesh of a ``VertexOnlyMesh``")
+
+    if subdomain_id == "on_boundary":
+        if subdim is None:
+            subdim = mesh.topological_dimension - 1
+        elif subdim != mesh.topological_dimension - 1:
+            raise ValueError('subdomain_id="on_boundary" requires subdim=dim-1')
+        if label_name is None:
+            label_name = "exterior_facets"
+        elif label_name != "exterior_facets":
+            raise ValueError('subdomain_id="on_boundary" requires label_name="exterior_facets"')
+        subdomain_id = 1
+
     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}:
@@ -4883,13 +4952,9 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
 
     # Parse non-integer subdomain_id
     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
+        raise ValueError(f"Submesh construction got invalid subdomain_id {subdomain_id}.")
+    elif isinstance(subdomain_id, Sequence):
+        # 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):
@@ -4898,6 +4963,7 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
             else:
                 cur = plex.getStratumIS(label_name, sub)
             iset = iset.union(cur)
+        # Create a temporary label
         label_name = "temp_label"
         subdomain_id = 1
         plex.createLabel(label_name)