Rotated UK Biobank Atlas Geometry#
This demo shows how to load a patient-specific geometry from the UK Biobank Atlas, rotate it to align the basal plane with a coordinate axis, and perform a cardiac cycle simulation.
Motivation: Geometry Alignment#
Patient-specific meshes often come in arbitrary orientations based on the imaging coordinate system. However, for defining boundary conditions—especially at the base—it is often convenient to align the geometry such that the base normal points along a principal axis (e.g., the X-axis).
In this example:
We generate a Bi-Ventricular (BiV) mesh from the UK Biobank atlas using
cardiac-geometries.We rotate the mesh so the base normal aligns with \(\mathbf{e}_1 = (1, 0, 0)\).
We apply a Dirichlet boundary condition fixing \(u_x = 0\) on the base. Since the base is now perpendicular to X, this acts as a “sliding” condition, allowing the base to expand/contract in the Y-Z plane while preventing longitudinal rigid body motion.
Imports#
from pathlib import Path
from mpi4py import MPI
import dolfinx
from dolfinx import log
import cardiac_geometries
import cardiac_geometries.geometry
import pulse
1. Geometry Generation and Rotation#
We generate the mesh and then apply a rotation. The rotate method in cardiac_geometries
takes a target_normal vector. It computes the rotation matrix required to align the
average normal of the surface marked by base_marker with this target vector.
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=True,
char_length_max=char_length,
char_length_min=char_length,
clipped=True,
fiber_angle_endo=60,
fiber_angle_epi=-60,
fiber_space="Quadrature_6",
)
INFO:ukb.atlas:Downloading https://www.cardiacatlas.org/share/download.php?id=60&token=AR3JSoaxJ9Ev9n8QAkvV4BHJUniyttqm&download to /github/home/.ukb/UKBRVLV.zip. This may take a while.
INFO:ukb.atlas:Done downloading.
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.clip:Folder: ukb_atlas/geometry
INFO:ukb.clip:Case: ED
INFO:ukb.clip:Origin: [-13.612554383622273, 18.55767189380559, 15.135103714006394]
INFO:ukb.clip:Normal: [-0.7160843664428893, 0.544394641424108, 0.4368725838557541]
INFO:ukb.clip:Reading ukb_atlas/geometry/LV_ED.stl
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.clip:Saved ukb_atlas/geometry/lv_clipped.ply
INFO:ukb.clip:Reading ukb_atlas/geometry/RV_ED.stl
INFO:ukb.clip:Reading ukb_atlas/geometry/RVFW_ED.stl
INFO:ukb.clip:Merging RV and RVFW
INFO:ukb.clip:Smoothing RV
INFO:ukb.clip:Saving ukb_atlas/geometry/rv_clipped.ply
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.clip:Reading ukb_atlas/geometry/EPI_ED.stl
INFO:ukb.clip:Saving ukb_atlas/geometry/epi_clipped.ply
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
INFO:ukb.mesh:Creating clipped mesh for ED with char_length_max=10.0, char_length_min=10.0
0
INFO:ukb.mesh:Created mesh ukb_atlas/geometry/ED_clipped.msh
2025-12-31 10:14:33 [debug ] Convert file ukb_atlas/geometry/ED_clipped.msh to dolfin
Info : Reading 'ukb_atlas/geometry/ED_clipped.msh'...
Info : 11 entities
Info : 697 nodes
Info : 3440 elements
Info : 3 parametrizations
Info : [ 0%] Processing parametrizations
Info : [ 10%] Processing parametrizations
Info : [ 40%] Processing parametrizations
Info : Done reading 'ukb_atlas/geometry/ED_clipped.msh'
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
lv: [1]
rv: [2]
epi: [3]
base: [4]
INFO:ldrb.ldrb: Num vertices: 697
INFO:ldrb.ldrb: Num cells: 2156
INFO:ldrb.ldrb: Apex coord: (50.00, -13.53, -29.56)
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: 60
epi_rv: -60
INFO:ldrb.ldrb:beta:
endo_lv: 0
epi_lv: 0
endo_septum: 0
epi_septum: 0
endo_rv: 0
epi_rv: 0
INFO:cardiac_geometries.geometry:Reading geometry from ukb_atlas/geometry
# Rotate the geometry so the Base normal points in the X-direction (1, 0, 0)
rotated_geo = geo.rotate(target_normal=[1.0, 0.0, 0.0], base_marker="BASE")
rot_geodir = outdir / "geometry_rotated"
rotated_geo.save_folder(folder=rot_geodir)
2025-12-31 10:14:40 [info ] Rotated geometry. Base normal [ 0.71608436 -0.54439464 -0.43687259] aligned to [1.0, 0.0, 0.0]
2025-12-31 10:14:40 [debug ] Rotation matrix:
[[ 0.71608436 -0.54439464 -0.43687259]
[ 0.54439464 0.82730131 -0.1385894 ]
[ 0.43687259 -0.1385894 0.88878305]]
INFO:cardiac_geometries.geometry:Geometry saved to ukb_atlas/geometry_rotated
2. Load Geometry and Visualization#
We load the original and rotated geometries to compare them.
geo_orig = cardiac_geometries.geometry.Geometry.from_folder(
comm=MPI.COMM_WORLD,
folder=geodir,
)
geo_rot = cardiac_geometries.geometry.Geometry.from_folder(
comm=MPI.COMM_WORLD,
folder=rot_geodir,
)
INFO:cardiac_geometries.geometry:Reading geometry from ukb_atlas/geometry
INFO:cardiac_geometries.geometry:Reading geometry from ukb_atlas/geometry_rotated
# Visualize using PyVista (Optional)
try:
import pyvista
except ImportError:
print("Pyvista is not installed")
else:
p = pyvista.Plotter()
# Original Mesh (Red)
topology, cell_types, geometry_pts = dolfinx.plot.vtk_mesh(geo_orig.mesh)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_pts)
p.add_mesh(grid, show_edges=True, color="red", opacity=0.3, label="Original")
# Rotated Mesh (Yellow)
topology_rot, cell_types_rot, geometry_rot = dolfinx.plot.vtk_mesh(rotated_geo.mesh)
grid_rot = pyvista.UnstructuredGrid(topology_rot, cell_types_rot, geometry_rot)
p.add_mesh(grid_rot, show_edges=True, color="yellow", opacity=0.5, label="Rotated (Base -> X)")
p.add_legend()
axes = pyvista.CubeAxesActor(camera=p.camera)
axes.bounds = grid_rot.bounds
p.add_actor(axes)
p.view_zy() # View along the X-axis
if not pyvista.OFF_SCREEN:
p.show()
else:
p.screenshot("rotated_mesh.png")
2025-12-31 10:14:41.837 ( 9.155s) [ 7F6EED3F0140]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
# Convert the rotated geometry to a pulse.Geometry object for simulation
geometry = pulse.Geometry.from_cardiac_geometries(geo_rot, metadata={"quadrature_degree": 4})
3. Constitutive Model#
We use the Holzapfel-Ogden passive law and an Active Stress model.
Note that we use the fiber fields (f0, s0) from the rotated geometry geo_rot.
material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
material = pulse.HolzapfelOgden(f0=geo_rot.f0, s0=geo_rot.s0, **material_params)
Ta = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
# eta=0.3 adds 30% transverse active stress
active_model = pulse.ActiveStress(geo_rot.f0, activation=Ta, eta=0.3)
comp_model = pulse.Incompressible()
model = pulse.CardiacModel(
material=material,
active=active_model,
compressibility=comp_model,
)
4. 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 and Base to mimic the pericardium and surrounding tissue.
robin_epi = pulse.RobinBC(
value=pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
"Pa / m",
),
marker=geometry.markers["EPI"][0],
)
robin_base = pulse.RobinBC(
value=pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e6)),
"Pa / m",
),
marker=geometry.markers["BASE"][0],
)
Dirichlet BC (Sliding Base)#
Because we rotated the mesh, the basal plane is now perpendicular to the X-axis. We constrain the displacement in the x-direction (\(u_x = 0\)) on the Base.
Why? This prevents the heart from flying away in the X-direction (Rigid Body Motion).
Effect: Since the normal is X, fixing \(u_x\) allows the base to expand and contract freely in the Y and Z directions (in-plane), which is a common approximation for the AV-plane motion constraint.
def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
# Find facets for the BASE marker
facets = geometry.facet_tags.find(geometry.markers["BASE"][0])
# Locate degrees of freedom for the x-component (sub(0)) on these facets
dofs = dolfinx.fem.locate_dofs_topological(V.sub(0), 2, facets)
# Constrain u_x to 0.0
return [dolfinx.fem.dirichletbc(0.0, dofs, V.sub(0))]
bcs = pulse.BoundaryConditions(
neumann=neumann,
dirichlet=(dirichlet_bc,),
robin=(robin_epi, robin_base),
)
5. Solving the Problem#
We set up the StaticProblem. We do not use the automatic base_bc parameter because we
explicitly provided our own Dirichlet condition in bcs.
problem = pulse.StaticProblem(
model=model,
geometry=geometry,
bcs=bcs,
)
log.set_log_level(log.LogLevel.INFO)
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))
INFO:scifem.solvers:Newton iteration 1: r (abs) = 68278.83879351504 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
Solving for pressure: 0.5 kPa
INFO:scifem.solvers:Newton iteration 2: r (abs) = 4829.736586642113 (tol=1e-06), r (rel) = 0.07073548220771485 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 3354.2145836846 (tol=1e-06), r (rel) = 0.04912524353011075 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 526.6951589154152 (tol=1e-06), r (rel) = 0.007713885710742922 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 19.581536726375578 (tol=1e-06), r (rel) = 0.00028678778187181744 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 0.08349779315314725 (tol=1e-06), r (rel) = 1.2228941591355486e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 7: r (abs) = 4.627378945738831e-06 (tol=1e-06), r (rel) = 6.777178738690455e-11 (tol=1e-10)
Solving for pressure: 1.0 kPa
INFO:scifem.solvers:Newton iteration 1: r (abs) = 16432.724287441222 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 1819.2706642969529 (tol=1e-06), r (rel) = 0.11071022871644831 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 399.69950695632167 (tol=1e-06), r (rel) = 0.024323386674344295 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 13.669527953941104 (tol=1e-06), r (rel) = 0.0008318479465019746 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.1007548247168061 (tol=1e-06), r (rel) = 6.131352474148695e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 2.874325629553768e-06 (tol=1e-06), r (rel) = 1.7491473594250489e-10 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 7: r (abs) = 1.6910766425855706e-10 (tol=1e-06), r (rel) = 1.0290908634534707e-14 (tol=1e-10)
Solving for pressure: 1.5 kPa
INFO:scifem.solvers:Newton iteration 1: r (abs) = 18538.97866409216 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 1640.959774677168 (tol=1e-06), r (rel) = 0.08851403329221776 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 112.0549190153566 (tol=1e-06), r (rel) = 0.006044287608593774 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 1.5728254930692902 (tol=1e-06), r (rel) = 8.483884261195412e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.0009263582459760627 (tol=1e-06), r (rel) = 4.9968138092219215e-08 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 3.173674569658301e-10 (tol=1e-06), r (rel) = 1.7118928864216983e-14 (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=True, color="firebrick", label="Inflated")
p.add_legend()
if not pyvista.OFF_SCREEN:
p.show()
else:
p.screenshot(outdir / "pressure.png")
[2025-12-31 10:15:05.609] [info] Checking required entities per dimension
[2025-12-31 10:15:05.610] [info] Cell type: 0 dofmap: 2156x4
[2025-12-31 10:15:05.610] [info] Global index computation
[2025-12-31 10:15:05.610] [info] Got 1 index_maps
[2025-12-31 10:15:05.610] [info] Get global indices
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) = 19465.414596699968 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 6230.433034710484 (tol=1e-06), r (rel) = 0.3200770784387376 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 436.76122175072663 (tol=1e-06), r (rel) = 0.022437807300789377 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 14.537283134610622 (tol=1e-06), r (rel) = 0.0007468262780837544 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.02736998605999696 (tol=1e-06), r (rel) = 1.4060828719588167e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 1.3207282188397285e-07 (tol=1e-06), r (rel) = 6.78499917008516e-12 (tol=1e-10)
Solving for active tension: 5.0 kPa
INFO:scifem.solvers:Newton iteration 1: r (abs) = 16230.883333739763 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 1757.830131265928 (tol=1e-06), r (rel) = 0.10830156899790283 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 67.09529069496023 (tol=1e-06), r (rel) = 0.004133804015181765 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.35642960262452195 (tol=1e-06), r (rel) = 2.1959963317806493e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.2023615776393665e-05 (tol=1e-06), r (rel) = 7.407862855744709e-10 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 1.62505710834386e-10 (tol=1e-06), r (rel) = 1.0012129807906334e-14 (tol=1e-10)
Solving for active tension: 10.0 kPa
INFO:scifem.solvers:Newton iteration 1: r (abs) = 43543.28653398155 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 7008.055214883029 (tol=1e-06), r (rel) = 0.1609445628182862 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 490.05138629583 (tol=1e-06), r (rel) = 0.011254349988336085 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 9.917170067306438 (tol=1e-06), r (rel) = 0.00022775428445363202 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.007696114157507397 (tol=1e-06), r (rel) = 1.7674628559562874e-07 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 7.702044716506979e-09 (tol=1e-06), r (rel) = 1.768824847544808e-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=True, color="red", label="Contracted")
p.add_legend()
if not pyvista.OFF_SCREEN:
p.show()
else:
p.screenshot(outdir / "active.png")