Pre-stressing of a Bi-Ventricular Geometry

Pre-stressing of a Bi-Ventricular Geometry#

In cardiac mechanics simulations, we often start from a geometry acquired from medical imaging (e.g., MRI or CT) at a specific point in the cardiac cycle, typically end-diastole. At this point, the ventricles are pressurized and thus already deformed. However, standard finite element mechanics assumes the initial geometry is stress-free.

To correctly simulate the mechanics, we need to find the unloaded reference configuration corresponding to the acquired geometry and the known end-diastolic pressures. This process is often called “pre-stressing” or “inverse mechanics”.

In this demo, we solve the Inverse Elasticity Problem (IEP) as formulated in [BRR24] for a Bi-Ventricular (BiV) geometry. We formulate the equilibrium equations directly on the known target configuration and solve for the “inverse displacement” that maps points back to the stress-free state.

Mathematical Formulation#

Let \(\Omega_t\) be the known target (loaded) configuration and \(\Omega_0\) be the unknown reference (unloaded) configuration. We seek a mapping \(\boldsymbol{\chi}^{-1}: \Omega_t \to \Omega_0\).

We define the inverse displacement field \(\mathbf{u}\) on \(\Omega_t\) such that:

\[ \mathbf{X} = \mathbf{x} + \mathbf{u}(\mathbf{x}) \]

where \(\mathbf{x} \in \Omega_t\) are the current coordinates and \(\mathbf{X} \in \Omega_0\) are the reference coordinates.

Kinematics#

The inverse deformation gradient \(\mathbf{f}\) is defined as:

\[ \mathbf{f} = \frac{\partial \mathbf{X}}{\partial \mathbf{x}} = \mathbf{I} + \nabla_{\mathbf{x}} \mathbf{u} \]

The physical deformation gradient \(\mathbf{F}\) (mapping reference to target) is the inverse of \(\mathbf{f}\):

\[ \mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \mathbf{f}^{-1} = (\mathbf{I} + \nabla_{\mathbf{x}} \mathbf{u})^{-1} \]

The Jacobian is \(J = \det \mathbf{F} = (\det \mathbf{f})^{-1}\).

Equilibrium#

We solve the balance of linear momentum. The weak form is pulled back from the reference configuration to the target configuration \(\Omega_t\) (where the mesh is defined).

\[ \int_{\Omega_t} \sigma : \nabla_{\mathbf{x}} \mathbf{v} \, dx - \int_{\partial \Omega_t} \mathbf{t} \cdot \mathbf{v} \, ds = 0 \]

Here \(\sigma = J^{-1} \mathbf{P} \mathbf{F}^T\) is the Cauchy stress, and \(\mathbf{P}\) is the First Piola-Kirchhoff stress computed from the material model using \(\mathbf{F}\).

In fenicsx-pulse, the PrestressProblem class automates this specific formulation.

Imports#

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

We set up logging to monitor the process.

comm = MPI.COMM_WORLD
logging.basicConfig(level=logging.INFO)
logging.getLogger("scifem").setLevel(logging.WARNING)

1. Geometry Generation (Target Configuration)#

We generate an Bi-Ventricular (BiV) geometry using cardiac-geometries which represents our target geometry (e.g., the end-diastolic state). This geometry is generated from the mean shape of an atlas from the UK Biobank. We also generate the fiber architecture using a rule-based method.

mode = -1
std = 0
char_length = 10.0

outdir = Path("biv-prestress")
outdir.mkdir(parents=True, exist_ok=True)
geodir = outdir / "geometry"
if not geodir.exists():
    comm.barrier()
    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:Generating points from /github/home/.ukb/UKBRVLV.h5 using mode -1 and std 0.0

INFO:ukb.surface:Saved biv-prestress/geometry/EPI_ED.stl

INFO:ukb.surface:Saved biv-prestress/geometry/MV_ED.stl

INFO:ukb.surface:Saved biv-prestress/geometry/AV_ED.stl

