Interface markers and mesh refinement

We start by reading in the mesh generated in Generating brain meshes.

from mpi4py import MPI
import dolfinx
from pathlib import Path
import os

folder = Path(os.environ["WILDFENICS_DATA_PATH"])
assert folder.exists(), "Could not find surface files"

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, folder / "brain.xdmf", "r") as xdmf:
    brain_with_csf = xdmf.read_mesh(dolfinx.cpp.mesh.GhostMode.none)
    region_markers = xdmf.read_meshtags(brain_with_csf, name="mesh_tags")

The meshes generated with Wildmeshing has a more or less uniform cell size. In some cases this is ok, but for fluid flow simulations, we need a finer mesh in the very thin regions, such as the cerebral aquaduct and in the subarachnoid spaces.

The following maps the various regions to an integer marker

subdomain_map = {
    "PAR": (2,),
    "SAS": (1,),
    "LV": (3,),
    "V34": (4,),
}

import pyvista

pyvista.global_theme.allow_empty_mesh = True
pv_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(brain_with_csf))
pv_grid.cell_data["marker"] = region_markers.values
clipped = pv_grid.clip("y")
pyvista.start_xvfb(1.0)
plotter = pyvista.Plotter()
plotter.add_mesh(clipped, categories=True)
plotter.show()

The fluid regions can be gathered in a sequence of integers as

fluid_domains = subdomain_map["LV"] + subdomain_map["SAS"] + subdomain_map["V34"]

Marking interfaces

With the refined mesh, we are ready to create markers between the various interfaces in the brain. As in [CKMR25], we will divide the brain parenchyma into the following interfaces:

• Lateral ventricle - brain parenchyma ("LV_PAR"): The interface between the CSF filled regions called the lateral ventricles and the brain tissue.

• 3rd/4th ventricle - brain parenchyma ("V34_PAR"): The interface between the CSF filled regions called the 3rd and 4th ventricles and the brain tissue

• Brain parenchyma - Subarachnoid space ("PAR_SAS"): The interface between the remaining CSF spaces in the brain and the brain tissue.

• Subarachnoid space - Upper dura matter ("AM_U"): The interface between the SAS and the upper part of the dura matter, called the arachnoid matter.

• Subarachnoid space - Lower dura matter ("AM_L"): The interface between the SAS and the lower part of the dura matter.

• Spinal subarachnoid space ("SSAS"): The interface from the CSF towards the spinal subarachnoid space. We define the following map that assigns an integer tag to each of these regions

interface_map = {
    "LV_PAR": 1,
    "V34_PAR": 2,
    "PAR_SAS": 5,
    "AM_U": 3,
    "AM_L": 4,
    "SSAS": 6,
}

We define the divide between the upper and lower arachnoid matter with the following (patient specific) geometrical function

def AM_U_function(x, upper_skull_z=0.027):
    """Indicator function that returns True whenever the point x is in the AM_U region"""
    return x[2] - 0.8 * x[1] > upper_skull_z

To obtain the SSAS, which is not marked by wildmeshing, we use that we want this surface to satisfy the following constraints:

\[\begin{split} \begin{align*} \mathbf{x}_M&:=(x_m, y_m, z_m)\\ \mathbf{n}(\mathbf{x}_M)\cdot \mathbf{v} &> \epsilon\\ z_m + \gamma &> H \end{align*} \end{split}\]

where \(H\) is a prescribed scalar value and \(\mathbf{v}\) a prescribed vector. In our example above, we choose \(\epsilon=10^{-16}\), \(\gamma=10^{-3}\), \(\mathbf{v}=(0,0-1)\) and \(H=\min_{(x,y,z)\in\Omega}z\).

import ufl
import basix
import numpy as np
import numpy.typing as npt
import scifem
import typing


def locate_ssas(
    mesh: dolfinx.mesh.Mesh,
    cell_marker: dolfinx.mesh.MeshTags,
    tags: tuple[int, ...],
    condition: typing.Callable[[npt.NDArray[np.floating]], npt.NDArray[np.bool_]],
    v: tuple[float, ...],
    epsilon: float = 1.0e-16,
) -> npt.NDArray[np.floating]:
    """Given a `mesh`, and a normal vector `v`, find all exterior facets whose vertices
    all satisfies `condition([x,y,z])` and whose midpoint $\mathbf{x}_m=(x_m, y_m, z_m)$ satisfies
    $v \cdot n(x_m, y_m, z_m) > \epsilon$.

    Args:
        mesh: The mesh
        cell_marker: A set of tagged cells
        tags: The tags to consider as a single domain
        v: The vector to compare the normal to.
        epsilon: Threshold for the dot product.

    Returns:
        List of facets (index local to process) that satisfy the condition.
    """
    tdim = mesh.topology.dim
    fdim = tdim - 1
    if len(v) != (gdim := mesh.geometry.dim):
        raise ValueError(f"n_ssas must be of length {gdim}, but got {len(v)}")
    mesh.topology.create_connectivity(tdim - 1, tdim)

    subdomain_exterior_facets = scifem.mesh.compute_subdomain_exterior_facets(
        mesh, cell_marker, tags
    )

    # Filter with respect to midpoint
    facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, condition)

    # Get intersection of subdomain exterior facets and those satisfying the midpoint condition
    conditional_facets = np.intersect1d(subdomain_exterior_facets, facets)
    try:
        integration_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet,
            mesh.topology,
            conditional_facets,
            fdim,
        )
    except TypeError:
        integration_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet, mesh.topology, facets
        )
    facet_types = set(basix.cell.subentity_types(mesh.basix_cell())[fdim])
    assert len(facet_types) == 1, "All facets must have the same topology"

    reference_facet_point, _ = basix.make_quadrature(facet_types.pop(), 0)
    # Create expression for normal and evaluate
    nh = ufl.FacetNormal(mesh)
    normal_expr = dolfinx.fem.Expression(nh, reference_facet_point)
    n_evaluated = normal_expr.eval(mesh, integration_entities)

    # Compute dot product and filter
    n_dot = np.dot(n_evaluated, v)
    return conditional_facets[np.flatnonzero(n_dot > epsilon)]

