BiV Mesh from UK Biobank Atlas

This demo illustrates how to set up a simulation using a patient-specific Bi-Ventricular (BiV) geometry from the UK Biobank Atlas.

Key features of this example:

1. Geometry: Loading a BiV mesh from the UK Biobank statistical atlas.

2. Microstructure: Generating fiber and sheet orientations using the Laplace-Dirichlet Rule-Based (LDRB) algorithm via the external library fenicsx-ldrb. This provides more flexibility than standard built-in methods, allowing for distinct fiber angles in the LV and RV.

3. Boundary Conditions: Fixing the outflow tracts (Aortic, Pulmonary, Mitral, and Tricuspid valves) to prevent rigid body motion, instead of fixing the entire base.

---

Imports

from pathlib import Path
from mpi4py import MPI
import dolfinx
import ldrb
import numpy as np
import cardiac_geometries
import cardiac_geometries.geometry
import pulse

1. Geometry Generation

We generate the mesh using cardiac_geometries. We set create_fibers=False because we will generate a more advanced fiber architecture using ldrb in the next step.

outdir = Path("ukb_atlas")
outdir.mkdir(parents=True, exist_ok=True)
geodir = outdir / "geometry"
comm = MPI.COMM_WORLD
mode = -1
std = 0
char_length = 10.0

Generate base UKB mesh

geo = cardiac_geometries.mesh.ukb(
    outdir=geodir,
    comm=comm,
    mode=mode,
    std=std,
    case="ED",
    create_fibers=False,
    char_length_max=char_length,
    char_length_min=char_length,
    clipped=False,
)

INFO:ukb.atlas:Generating points from /github/home/.ukb/UKBRVLV.h5 using mode -1 and std 0.0

INFO:ukb.surface:Saved ukb_atlas/geometry/EPI_ED.stl

INFO:ukb.surface:Saved ukb_atlas/geometry/MV_ED.stl

INFO:ukb.surface:Saved ukb_atlas/geometry/AV_ED.stl

INFO:ukb.surface:Saved ukb_atlas/geometry/TV_ED.stl

INFO:ukb.surface:Saved ukb_atlas/geometry/PV_ED.stl

INFO:ukb.surface:Saved ukb_atlas/geometry/LV_ED.stl

INFO:ukb.surface:Saved ukb_atlas/geometry/RV_ED.stl

INFO:ukb.surface:Saved ukb_atlas/geometry/RVFW_ED.stl

INFO:ukb.mesh:Creating mesh for ED with char_length_max=10.0, char_length_min=10.0

0

INFO:ukb.mesh:Created mesh ukb_atlas/geometry/ED.msh

2025-12-31 10:16:42 [debug    ] Convert file ukb_atlas/geometry/ED.msh to dolfin

INFO:cardiac_geometries.geometry:Reading geometry from ukb_atlas/geometry

Info    : Reading 'ukb_atlas/geometry/ED.msh'...
Info    : 27 entities
Info    : 2562 nodes
Info    : 12699 elements
Info    : 8 parametrizations
Info    : [  0%] Processing parametrizations                                                                                
Info    : [ 10%] Processing parametrizations                                                                                
Info    : [ 20%] Processing parametrizations                                                                                
Info    : [ 30%] Processing parametrizations                                                                                
Info    : [ 40%] Processing parametrizations                                                                                
Info    : [ 60%] Processing parametrizations                                                                                
Info    : [ 70%] Processing parametrizations                                                                                
Info    : [ 80%] Processing parametrizations                                                                                
Info    : Done reading 'ukb_atlas/geometry/ED.msh'

2. Fiber Generation with LDRB

We use fenicsx-ldrb to generate the fiber field. This algorithm solves Laplace problems to define transmural and apicobasal directions.

We specify distinct angles for the Left Ventricle (LV) and Right Ventricle (RV), which allows for a more physiologically accurate representation (e.g., \(\alpha_{endo}^{RV} = 90^\circ\)).

system = ldrb.dolfinx_ldrb(
    mesh=geo.mesh,
    ffun=geo.ffun,
    markers=cardiac_geometries.mesh.transform_markers(geo.markers, clipped=True),
    alpha_endo_lv=60,
    alpha_epi_lv=-60,
    alpha_endo_rv=90,
    alpha_epi_rv=-25,
    beta_endo_lv=-20,
    beta_epi_lv=20,
    beta_endo_rv=0,
    beta_epi_rv=20,
    fiber_space="Quadrature_6",
)

INFO:ldrb.ldrb:Calculating scalar fields

INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers: 
lv: [1]
rv: [2]
epi: [7]
base: [5, 6, 4, 3]

INFO:ldrb.ldrb:  Num vertices: 2562

INFO:ldrb.ldrb:  Num cells: 7787

INFO:ldrb.ldrb:  Apex coord: (46.59, -17.77, -30.95)

INFO:ldrb.ldrb:
Calculating gradients

INFO:ldrb.ldrb:Compute fiber-sheet system

INFO:ldrb.ldrb:Angles: 

INFO:ldrb.ldrb:alpha: 
 endo_lv: 60
 epi_lv: -60
 endo_septum: 60
 epi_septum: -60
 endo_rv: 90
 epi_rv: -25

INFO:ldrb.ldrb:beta: 
 endo_lv: -20
 epi_lv: 20
 endo_septum: -20
 epi_septum: 20
 endo_rv: -20
 epi_rv: 20

We save the generated microstructure to the geometry folder so it can be automatically loaded by the cardiac_geometries.geometry.Geometry class later.

cardiac_geometries.fibers.utils.save_microstructure(
    geo.mesh,
    [system.f0, system.s0, system.n0],
    path=geodir / "geometry.bp",
)

3. Load Geometry and Visualization

We reload the geometry (now including the fibers we just saved).

geo = cardiac_geometries.geometry.Geometry.from_folder(
    comm=MPI.COMM_WORLD,
    folder=geodir,
)

INFO:cardiac_geometries.geometry:Reading geometry from ukb_atlas/geometry

Visualizing the Mesh

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    p = pyvista.Plotter()
    topology, cell_types, geometry_pts = dolfinx.plot.vtk_mesh(geo.mesh)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_pts)
    p.add_mesh(grid, show_edges=True, color="lightblue")
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot("ukb_mesh.png")

2025-12-31 10:16:46.772 (   1.008s) [    7F7CD646D140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=:99.0
INFO:trame_server.utils.namespace:Translator(prefix=None)

INFO:wslink.backends.aiohttp:awaiting runner setup

INFO:wslink.backends.aiohttp:awaiting site startup

INFO:wslink.backends.aiohttp:Print WSLINK_READY_MSG

INFO:wslink.backends.aiohttp:Schedule auto shutdown with timout 0

INFO:wslink.backends.aiohttp:awaiting running future

Visualizing Fibers

We create a glyph plot to visualize the generated fiber vectors \(\mathbf{f}_0\).

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    plotter = pyvista.Plotter()
    N = 10 # Only show every 10th point to avoid clutter
    points = geo.f0.function_space.tabulate_dof_coordinates()
    point_cloud = pyvista.PolyData(points[::N, :])
    f0_arr = geo.f0.x.array.reshape((points.shape[0], 3))
    point_cloud["fibers"] = f0_arr[::N, :]
    fibers = point_cloud.glyph(
        orient="fibers",
        scale=False,
        factor=5.0,
    )
    plotter.add_mesh(fibers, color="red")

    if not pyvista.OFF_SCREEN:
        plotter.show()
    else:
        fig_as_array = plotter.screenshot("fiber_ukb.png")

Visualizing Boundary Markers

The UK Biobank mesh comes with detailed surface markers, including specific tags for the outflow tracts (valves). We can visualize these to verify their location.

print("Available markers:", geo.markers)

Available markers: {'LV': [1, 2], 'RV': [2, 2], 'MV': [3, 2], 'AV': [4, 2], 'PV': [5, 2], 'TV': [6, 2], 'EPI': [7, 2], 'Wall': [9, 3]}

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    p = pyvista.Plotter()
    vtk_bmesh = dolfinx.plot.vtk_mesh(
        geo.mesh, geo.ffun.dim, geo.ffun.indices,
    )
    bgrid = pyvista.UnstructuredGrid(*vtk_bmesh)
    bgrid.cell_data["Facet tags"] = geo.ffun.values
    bgrid.set_active_scalars("Facet tags")
    p = pyvista.Plotter(window_size=[800, 800])
    p.add_mesh(bgrid, show_edges=True, cmap="tab20", show_scalar_bar=False)
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        figure = p.screenshot("facet_tags_ukb.png")

Convert to pulse.Geometry for simulation.

geometry = pulse.Geometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})

4. Constitutive Model

We use the Holzapfel-Ogden passive law and an Active Stress model.

material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
material = pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)

Ta = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
active_model = pulse.ActiveStress(geo.f0, activation=Ta, eta=0.3)
comp_model = pulse.Incompressible()