INFO:ukb.surface:Saved biv-prestress/geometry/TV_ED.stl

INFO:ukb.surface:Saved biv-prestress/geometry/PV_ED.stl

INFO:ukb.surface:Saved biv-prestress/geometry/LV_ED.stl

INFO:ukb.surface:Saved biv-prestress/geometry/RV_ED.stl

INFO:ukb.surface:Saved biv-prestress/geometry/RVFW_ED.stl

INFO:ukb.clip:Folder: biv-prestress/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 biv-prestress/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 biv-prestress/geometry/lv_clipped.ply

INFO:ukb.clip:Reading biv-prestress/geometry/RV_ED.stl

INFO:ukb.clip:Reading biv-prestress/geometry/RVFW_ED.stl

INFO:ukb.clip:Merging RV and RVFW

INFO:ukb.clip:Smoothing RV

INFO:ukb.clip:Saving biv-prestress/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 biv-prestress/geometry/EPI_ED.stl

INFO:ukb.clip:Saving biv-prestress/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

INFO:ukb.mesh:Created mesh biv-prestress/geometry/ED_clipped.msh

2025-12-31 10:19:23 [debug    ] Convert file biv-prestress/geometry/ED_clipped.msh to dolfin

Info    : Reading 'biv-prestress/geometry/ED_clipped.msh'...
Info    : 11 entities
Info    : 690 nodes
Info    : 3406 elements
Info    : 3 parametrizations
Info    : [  0%] Processing parametrizations                                                                                
Info    : [ 10%] Processing parametrizations                                                                                
Info    : [ 40%] Processing parametrizations                                                                                
Info    : Done reading 'biv-prestress/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: 690

INFO:ldrb.ldrb:  Num cells: 2130

INFO:ldrb.ldrb:  Apex coord: (45.87, -20.08, -29.44)

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 biv-prestress/geometry

Load the geometry and convert it to pulse.HeartGeometry.

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

INFO:cardiac_geometries.geometry:Reading geometry from biv-prestress/geometry

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

2. Constitutive Model#

We define the material properties. For this example, we use the Usyk model [ULM02], which is a transversely isotropic hyperelastic model.

We use an incompressible formulation.

3. Boundary Conditions#

We apply the loading and constraints that define the target state.

Neumann (Pressure): The target end-diastolic pressure is applied to the endocardial surfaces. We apply \(P_{LV} = 2000\) Pa to the LV and assume \(P_{RV} = 0.5 P_{LV}\).
Robin (Springs): To prevent rigid body motion and mimic the pericardial constraint, we apply spring conditions to the epicardium and the base.

def setup_problem(geometry, f0, s0, target_pressure=2000.0):

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

    comp = pulse.compressibility.Incompressible()

    model = pulse.CardiacModel(
        material=material,
        compressibility=comp,
        active=pulse.active_model.Passive(),
    )
    pressure_lv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, 0.0), "Pa")
    pressure_rv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, 0.0), "Pa")

    # Epicardial and Basal springs
    alpha_epi = pulse.Variable(
        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
        "Pa / m",
    )
    robin_epi = pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
    alpha_base = pulse.Variable(
        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(2e5)),
        "Pa / m",
    )
    robin_base = pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])


    # Endocardial Pressures
    neumann_lv = pulse.NeumannBC(traction=pressure_lv, marker=geometry.markers["LV"][0])
    neumann_rv = pulse.NeumannBC(traction=pressure_rv, marker=geometry.markers["RV"][0])

    bcs = pulse.BoundaryConditions(neumann=(neumann_lv, neumann_rv), robin=(robin_epi, robin_base))
    return model, bcs, pressure_lv, pressure_rv, target_pressure

model, bcs, pressure_lv, pressure_rv, target_pressure = setup_problem(geometry, geo.f0, geo.s0)

4. Solving the Inverse Elasticity Problem (IEP)#

The PrestressProblem class sets up the inverse formulation.

We solve the problem incrementally using continuation (load stepping). We ramp the pressures from 0 to the target pressures.

