Support string subdomain_id in ufl.Measure, DirichletBC and Submesh

.github/actions/install/action.yml

@@ -152,6 +152,7 @@ runs:
           --extra-index-url https://download.pytorch.org/whl/cpu \
           "./firedrake-repo[${{ inputs.deps }}]"
 
+        pip install -v --no-deps --ignore-installed git+https://github.com/firedrakeproject/ufl.git@pbrubeck/string-subdomain-id
         firedrake-clean
         pip list
 

demos/boussinesq/boussinesq.py.rst

@@ -83,9 +83,6 @@ to prevent hydrostatic equilibrium and allow for a non-trivial velocity solution
 
     ngmesh = ngocc.OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.1)
 
-    left_id = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(dim=1)) if name == "left"]
-    right_id = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(dim=1)) if name == "right"]
-
     mesh = Mesh(ngmesh)
     x, y = SpatialCoordinate(mesh)
 
@@ -133,8 +130,8 @@ The nonlinear form for the problem is::
      - inner(div(u), q) * dx                            # Incompressibility constraint
      + inner(dot(u, grad(T)), w) * dx                   # Temperature advection term
      + inner(grad(T), grad(w)) * dx                     # Temperature diffusion term
-     - inner(g_left, w) * ds(tuple(left_id))            # Weakly imposed Neumann BC terms
-     - inner(g_right, w) * ds(tuple(right_id))          # Weakly imposed Neumann BC terms
+     - inner(g_left, w) * ds("left")                    # Weakly imposed Neumann BC terms
+     - inner(g_right, w) * ds("right")                  # Weakly imposed Neumann BC terms
      + inner(T - T0, s) * dx                            # Integral constraint on T
      )
 

demos/netgen/netgen_mesh.py.rst

@@ -89,7 +89,7 @@ Example: Poisson Problem
 -------------------------
 Let's now have a look at some features that can be useful when solving a variational problem using a Netgen mesh.
 We will consider as example the Poisson problem with constant source data and homogeneous boundary conditions on the boundary of the PacMan.
-The only method new to the reader should be the ``GetRegionNames`` which allows to find the IDs of the boundary we have labeled in Netgen. As usual we begin defining the ``FunctionSpace`` that will be used in our discretisation, defining trial and test functions and the source data::
+As usual we begin defining the ``FunctionSpace`` that will be used in our discretisation, defining trial and test functions and the source data::
 
    V = FunctionSpace(msh, "CG", 1)
    u = TrialFunction(V)
@@ -110,10 +110,9 @@ In code this becomes: ::
    a = inner(grad(u), grad(v))*dx
    L = inner(f, v) * dx
 
-Now we are ready to assemble the stiffness matrix for the problem. Since we want to enforce Dirichlet boundary conditions we construct a ``DirichletBC`` object and we use the ``GetRegionNames`` method from the Netgen mesh in order to map the label we have given when describing the geometry to the PETSc ``DMPLEX`` IDs. In particular if we look for the IDs of boundary element labeled either "line" or "curve" we would get::
+Now we are ready to assemble the stiffness matrix for the problem. Since we want to enforce Dirichlet boundary conditions we construct a ``DirichletBC`` object with the same string labels from the Netgen mesh. Here we specify homogenous boundary conditions on boundary elements labeled as ``"line"`` or ``"curve"``::
 
-   labels = [i+1 for i, name in enumerate(ngmsh.GetRegionNames(codim=1)) if name in ["line", "curve"]]
-   bc = DirichletBC(V, 0, labels)
+   bc = DirichletBC(V, 0, ["line", "curve"])
    print(labels)
 
 We then proceed to solve the problem::

firedrake/bcs.py

@@ -43,6 +43,9 @@ class BCBase(object):
     def __init__(self, V, sub_domain):
 
         self._function_space = V
+        mesh = V.mesh()
+        if hasattr(mesh, "parse_subdomain_id"):
+            sub_domain = mesh.parse_subdomain_id(mesh.topological_dimension - 1, sub_domain)
         self.sub_domain = (sub_domain, ) if isinstance(sub_domain, str) else as_tuple(sub_domain)
         # If this BC is defined on a subspace (IndexedFunctionSpace or
         # ComponentFunctionSpace, possibly recursively), pull out the appropriate

firedrake/functionspace.py

@@ -342,10 +342,12 @@ def RestrictedFunctionSpace(function_space, boundary_set=[], name=None):
     make_space = type(function_space)
     mesh = function_space.mesh()
 
+    if isinstance(boundary_set, str):
+        boundary_set = (boundary_set,)
     if function_space.boundary_set:
         boundary_set = frozenset(boundary_set) | function_space.boundary_set
         function_space = function_space.function_space