model = pulse.CardiacModel(
    material=material,
    active=active_model,
    compressibility=comp_model,
)

5. Boundary Conditions

Neumann BCs (Cavity Pressures)

We apply pressure loads to the LV and RV endocardial surfaces.

traction_lv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
neumann_lv = pulse.NeumannBC(traction=traction_lv, marker=geometry.markers["LV"][0])

traction_rv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
neumann_rv = pulse.NeumannBC(traction=traction_rv, marker=geometry.markers["RV"][0])

neumann = [neumann_lv, neumann_rv]

Robin BCs (Spring Support)

Elastic springs on the Epicardium to mimic the pericardial constraint.

robin_epi = pulse.RobinBC(
    value=pulse.Variable(
        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
        "Pa / m",
    ),
    marker=geometry.markers["EPI"][0],
)

Dirichlet BC (Fixed Outflow Tracts)

Instead of fixing the entire base plane (which might be unphysical if the base moves), we fix the displacement at the four valves:

• PV: Pulmonary Valve

• TV: Tricuspid Valve

• AV: Aortic Valve

• MV: Mitral Valve

This mimics the anatomical anchoring of the heart to the great vessels.

def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
    # Find facets for the outflow tracts
    bcs = []
    u_bc = dolfinx.fem.Function(V)
    u_bc.x.array[:] = 0
    # Iterate over the valve markers provided by the UKB mesh
    for marker in ["PV", "TV", "AV", "MV"]:
        facets = geometry.facet_tags.find(geometry.markers[marker][0])
        dofs = dolfinx.fem.locate_dofs_topological(V, 2, facets)
        bcs.append(dolfinx.fem.dirichletbc(u_bc, dofs))
    return bcs

bcs = pulse.BoundaryConditions(
    neumann=neumann,
    dirichlet=(dirichlet_bc,),
    robin=(robin_epi,),
)

6. Solving the Problem

We set up the StaticProblem and run the simulation phases.

problem = pulse.StaticProblem(
    model=model,
    geometry=geometry,
    bcs=bcs,
)

Phase 1: Passive Inflation

We ramp up the pressure in both ventricles.

vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, outdir / "displacement.bp", [problem.u], engine="BP4")
vtx.write(0.0)
pressures = [0.5, 1.0, 1.5]  #  kPa
for i, plv in enumerate(pressures, start=1):
    print(f"Solving for pressure: {plv} kPa")
    traction_lv.assign(plv)
    traction_rv.assign(plv * 0.5)  # Assume RV pressure is half of LV
    problem.solve()
    vtx.write(float(i))

Solving for pressure: 0.5 kPa

INFO:scifem.solvers:Newton iteration 1: r (abs) = 131706.29379412797 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 22890.33237766097 (tol=1e-06), r (rel) = 0.17379831835099074 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 12844.138552130093 (tol=1e-06), r (rel) = 0.09752106890356321 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 3713.256511185282 (tol=1e-06), r (rel) = 0.02819346292584565 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 6418.689434539904 (tol=1e-06), r (rel) = 0.04873487249267716 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 6060.284305099953 (tol=1e-06), r (rel) = 0.046013627219462054 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 7: r (abs) = 3531.130659936827 (tol=1e-06), r (rel) = 0.026810644793151563 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 8: r (abs) = 1341.1705604197723 (tol=1e-06), r (rel) = 0.01018304077796142 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 9: r (abs) = 486.97753107488876 (tol=1e-06), r (rel) = 0.0036974507219532764 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 10: r (abs) = 636.708858264831 (tol=1e-06), r (rel) = 0.004834308520290457 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 11: r (abs) = 225.51996223022877 (tol=1e-06), r (rel) = 0.0017122944981106394 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 12: r (abs) = 11.574098148950824 (tol=1e-06), r (rel) = 8.787809462653672e-05 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 13: r (abs) = 0.018108278733581108 (tol=1e-06), r (rel) = 1.374898511827113e-07 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 14: r (abs) = 3.9997058398017985e-08 (tol=1e-06), r (rel) = 3.036837287406931e-13 (tol=1e-10)

Solving for pressure: 1.0 kPa

INFO:scifem.solvers:Newton iteration 1: r (abs) = 34504.94276080669 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 4605.563541382535 (tol=1e-06), r (rel) = 0.13347547258110162 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 1597.0259713501912 (tol=1e-06), r (rel) = 0.04628397683256595 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 89.73981945112685 (tol=1e-06), r (rel) = 0.0026007815770979943 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.2990810598656954 (tol=1e-06), r (rel) = 3.764912954271842e-05 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 0.00031223869397991715 (tol=1e-06), r (rel) = 9.04910047654336e-09 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 7: r (abs) = 2.2108438212418544e-10 (tol=1e-06), r (rel) = 6.4073249927343655e-15 (tol=1e-10)