Note: The geometry geometry.mesh does not change during this loop. It remains the target configuration. The solution prestress_problem.u is updated to reflect the inverse displacement required to map from this fixed target back to the corresponding reference state for the current pressure load.

prestress_fname = outdir / "prestress_biv_inverse.bp"
if not prestress_fname.exists():
    prestress_problem = pulse.unloading.PrestressProblem(
        geometry=geometry,
        model=model,
        bcs=bcs,
        parameters={"u_space": "P_2"},
        targets=[
            pulse.unloading.TargetPressure(traction=pressure_lv, target=target_pressure, name="LV"),
            pulse.unloading.TargetPressure(traction=pressure_rv, target=target_pressure * 0.5, name="RV"),
        ],
        ramp_steps=20,
    )

    u_pre = prestress_problem.unload()

    # We save the computed inverse displacement field $\mathbf{u}_{pre}$.
    # This field maps: **Target Geometry** ($\mathbf{x}$) $\to$ **Reference Geometry** ($\mathbf{X}$).
    adios4dolfinx.write_function_on_input_mesh(
        prestress_fname,
        u_pre,
        time=0.0,
        name="u_pre",
    )

    with dolfinx.io.VTXWriter(
        comm, outdir / "prestress_biv_backward.bp", [u_pre], engine="BP4",
    ) as vtx:
        vtx.write(0.0)

    # We can visualize the inverse displacement field.
    try:
        import pyvista
    except ImportError:
        print("Pyvista is not installed")
    else:
        # Create plotter and pyvista grid
        p = pyvista.Plotter()

        topology, cell_types, vtk_geometry = dolfinx.plot.vtk_mesh(u_pre.function_space)
        grid = pyvista.UnstructuredGrid(topology, cell_types, vtk_geometry)

        # Attach vector values to grid
        grid["u"] = u_pre.x.array.reshape((vtk_geometry.shape[0], 3))
        actor_0 = p.add_mesh(grid, style="wireframe", color="k")

        # Warp the mesh by the inverse displacement vector to visualize the unloaded Reference Configuration
        warped = grid.warp_by_vector("u", factor=1.0)
        actor_1 = p.add_mesh(warped, color="red", opacity=0.8, label="Recovered Reference")

        p.add_legend()
        p.show_axes()
        if not pyvista.OFF_SCREEN:
            p.show()
        else:
            figure_as_array = p.screenshot("prestress_inverse_displacement.png")

INFO:pulse.unloading:Ramping LV traction to 0.0000

INFO:pulse.unloading:Ramping RV traction to 0.0000

INFO:pulse.unloading:Ramping LV traction to 105.2632

INFO:pulse.unloading:Ramping RV traction to 52.6316

INFO:pulse.unloading:Ramping LV traction to 210.5263

INFO:pulse.unloading:Ramping RV traction to 105.2632

INFO:pulse.unloading:Ramping LV traction to 315.7895

INFO:pulse.unloading:Ramping RV traction to 157.8947

INFO:pulse.unloading:Ramping LV traction to 421.0526

INFO:pulse.unloading:Ramping RV traction to 210.5263

INFO:pulse.unloading:Ramping LV traction to 526.3158

INFO:pulse.unloading:Ramping RV traction to 263.1579

INFO:pulse.unloading:Ramping LV traction to 631.5789

INFO:pulse.unloading:Ramping RV traction to 315.7895

INFO:pulse.unloading:Ramping LV traction to 736.8421

INFO:pulse.unloading:Ramping RV traction to 368.4211

INFO:pulse.unloading:Ramping LV traction to 842.1053

INFO:pulse.unloading:Ramping RV traction to 421.0526

INFO:pulse.unloading:Ramping LV traction to 947.3684

INFO:pulse.unloading:Ramping RV traction to 473.6842

INFO:pulse.unloading:Ramping LV traction to 1052.6316

INFO:pulse.unloading:Ramping RV traction to 526.3158

INFO:pulse.unloading:Ramping LV traction to 1157.8947

