Pre-stressing of a Left Ventricle Geometry using Fixed-Point Unloading

Pre-stressing of a Left Ventricle Geometry using Fixed-Point Unloading#

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 ventricle is 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 pressure. This process is often called “pre-stressing” or “inverse mechanics”.

In this demo, we solve the inverse problem using a Fixed-Point Iteration (also known as the Backward Displacement Method o Sellier’s method) [Sel11]. Unlike the Inverse Elasticity Problem (IEP) which formulates equilibrium on the target configuration, this method iteratively updates the reference coordinates \(\mathbf{X}\) by subtracting the displacement \(\mathbf{u}\) computed from a forward solve.

Mathematical Formulation#

Let \(\Omega_t\) be the known target (loaded) configuration with coordinates \(\mathbf{x}_{target}\). We seek the reference (unloaded) configuration \(\Omega_0\) with coordinates \(\mathbf{X}\).

The algorithm proceeds as follows:

  1. Initialize reference geometry: \(\mathbf{X}_0 = \mathbf{x}_{target}\).

  2. For iteration \(k=0, 1, \dots\): a. Solve the Forward mechanics problem on the geometry defined by \(\mathbf{X}_k\) to get displacement \(\mathbf{u}_k\). b. Update the reference geometry: \( \mathbf{X}_{k+1} = \mathbf{x}_{target} - \mathbf{u}_k \) c. Check convergence: \(||\mathbf{X}_{k+1} - \mathbf{X}_k|| < \text{tol}\).

In fenicsx-pulse, the FixedPointUnloader class automates this iterative process.

Imports#

from pathlib import Path
import math
from mpi4py import MPI
import dolfinx
import logging
import numpy as np
import pulse
import pulse.unloading
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 idealized Left Ventricular (LV) geometry using cardiac-geometries which represents our target geometry (e.g., the end-diastolic state). We also generate the fiber architecture.

outdir = Path("lv_fixedpoint_unloader")
outdir.mkdir(parents=True, exist_ok=True)
geodir = outdir / "geometry"
if not geodir.exists():
    comm.barrier()
    cardiac_geometries.mesh.lv_ellipsoid(
        outdir=geodir,
        create_fibers=True,
        fiber_space="Quadrature_6",
        r_short_endo=0.025,
        r_short_epi=0.035,
        r_long_endo=0.09,
        r_long_epi=0.097,
        psize_ref=0.03,
        mu_apex_endo=-math.pi,
        mu_base_endo=-math.acos(5 / 17),
        mu_apex_epi=-math.pi,
        mu_base_epi=-math.acos(5 / 20),
        comm=comm,
        fiber_angle_epi=-60,
        fiber_angle_endo=60,
    )

2025-12-31 10:21:52 [debug    ] Convert file lv_fixedpoint_unloader/geometry/lv_ellipsoid.msh to dolfin

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

Info    : Reading 'lv_fixedpoint_unloader/geometry/lv_ellipsoid.msh'...
Info    : 54 entities
Info    : 191 nodes
Info    : 990 elements
Info    : Done reading 'lv_fixedpoint_unloader/geometry/lv_ellipsoid.msh'

# 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 lv_fixedpoint_unloader/geometry

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

2. Constitutive Model & Boundary Conditions#

We define a helper function to setup the mechanical model.

  • Material: Usyk model (transversely isotropic).

  • Compressibility: Compressible formulation with high bulk modulus.

  • BCs:

    • Neumann: Endocardial pressure (Target = 2000 Pa).

    • Robin: Epicardial and Basal springs to mimic tissue support.

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

    material = pulse.material_models.Usyk(f0=f0, s0=s0, n0=n0)
    comp = pulse.compressibility.Compressible3(kappa=pulse.Variable(5e4, "Pa"))

    model = pulse.CardiacModel(
        material=material,
        compressibility=comp,
        active=pulse.active_model.Passive(),
    )
    pressure = 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(2e5)),
        "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 Pressure
    neumann = pulse.NeumannBC(traction=pressure, marker=geometry.markers["ENDO"][0])

    bcs = pulse.BoundaryConditions(neumann=(neumann,), robin=(robin_epi, robin_base))
    return model, bcs, pressure, target_pressure

Initialize model with fibers from the target geometry

model, bcs, pressure, target_pressure = setup_problem(geometry, geo.f0, geo.s0, geo.n0)

# Define the loading target for the unloader
target_pressure_obj = pulse.unloading.TargetPressure(traction=pressure, target=target_pressure)

3. Solving the Inverse Problem (Fixed Point Iteration)#

The FixedPointUnloader automates the iterative process of finding the stress-free reference configuration.

  • ramp_steps: Number of load steps for the forward solver within each iteration (improves stability).

  • max_iter: Maximum number of fixed-point iterations.

  • tol: Geometric convergence tolerance.

prestress_fname = outdir / "prestress_lv_inverse.bp"