With this definition at hand, we can define the various interface markers. We start by finding all facets in the mesh that is exterior, i.e. is only connected to a single cell. These facets will be the union of "AM_U" and "AM_L".

tdim = brain_with_csf.topology.dim
brain_with_csf.topology.create_connectivity(tdim - 1, tdim)
AM_UL = dolfinx.mesh.exterior_facet_indices(brain_with_csf.topology)

Next, we can find the different interfaces between the various regions in the brain.

import numpy as np

LV_PAR = scifem.mesh.find_interface(
    region_markers, subdomain_map["PAR"], subdomain_map["LV"]
)
V34_PAR = scifem.mesh.find_interface(
    region_markers, subdomain_map["PAR"], subdomain_map["V34"]
)
PAR_SAS = scifem.mesh.find_interface(
    region_markers, subdomain_map["PAR"], subdomain_map["SAS"]
)
AM_U = dolfinx.mesh.locate_entities_boundary(brain_with_csf, tdim - 1, AM_U_function)


# Find smallest z coordinate of the mesh across all processes
def midpoint_condition(z_max, gamma):
    return lambda x: x[2] < z_max + gamma


Z_min = brain_with_csf.comm.allreduce(brain_with_csf.geometry.x[:, 2].min(), op=MPI.MIN)
SSAS = locate_ssas(
    brain_with_csf,
    cell_marker=region_markers,
    condition=midpoint_condition(Z_min, 1e-3),
    tags=subdomain_map["LV"] + subdomain_map["SAS"] + subdomain_map["V34"],
    v=(0, 0, -1),
)

With all the various facets of interest located, we can store them in an distributed array.

facet_map = brain_with_csf.topology.index_map(tdim - 1)
facet_vec = dolfinx.la.vector(facet_map)
facet_marker = facet_vec.array
facet_marker[:] = -1
facet_marker[AM_UL] = interface_map["AM_L"]
# NOTE: Upper skull should always happen after lower skull
facet_marker[AM_U] = interface_map["AM_U"]
facet_marker[PAR_SAS] = interface_map["PAR_SAS"]
facet_marker[LV_PAR] = interface_map["LV_PAR"]
facet_marker[V34_PAR] = interface_map["V34_PAR"]
facet_marker[SSAS] = interface_map["SSAS"]
facet_vec.scatter_forward()

We filter out unmarked facets, which are those interior to each subdomain.

interface_marker = facet_marker.astype(np.int32)
facet_pos = interface_marker != -1
interface_tag = dolfinx.mesh.meshtags(
    brain_with_csf,
    tdim - 1,
    np.flatnonzero(facet_pos),
    interface_marker[facet_pos],
)
interface_tag.name = "interfaces_and_boundaries"

filename = folder / "marked_brain.xdmf"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "w") as xdmf:
    xdmf.write_mesh(brain_with_csf)
    xdmf.write_meshtags(region_markers, brain_with_csf.geometry)
    xdmf.write_meshtags(interface_tag, brain_with_csf.geometry)

We store the interface markers along with the mesh

Mesh refinement

For realistic fluid flow modelling, we need to refine the mesh, especially in the narrow aquaduct connecting the venctriles with the subarachnoid space. We will refine all cells (tetrahedra) whose facet is considered external, i.e. it does not connect to another fluid cell. We create the following convenience function to locate these cells. We will transfer the interface markers from the coarse grid to the finer grid.