Solving for pressure: 1.5 kPa

INFO:scifem.solvers:Newton iteration 1: r (abs) = 38833.37853524775 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 3904.782245448592 (tol=1e-06), r (rel) = 0.10055221545826493 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 423.9886940349196 (tol=1e-06), r (rel) = 0.010918151086186831 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 9.952663452068716 (tol=1e-06), r (rel) = 0.00025629146439151614 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.008251729527635632 (tol=1e-06), r (rel) = 2.1249064178502562e-07 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 9.851196345993435e-09 (tol=1e-06), r (rel) = 2.536785805811832e-13 (tol=1e-10)

Visualize Passive Inflation

try:
    import pyvista
except ImportError:
    pass
else:
    V = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 1, (geometry.mesh.geometry.dim,)))
    uh = dolfinx.fem.Function(V)
    uh.interpolate(problem.u)

    p = pyvista.Plotter()
    topology, cell_types, geometry_data = dolfinx.plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_data)
    grid["u"] = uh.x.array.reshape((geometry_data.shape[0], 3))

    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")
    warped = grid.warp_by_vector("u", factor=1.0)
    p.add_mesh(warped, show_edges=False, opacity=0.5, color="firebrick", label="Inflated")
    p.add_legend()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot(outdir / "pressure.png")

Phase 2: Active Contraction

We hold the pressure constant and ramp up the active tension \(T_a\).

active_tensions = [3.0, 5.0, 10.0]  # kPa
for i, ta in enumerate(active_tensions, start=len(pressures) + 1):
    print(f"Solving for active tension: {ta} kPa")
    Ta.assign(ta)
    problem.solve()
    vtx.write(float(i))

Solving for active tension: 3.0 kPa

INFO:scifem.solvers:Newton iteration 1: r (abs) = 45714.40016017963 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 13045.148172979134 (tol=1e-06), r (rel) = 0.2853619018792759 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 1292.8517702603283 (tol=1e-06), r (rel) = 0.028281061672695652 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 58.5636129963535 (tol=1e-06), r (rel) = 0.0012810758271168657 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.20987852234164225 (tol=1e-06), r (rel) = 4.591081182433643e-06 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 1.0516554262196463e-05 (tol=1e-06), r (rel) = 2.3004904855685059e-10 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 7: r (abs) = 2.2697237115200444e-10 (tol=1e-06), r (rel) = 4.965008189032586e-15 (tol=1e-10)

Solving for active tension: 5.0 kPa

INFO:scifem.solvers:Newton iteration 1: r (abs) = 35476.677161589854 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 3899.0798069886455 (tol=1e-06), r (rel) = 0.10990543982541097 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 192.33672268797912 (tol=1e-06), r (rel) = 0.0054214976733001265 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 1.284932692520565 (tol=1e-06), r (rel) = 3.621908237538507e-05 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 8.876500286723684e-05 (tol=1e-06), r (rel) = 2.5020664269916908e-09 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 2.2221870395309057e-10 (tol=1e-06), r (rel) = 6.263796999389896e-15 (tol=1e-10)

Solving for active tension: 10.0 kPa

INFO:scifem.solvers:Newton iteration 1: r (abs) = 91025.218713517 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 14517.385222732728 (tol=1e-06), r (rel) = 0.1594875071756014 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 1217.6769347795 (tol=1e-06), r (rel) = 0.013377357967267134 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 24.5127518895099 (tol=1e-06), r (rel) = 0.0002692962701540845 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.023371833177595468 (tol=1e-06), r (rel) = 2.5676217544891013e-07 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 3.6644060633248455e-08 (tol=1e-06), r (rel) = 4.0257042115524095e-13 (tol=1e-10)

vtx.close()

# Visualize Active Contraction
try:
    import pyvista
except ImportError:
    pass
else:
    uh.interpolate(problem.u)
    grid["u"] = uh.x.array.reshape((geometry_data.shape[0], 3))

    p = pyvista.Plotter()
    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")
    warped = grid.warp_by_vector("u", factor=1.0)
    p.add_mesh(warped, show_edges=False, opacity=0.5, color="red", label="Contracted")
    p.add_legend()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot(outdir / "active.png")