INFO:pulse.unloading:Ramping RV traction to 578.9474

INFO:pulse.unloading:Ramping LV traction to 1263.1579

INFO:pulse.unloading:Ramping RV traction to 631.5789

INFO:pulse.unloading:Ramping LV traction to 1368.4211

INFO:pulse.unloading:Ramping RV traction to 684.2105

INFO:pulse.unloading:Ramping LV traction to 1473.6842

INFO:pulse.unloading:Ramping RV traction to 736.8421

INFO:pulse.unloading:Ramping LV traction to 1578.9474

INFO:pulse.unloading:Ramping RV traction to 789.4737

INFO:pulse.unloading:Ramping LV traction to 1684.2105

INFO:pulse.unloading:Ramping RV traction to 842.1053

INFO:pulse.unloading:Ramping LV traction to 1789.4737

INFO:pulse.unloading:Ramping RV traction to 894.7368

INFO:pulse.unloading:Ramping LV traction to 1894.7368

INFO:pulse.unloading:Ramping RV traction to 947.3684

INFO:pulse.unloading:Ramping LV traction to 2000.0000

INFO:pulse.unloading:Ramping RV traction to 1000.0000

2025-12-31 10:20:39.402 (  76.657s) [    7FBDE0694140]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

V = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 2, (3,)))
u_pre = dolfinx.fem.Function(V)
adios4dolfinx.read_function(
    prestress_fname,
    u_pre,
    time=0.0,
    name="u_pre",
)

5. Verification (Forward Problem)#

To verify the result, we perform a full explicit deformation to the reference configuration.

Deform Mesh: We apply the inverse displacement \(\mathbf{u}_{pre}\) to the mesh nodes: \(\mathbf{X} = \mathbf{x} + \mathbf{u}_{pre}\). The mesh now represents the Reference (Unloaded) configuration.
Map Fibers: The fiber fields must be mapped (pulled back) to this new reference configuration to ensure correct material orientation.
Forward Solve: We assume this new geometry is stress-free (pressure=0). We then ramp the pressure back up to the target value.

If the prestressing was successful, the deformed geometry at full load should match the original target mesh.

print("\nDeforming mesh to recovered reference configuration...")
geometry.deform(u_pre)

Deforming mesh to recovered reference configuration...

# Map fiber fields to the new configuration
print("Mapping fiber fields...")
f0 = pulse.utils.map_vector_field(f=geo.f0, u=u_pre, normalize=True, name="f0_unloaded")
s0 = pulse.utils.map_vector_field(f=geo.s0, u=u_pre, normalize=True, name="s0_unloaded")

Mapping fiber fields...

# Setup the problem on the REFERENCE geometry with MAPPED fibers
model_unloaded, bcs_unloaded, pressure_lv_unloaded, pressure_biv_unloaded, target_pressure_unloaded = setup_problem(geometry, f0, s0)

# Initialize Forward Problem using the UNLOADED model and BCs
forward_problem = pulse.StaticProblem(
    model=model_unloaded,
    geometry=geometry,
    bcs=bcs_unloaded,
    parameters={"u_space": "P_2"},
)

import shutil
shutil.rmtree(outdir / "prestress_biv.bp", ignore_errors=True)
vtx = dolfinx.io.VTXWriter(comm, outdir / "prestress_biv.bp", [forward_problem.u], engine="BP4")

# Step 1: Initial State (Zero Pressure)
# Since the mesh is now the Reference configuration, zero pressure means zero displacement.
print("\nSolving forward problem: Initial state (Reference)...")
pressure_lv_unloaded.assign(0.0)
pressure_biv_unloaded.assign(0.0)
forward_problem.solve()
vtx.write(0.0)

Solving forward problem: Initial state (Reference)...