-
+    boundary_set = mesh.parse_subdomain_id(mesh.topological_dimension - 1, boundary_set)
     return make_space(impl.RestrictedFunctionSpace(function_space.topological,
                                                    boundary_set=boundary_set,
                                                    name=name),

firedrake/mesh.py

@@ -163,6 +163,74 @@ def _generate_default_mesh_topology_permutation_name(reorder):
     return "_".join(["firedrake", "default", str(reorder)])
 
 
+def _as_subdomain_id_list(subdomain_id):
+    if isinstance(subdomain_id, str) or not isinstance(subdomain_id, (list, tuple, set, frozenset, np.ndarray)):
+        subdomain_ids = (subdomain_id,)
+    else:
+        subdomain_ids = subdomain_id
+
+    ids = []
+    for sid in subdomain_ids:
+        if not isinstance(sid, numbers.Integral):
+            raise TypeError(f"Expected integer subdomain id, got {sid!r}")
+        ids.append(int(sid))
+    return ids
+
+
+def _parse_gmsh_physical_names(filename):
+    """Parse the ``$PhysicalNames`` block from a Gmsh ``.msh`` file."""
+    region_names = []
+    try:
+        with open(filename, encoding="utf-8", errors="replace") as f:
+            in_physical_names = False
+            for line in f:
+                line = line.strip()
+                if line == "$PhysicalNames":
+                    in_physical_names = True
+                    continue
+                if line == "$EndPhysicalNames":
+                    break
+                if not in_physical_names:
+                    continue
+                parts = line.split(maxsplit=2)
+                if len(parts) < 3:
+                    continue
+                try:
+                    dim = int(parts[0])
+                    tag = int(parts[1])
+                except ValueError:
+                    continue
+                name = parts[2].strip().strip('"')
+                if name:
+                    region_names.append((dim, tag, name))
+    except OSError:
+        pass
+    return region_names
+
+
+def _extract_netgen_region_names(ngmesh, mesh_dim):
+    """Extract region names from Netgen, converting codimension to dimension."""
+    region_names = []
+    for dim in range(mesh_dim + 1):
+        names = ngmesh.GetRegionNames(dim=dim)
+        for subdomain_id, name in enumerate(names, start=1):
+            if name != "":
+                region_names.append((dim, subdomain_id, name))
+    return region_names
+
+
+def _register_region_names(mesh, region_names):
+    all_region_names = mesh.comm.allgather(region_names)
+    seen = set()
+    for rank_region_names in all_region_names:
+        for dim, subdomain_id, name in rank_region_names:
+            key = (dim, subdomain_id, name)
+            if key in seen:
+                continue
+            seen.add(key)
+            mesh.rename_subdomain(dim, subdomain_id, name)
+
+
 class _Facets(object):
     """Wrapper class for facet interation information on a :func:`Mesh`
 
@@ -599,6 +667,8 @@ def __init__(self, topology_dm, name, reorder, sfXB, perm_is, distribution_name,
         self._shared_data_cache = defaultdict(dict)
         # Cell subsets for integration over subregions
         self._subsets = {}
+        # Human-readable region names, keyed by topological entity dimension.
+        self.region_names = defaultdict(dict)
         # A set of weakrefs to meshes that are explicitly labelled as being
         # parallel-compatible for interpolation/projection/supermeshing
         # To set, do e.g.
@@ -2876,6 +2946,75 @@ def __dir__(self):
         current = super(MeshGeometry, self).__dir__()
         return list(OrderedDict.fromkeys(dir(self.topology) + current))
 
+    _reserved_subdomain_names = frozenset({"on_boundary", "everywhere", "otherwise", "top", "bottom"})
+
+    def rename_subdomain(self, dim, subdomain_id, subdomain_name):
+        """Bind a string name to one or more integer subdomain ids."""
+        if not isinstance(subdomain_name, str):
+            raise TypeError("subdomain_name must be a string")
+        dim = int(dim)
+        mesh_dim = self.topological_dimension
+        if not 0 <= dim <= mesh_dim:
+            raise ValueError(f"Invalid entity dimension {dim} for mesh of dimension {mesh_dim}")
+
+        dim_region_names = self.topology.region_names[dim]
+        dim_region_names.setdefault(subdomain_name, []).extend(_as_subdomain_id_list(subdomain_id))
+
+    def parse_subdomain_id(self, dim, subdomain_id):
+        """Resolve string region names to integer subdomain ids for an entity dimension."""
+        dim = int(dim)
+        mesh_dim = self.topological_dimension
+        if not 0 <= dim <= mesh_dim:
+            raise ValueError(f"Invalid entity dimension {dim} for mesh of dimension {mesh_dim}")
+
+        dim_region_names = self.topology.region_names[dim]
+
+        def is_collection(value):
+            return not isinstance(value, str) and isinstance(value, Sequence)
+
+        def resolve_string(value):
+            if value in self._reserved_subdomain_names:
+                return (value,)
+            try:
+                return tuple(dim_region_names[value])
+            except KeyError:
+                available = sorted(dim_region_names)
+                raise ValueError(
+                    f"Subdomain region {value!r} not found in dimension {dim}. "
+                    f"Available names: {available}"
+                )
+
+        def resolve(value, flatten_named):
+            if isinstance(value, str):
+                resolved = resolve_string(value)
+                if flatten_named:
+                    return resolved
+                if len(resolved) == 1:
+                    return resolved[0]
+                return resolved
+            if is_collection(value):
+                resolved = []
+                for entry in value:
+                    if is_collection(entry):
+                        resolved.append(resolve(entry, flatten_named=False))
+                    else:
+                        entry_resolved = resolve(entry, flatten_named=flatten_named)
+                        if flatten_named and isinstance(entry_resolved, Sequence):
+                            resolved.extend(entry_resolved)
+                        else:
+                            resolved.append(entry_resolved)
+                return tuple(resolved)
+            return value
+
+        if isinstance(subdomain_id, str):
+            resolved = resolve_string(subdomain_id)
+            if resolved == (subdomain_id,):
+                return subdomain_id
+            return resolved
+        if is_collection(subdomain_id):
+            return resolve(subdomain_id, flatten_named=True)
+        return subdomain_id
+
     def mark_entities(self, f, label_value, label_name=None):
         """Mark selected entities.
 
@@ -3430,6 +3569,13 @@ def Mesh(meshfile, **kwargs):
             temp._did_reordering = mesh._did_reordering
             mesh = temp
 
+    if from_netgen:
+        region_names = _extract_netgen_region_names(meshfile, mesh.topological_dimension)
+        _register_region_names(mesh, region_names)
+    elif isinstance(meshfile, str) and meshfile.lower().endswith(".msh"):
+        region_names = _parse_gmsh_physical_names(meshfile)
+        _register_region_names(mesh, region_names)
+
     mesh.submesh_parent = submesh_parent
     mesh._tolerance = tolerance
     return mesh
@@ -4985,6 +5131,8 @@ def Submesh(mesh, subdim=None, subdomain_id=None, label_name=None, name=None, ig
     dim = plex.getDimension()
     if subdim not in {dim, dim - 1}:
         raise NotImplementedError(f"Found submesh dim ({subdim}) and parent dim ({dim})")
+    if subdomain_id is not None:
+        subdomain_id = mesh.parse_subdomain_id(subdim, subdomain_id)
     if subdomain_id is None:
         if label_name is not None:
             raise ValueError("subdomain_id=None requires label_name=None.")

firedrake/tsfc_interface.py

@@ -11,6 +11,7 @@
 import ufl
 import finat.ufl
 from ufl import conj, Form, ZeroBaseForm
+from ufl.measure import as_integral_type
 from .ufl_expr import TestFunction, extract_domains
 
 from tsfc import compile_form as original_tsfc_compile_form
@@ -150,7 +151,42 @@ def __init__(
 SplitKernel = collections.namedtuple("SplitKernel", ["indices", "kinfo"])
 
 
+def entity_dimension_from_integral_type(domain, integral_type):
+    integral_type = as_integral_type(integral_type)
+    topological_dimension = domain.topological_dimension
+    if integral_type == "cell":
+        return topological_dimension
+    if integral_type.startswith("exterior_facet") or integral_type.startswith("interior_facet"):
+        return topological_dimension - 1
+    if integral_type == "ridge":
+        return topological_dimension - 2
+    if integral_type == "vertex":
+        return 0
+    return topological_dimension
+
+
+def _resolve_region_names(form):
+    if not isinstance(form, Form):
+        return form
+
+    integrals = []
+    changed = False
+    for integral in form.integrals():
+        domain = integral.ufl_domain()
+        parse_subdomain_id = getattr(domain, "parse_subdomain_id", None)
+        if parse_subdomain_id is not None:
+            dim = entity_dimension_from_integral_type(domain, integral.integral_type())
+            subdomain_id = integral.subdomain_id()
+            new_subdomain_id = parse_subdomain_id(dim, subdomain_id)
+            if new_subdomain_id != subdomain_id:
+                integral = integral.reconstruct(subdomain_id=new_subdomain_id)
+                changed = True
+        integrals.append(integral)
+    return Form(integrals) if changed else form
+
+
 def _compile_form_hashkey(form, name, parameters=None, split=True, dont_split=(), diagonal=False):
+    form = _resolve_region_names(form)
     return (
         form.signature(),
         name,
@@ -210,6 +246,8 @@ def compile_form(form, name, parameters=None, split=True, dont_split=(), diagona
     if not isinstance(form, Form):
         raise RuntimeError("Unable to convert object to a UFL form: %s" % repr(form))
 
+    form = _resolve_region_names(form)
+
     if parameters is None:
         parameters = default_parameters["form_compiler"].copy()
     else: