Native Region Name Support in Firedrake

[pbrubeck]

Disclaimer: This architectural blueprint was co-authored using Gemini in collaboration with a core Firedrake developer.

---

Architecture Plan: Native Region Name Support in Firedrake

This development plan introduces native, human-readable region naming capabilities across Firedrake. By allowing users to seamlessly use string identifiers preserved from external mesh generators or explicitly assign custom semantic names to subdomains at runtime, it eliminates the cognitive overhead of tracking opaque numerical tags when defining boundary conditions, extracting submeshes, or specifying UFL variational forms.

---

Motivation: The Target User Experience

Firedrake users can use names natively within their mesh generators, which are automatically preserved and parsed when defining boundary conditions, submesh extractions, or integration domains in modern UFL measures.

from firedrake import *
import netgen.occ as ngocc

# 1. Generate geometry & name boundaries natively in Netgen
wp = ngocc.WorkPlane()
box = wp.Rectangle(1, 1).Face()
box.edges.Max(ngocc.Y).name = "inlet"
box.faces.name = "fluid_domain"
ngmesh = ngocc.OCCGeometry(box, dim=2).GenerateMesh(maxh=0.1)

# 2. Native Firedrake Mesh captures the names automatically
mesh = Mesh(ngmesh)

# 3. Explicitly map a single string name to multiple disjoint entity IDs

# 4. Clean, human-readable API usage
V = FunctionSpace(mesh, "CG", 1)

# DirichletBC automatically maps "inlet" to the correct face ID
bc = DirichletBC(V, 0.0, "inlet")  

# SubMesh automatically maps "fluid_domain" using the correct keyword argument
submesh = SubMesh(mesh, subdomain_id="fluid_domain")

# Form language integration: UFL Measures accept string region names directly
dx = Measure("dx", domain=mesh, subdomain_id="fluid_domain")
f = Function(V)
L = f * dx

---

Phase 1: Implementing the Naming & Registration Engine

Mesh.rename_subdomain serves as the sole gateway for mutating the region_names store. It handles both scalar IDs and iterables, treating the underlying values as flat list structures.

Expected Storage Structure

# Each integer dimension key maps to a dictionary of string names pointing to integer lists
{
    3: {"fluid_domain": [1, 2]}, 
    2: {"inlet": [10], "internal_obstacles": [12, 14, 15, 42]},
    1: {},
    0: {}
}

1. Naming and Registration API

# Inside firedrake/mesh.py -> Class MeshGeometry / Mesh

def rename_subdomain(self, dim, subdomain_id, subdomain_name):
    """
    Bind a custom string name to an integer subdomain ID (or multiple IDs) for a given dimension.
    """
    if not hasattr(self.topology, "region_names"):
        self.topology.region_names = {}
        
    mesh_dim = self.topological_dimension

    if not (0 <= dim <= mesh_dim):
        raise ValueError(f"Invalid entity dimension {dim} for mesh of dimension {mesh_dim}")

    if dim not in self.topology.region_names:
        self.topology.region_names[dim] = {}

    if subdomain_name not in self.topology.region_names[dim]:
        self.topology.region_names[dim][subdomain_name] = []

    # Dynamically extend or append depending on whether an iterable container is provided
    if isinstance(subdomain_id, (list, tuple, set)):
        self.topology.region_names[dim][subdomain_name].extend(subdomain_id)
    else:
        self.topology.region_names[dim][subdomain_name].append(subdomain_id)

2. Streamlined File Parsers

# Inside firedrake/mesh.py

def extract_netgen_region_names(ngmesh, mesh_dim):
    """Query a live Netgen mesh object using its native region name API."""
    extracted = []
    get_names_fn = getattr(ngmesh, "GetRegionNames", getattr(ngmesh, "getRegionNames", None))
    if not get_names_fn:
        return extracted

    # Sweep sequentially through all available dimensions
    for target_dim in range(mesh_dim + 1):
        names = get_names_fn(dim=target_dim)
        for i, name in enumerate(names, 1):
            if name:
                extracted.append((target_dim, i, name))
    return extracted


def parse_gmsh_physical_names(filename):
    """Parses the $PhysicalNames block from a Gmsh .msh file."""
    extracted = []
    try:
        with open(filename, 'r') as f:
            in_physical_names = False
            for line in f:
                line = line.strip()
                if line == "$PhysicalNames":
                    in_physical_names = True
                    continue
                elif line == "$EndPhysicalNames":
                    break
                if in_physical_names:
                    parts = line.split()
                    if len(parts) >= 3: 
                        dim = int(parts[0])
                        tag = int(parts[1])
                        name = parts[2].strip('"') 
                        extracted.append((dim, tag, name))
    except IOError:
        pass
    return extracted

3. Mesh Constructor Integration