if not prestress_fname.exists():
    prestress_problem = pulse.unloading.FixedPointUnloader(
        model=model,
        geometry=geometry,
        bcs=bcs,
        problem_parameters={"u_space": "P_2"},
        targets=[target_pressure_obj],
        unload_parameters={"ramp_steps": 20, "max_iter": 15, "tol": 1e-4},
    )

    # Execute unloading
    # This returns the inverse displacement u_pre (mapping Target -> Reference)
    u_pre = prestress_problem.unload()

    # Save the result
    adios4dolfinx.write_function_on_input_mesh(
        prestress_fname,
        u_pre,
        time=0.0,
        name="u_pre",
    )

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

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

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

        grid["u"] = u_pre.x.array.reshape((vtk_geometry.shape[0], 3))
        actor_0 = p.add_mesh(grid, style="wireframe", color="k", label="Target (Loaded)")

        # Warp by u_pre to show the recovered Reference (Unloaded) configuration
        warped = grid.warp_by_vector("u", factor=1.0)
        actor_1 = p.add_mesh(warped, color="red", opacity=0.8, label="Reference (Unloaded)")

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

INFO:pulse.unloading:Starting FixedPointUnloader...

INFO:pulse.unloading:Unloading Iteration 1/15

INFO:pulse.unloading:  Geometric Error (rel): 2.04e-01 (abs: 1.48e-01)

INFO:pulse.unloading:Unloading Iteration 2/15

INFO:pulse.unloading:  Geometric Error (rel): 7.32e-02 (abs: 5.32e-02)

INFO:pulse.unloading:Unloading Iteration 3/15

INFO:pulse.unloading:  Geometric Error (rel): 8.68e-03 (abs: 6.31e-03)

INFO:pulse.unloading:Unloading Iteration 4/15

INFO:pulse.unloading:  Geometric Error (rel): 5.76e-03 (abs: 4.18e-03)

INFO:pulse.unloading:Unloading Iteration 5/15

INFO:pulse.unloading:  Geometric Error (rel): 2.79e-04 (abs: 2.03e-04)

INFO:pulse.unloading:Unloading Iteration 6/15

INFO:pulse.unloading:  Geometric Error (rel): 2.34e-04 (abs: 1.70e-04)

INFO:pulse.unloading:Unloading Iteration 7/15

INFO:pulse.unloading:  Geometric Error (rel): 4.48e-05 (abs: 3.25e-05)

INFO:pulse.unloading:FixedPointUnloader converged in 7 iterations.

2025-12-31 10:23:46.908 (   0.877s) [    7FD99616C140]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

4. Verification (Forward Problem)#

To verify the result, we perform a explicit deformation to the reference configuration and solve the forward problem.

  1. Deform Mesh: Apply u_pre to the mesh nodes. The mesh now represents the Reference configuration.

  2. Map Fibers: Pull back the fiber fields to the reference configuration.

  3. Forward Solve: Ramp pressure from 0 to target_pressure.

If successful, the final deformed geometry should match the original target mesh.

# Reload u_pre
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",
)

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

Deforming mesh to recovered reference configuration...

# Map fiber fields
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")
n0 = pulse.utils.map_vector_field(f=geo.n0, u=u_pre, normalize=True, name="n0_unloaded")

Mapping fiber fields...

# Setup Forward Problem on Reference Mesh
model_unloaded, bcs_unloaded, pressure_unloaded, target_pressure_unloaded = setup_problem(geometry, f0, s0, n0)

forward_problem = pulse.StaticProblem(
    model=model_unloaded,
    geometry=geometry,
    bcs=bcs_unloaded,
    parameters={"u_space": "P_2"},
)

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

print("\nSolving forward problem: Initial state (Reference)...")
pressure_unloaded.assign(0.0)
forward_problem.solve()
vtx.write(0.0)

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

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_unloaded.assign(current_p)

    print(f"Solving for pressure fraction: {ramp:.2f} (P = {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 = 0.00 Pa)
Solving for pressure fraction: 0.05 (P = 105.26 Pa)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

vtx.close()

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

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

# Visualization of Forward Result
try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    p = pyvista.Plotter()
    topology, cell_types, vtk_geometry = dolfinx.plot.vtk_mesh(forward_problem.u_space)
    grid = pyvista.UnstructuredGrid(topology, cell_types, vtk_geometry)

    grid["u"] = forward_problem.u.x.array.reshape((vtk_geometry.shape[0], 3))

    # Reference (Wireframe)
    actor_0 = p.add_mesh(grid, style="wireframe", color="k", label="Reference Config")

    # Recovered Target (Deformed)
    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("lv_prestress_forward_displacement.png")

References#

[Sel11]

M. Sellier. An iterative method for the inverse elasto-static problem. Journal of Fluids and Structures, 27(8):1461–1470, 2011. URL: https://www.sciencedirect.com/science/article/pii/S088997461100123X, doi:https://doi.org/10.1016/j.jfluidstructs.2011.08.002.

Contents