# Step 2: Reload to Target Pressure
# We apply the original target pressures. The mesh should deform from Reference -> Target.
print("Solving forward problem: Reloading to target pressure...")
ramp_steps = 20
for ramp in np.linspace(0.0, 1.0, ramp_steps):
    current_p = target_pressure_unloaded * ramp
    pressure_lv_unloaded.assign(current_p)
    pressure_biv_unloaded.assign(current_p * 0.5)

    print(f"Solving for pressure fraction: {ramp:.2f} (P_LV = {current_p:.2f} Pa)")
    forward_problem.solve()

    vtx.write(ramp * ramp_steps + 1)

Solving forward problem: Reloading to target pressure...
Solving for pressure fraction: 0.00 (P_LV = 0.00 Pa)
Solving for pressure fraction: 0.05 (P_LV = 105.26 Pa)

Solving for pressure fraction: 0.11 (P_LV = 210.53 Pa)

Solving for pressure fraction: 0.16 (P_LV = 315.79 Pa)

Solving for pressure fraction: 0.21 (P_LV = 421.05 Pa)

Solving for pressure fraction: 0.26 (P_LV = 526.32 Pa)

Solving for pressure fraction: 0.32 (P_LV = 631.58 Pa)

Solving for pressure fraction: 0.37 (P_LV = 736.84 Pa)

Solving for pressure fraction: 0.42 (P_LV = 842.11 Pa)

Solving for pressure fraction: 0.47 (P_LV = 947.37 Pa)

Solving for pressure fraction: 0.53 (P_LV = 1052.63 Pa)

Solving for pressure fraction: 0.58 (P_LV = 1157.89 Pa)

Solving for pressure fraction: 0.63 (P_LV = 1263.16 Pa)

Solving for pressure fraction: 0.68 (P_LV = 1368.42 Pa)

Solving for pressure fraction: 0.74 (P_LV = 1473.68 Pa)

Solving for pressure fraction: 0.79 (P_LV = 1578.95 Pa)

Solving for pressure fraction: 0.84 (P_LV = 1684.21 Pa)

Solving for pressure fraction: 0.89 (P_LV = 1789.47 Pa)

Solving for pressure fraction: 0.95 (P_LV = 1894.74 Pa)

Solving for pressure fraction: 1.00 (P_LV = 2000.00 Pa)

vtx.close()

# Save final state. This displacement maps **Reference** -> **Recovered Target**.
with dolfinx.io.VTXWriter(comm, outdir / "prestress_biv_forward.bp", [forward_problem.u], engine="BP4") as vtx:
    vtx.write(0.0)

print("Done. You can now verify that the final geometry matches the original target geometry in Paraview.")

Done. You can now verify that the final geometry matches the original target geometry in Paraview.

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    # Create plotter
    p = pyvista.Plotter()

    topology, cell_types, vtk_geometry = dolfinx.plot.vtk_mesh(forward_problem.u_space)
    grid = pyvista.UnstructuredGrid(topology, cell_types, vtk_geometry)

    # Attach the final displacement result
    grid["u"] = forward_problem.u.x.array.reshape((vtk_geometry.shape[0], 3))

    # Reference Configuration (Wireframe) - This is the mesh we started the forward solve with
    actor_0 = p.add_mesh(grid, style="wireframe", color="k", label="Reference Config")

    # Recovered Target (Deformed by final u)
    # Ideally, we should also load the original target mesh for comparison, but here we show the deformed result.
    warped = grid.warp_by_vector("u", factor=1.0)
    actor_1 = p.add_mesh(warped, color="blue", opacity=0.5, label="Recovered Target")

    p.add_legend()
    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        figure_as_array = p.screenshot("prestress_forward_displacement.png")

References#

[BRR24]

NA Barnafi, Francesco Regazzoni, and Davide Riccobelli. Reconstructing relaxed configurations in elastic bodies: mathematical formulations and numerical methods for cardiac modeling. Computer Methods in Applied Mechanics and Engineering, 423:116845, 2024.

[ULM02]

Taras P Usyk, Ian J LeGrice, and Andrew D McCulloch. Computational model of three-dimensional cardiac electromechanics. Computing and visualization in science, 4(4):249–257, 2002.