def compute_subdomain_exterior_cells(
    mesh: dolfinx.mesh.Mesh, ct: dolfinx.mesh.MeshTags, markers: tuple[int, ...]
):
    """Compute the exterior boundary of a set of subdomains.

    Args:
        mesh: Mesh to extract subdomains from
        ct: MeshTags object marking subdomains
        markers: The tags making up the "new" mesh
    Returns:
        Cells which has a facet on the exterior boundary of the subdomains.
    """
    # Find facets that are considered exterior
    subdomain_exterior_facets = scifem.mesh.compute_subdomain_exterior_facets(
        mesh, ct, markers
    )
    tdim = mesh.topology.dim
    assert ct.dim == tdim
    sub_cells = dolfinx.mesh.compute_incident_entities(
        mesh.topology,
        subdomain_exterior_facets,
        tdim - 1,
        tdim,
    )
    full_subdomain = ct.indices[
        np.isin(ct.values, np.asarray(markers, dtype=ct.values.dtype))
    ]
    cell_map = mesh.topology.index_map(tdim)
    return scifem.mesh.reverse_mark_entities(
        cell_map, np.intersect1d(full_subdomain, sub_cells)
    )

With this function at hand we can locally refine the mesh

# Refine parent mesh within ventricles
num_refinements = 2
if num_refinements == 0:
    refined_mesh = brain_with_csf
    refined_regions = region_markers
    num_cells = brain_with_csf.topology.index_map(3).size_global
    print(f"Mesh, #Cells={num_cells}")
for i in range(num_refinements):
    # Refine parent mesh within ventricles
    refine_cells = region_markers.indices[
        np.isin(
            region_markers.values,
            np.asarray(subdomain_map["V34"] + subdomain_map["LV"]),
        )
    ]

    # Find all cells associated with outer boundary (dura) and refine the cells they correspond to
    brain_with_csf.topology.create_connectivity(
        brain_with_csf.topology.dim - 1, brain_with_csf.topology.dim
    )
    fmap = brain_with_csf.topology.index_map(brain_with_csf.topology.dim - 1)
    exterior_facet_indices = scifem.mesh.reverse_mark_entities(
        fmap, dolfinx.mesh.exterior_facet_indices(brain_with_csf.topology)
    )
    boundary_cells = dolfinx.mesh.compute_incident_entities(
        brain_with_csf.topology,
        exterior_facet_indices,
        brain_with_csf.topology.dim - 1,
        brain_with_csf.topology.dim,
    )

    fluid_boundary_cells = compute_subdomain_exterior_cells(
        brain_with_csf, region_markers, fluid_domains
    )

    # For any further refinement, only refine the boundary of the fluid domains, not the interior
    if i < 1:
        cells_to_refine = np.unique(
            np.hstack([boundary_cells, fluid_boundary_cells, refine_cells])
        ).astype(np.int32)

    else:
        cells_to_refine = refine_cells

    edges_to_refine = dolfinx.mesh.compute_incident_entities(
        brain_with_csf.topology, cells_to_refine, brain_with_csf.topology.dim, 1
    )
    edge_map = brain_with_csf.topology.index_map(1)
    edges_to_refine = scifem.mesh.reverse_mark_entities(edge_map, edges_to_refine)
    refined_mesh, parent_cell, parent_facet = dolfinx.mesh.refine(
        brain_with_csf,
        edges_to_refine,
        partitioner=None,
        option=dolfinx.mesh.RefinementOption.parent_cell_and_facet,
    )
    refined_regions = dolfinx.mesh.transfer_meshtag(
        region_markers, refined_mesh, parent_cell
    )
    refined_mesh.topology.create_connectivity(
        refined_mesh.topology.dim - 1, refined_mesh.topology.dim
    )
    refined_interfaces = dolfinx.mesh.transfer_meshtag(
        interface_tag, refined_mesh, parent_cell, parent_facet
    )
    refined_regions.name = region_markers.name
    refined_interfaces.name = interface_tag.name
    interface_tag = refined_interfaces
    brain_with_csf = refined_mesh
    region_markers = refined_regions
    num_cells = brain_with_csf.topology.index_map(3).size_global
    print(f"Mesh refinement {i + 1}, #Cells={num_cells}")

Mesh refinement 1, #Cells=1431351

Mesh refinement 2, #Cells=1745754

We store the refined mesh and the corresponding cell and facet markers.

filename = folder / "refined_marked_brain.xdmf"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "w") as xdmf:
    xdmf.write_mesh(brain_with_csf)
    xdmf.write_meshtags(region_markers, brain_with_csf.geometry)
    xdmf.write_meshtags(interface_tag, brain_with_csf.geometry)

We visualize the surface markers

import pyvista

pv_grid = pyvista.UnstructuredGrid(
    *dolfinx.plot.vtk_mesh(brain_with_csf, tdim - 1, interface_tag.indices)
)
pv_grid.cell_data["marker"] = interface_tag.values
clipped = pv_grid.clip("y")
pyvista.start_xvfb(1.0)
plotter = pyvista.Plotter()
plotter.add_mesh(clipped, categories=True)
plotter.show()

References

[CKMR25]

Marius Causemann, Miroslav Kuchta, Rami Masri, and Marie E. Rognes. In-silico molecular enrichment and clearance of the human intracranial space. bioRxiv, 2025. doi:10.1101/2025.01.30.635680.