# Inside firedrake/mesh.py -> Mesh()

# ... [Standard DMPlex/Mesh topology construction happens here] ...

mesh.topology.region_names = {d: {} for d in range(dim + 1)}

extracted_names = []
if isinstance(filename_or_mesh, str) and filename_or_mesh.endswith('.msh'):
    extracted_names = parse_gmsh_physical_names(filename_or_mesh)
elif hasattr(filename_or_mesh, "GetRegionNames") or hasattr(filename_or_mesh, "getRegionNames"):
    extracted_names = extract_netgen_region_names(filename_or_mesh, mesh_dim=dim)

# Populate storage uniformly using integer dimensions; scalar IDs safely append
for dim_val, tag, name in extracted_names:
    mesh.rename_subdomain(dim_val, tag, name)

---

Phase 2: Implementing Mesh.parse_subdomain_id

Because names point directly to integer arrays, parse_subdomain_id automatically flattens lookups when evaluating collections containing mixed strings and integers. It includes an explicit pass-through guard for standard fallback strings like "everywhere".

# Inside firedrake/mesh.py -> Class MeshGeometry / Mesh

def parse_subdomain_id(self, dim, subdomain_id):
    """
    Translate a string region name (or collection of names) into the corresponding
    integer subdomain ID(s) for a given entity dimension. Always flattens named lists.
    """
    meta = getattr(self.topology, "region_names", {})
    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_dict = meta.get(dim, {})

    def resolve_single(sd):
        if isinstance(sd, str):
            if sd in ("on_boundary", "everywhere"):
                return [sd]
            if sd in dim_dict:
                return dim_dict[sd]
            else:
                available = list(dim_dict.keys())
                raise ValueError(
                    f"Subdomain region '{sd}' not found in dimension {dim}. "
                    f"Available names: {available}"
                )
        return [sd]

    # Process iterables and flatten underlying integer lists seamlessly
    if isinstance(subdomain_id, (list, tuple, set)):
        resolved = []
        for sd in subdomain_id:
            resolved.extend(resolve_single(sd))
        return type(subdomain_id)(resolved)
    else:
        if isinstance(subdomain_id, str):
            if subdomain_id in ("on_boundary", "everywhere"):
                return subdomain_id
            return resolve_single(subdomain_id)
        return subdomain_id

---

Phase 3: High-Level API Integration (DirichletBC & SubMesh)

Both high-level APIs route their internal checks completely through the new centralized parser method on the mesh object using integer definitions.

# Inside firedrake/bcs.py -> DirichletBC.__init__

mesh = V.mesh()
facet_dim = mesh.topological_dimension - 1
sub_domain = mesh.parse_subdomain_id(facet_dim, sub_domain)

# Inside firedrake/submesh.py -> SubMesh.__init__

subdomain_id = mesh.parse_subdomain_id(subdim, subdomain_id)

---

Phase 4: Form Compiler Integration (TSFC Driver)

To support string region names passed into ufl.Measure, we hook into the form compilation pipeline right as an individual integral block enters TSFC. We determine the entity dimension from the UFL integral type and swap out string subdomains for their actual integer IDs.

# Inside firedrake/tsfc/driver.py -> compile_integral

def compile_integral(integral, coordinate_mapping, ...):
    """Compiles a single UFL integral using TSFC."""
    
    # 1. Extract the underlying Firedrake mesh object and current metadata
    mesh = integral.ufl_domain()
    subdomain_id = integral.subdomain_id()
    integral_type = integral.integral_type()
    
    # 2. Derive the correct topological dimension based on the integral type
    if integral_type == "cell":
        dim = mesh.topological_dimension
    elif integral_type in ("exterior_facet", "interior_facet"):
        dim = mesh.topological_dimension - 1
    elif integral_type == "vertex":
        dim = 0
    else:
        dim = mesh.topological_dimension

    # 3. Resolve potential string names into their corresponding numerical IDs
    new_subdomain_id = mesh.parse_subdomain_id(dim, subdomain_id)
    
    # 4. Reconstruct the immutable UFL integral if structural updates occurred
    if new_subdomain_id != subdomain_id:
        integral = integral.reconstruct(subdomain_id=new_subdomain_id)
        
    # ... [Proceed with standard TSFC translation and code generation] ...

---

Phase 5: Verification and Edge Cases

• Sequential Appends: Ensure that executing rename_subdomain multiple times for the same name and dimension correctly appends to the internal value list instead of replacing it.
• String Pass-Through Integrity: Verify that hardcoded strings like "on_boundary" and "everywhere" evaluate properly both as standalone arguments and when combined within collections.
• Compiler Isolation: Ensure that updating the subdomain_id on the UFL expression level leaves the remainder of the active variational form definition mathematically unchanged.