Mesh refinement and interface markers

Mesh refinement and interface markers#

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"]

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.

import scifem
import numpy as np


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 our 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, _ = dolfinx.mesh.refine(
        brain_with_csf,
        edges_to_refine,
        partitioner=None,
        option=dolfinx.mesh.RefinementOption.parent_cell,
    )
    refined_regions = dolfinx.mesh.transfer_meshtag(
        region_markers, refined_mesh, parent_cell
    )
    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=1405474

Mesh refinement 2, #Cells=1724676

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.

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,
}

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

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.

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)

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_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"

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

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, folder / "refined_brain.xdmf", "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.

Mesh refinement and interface markers

Contents

Mesh refinement and interface markers#

References#