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 io4dolfinx
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,
    )

2026-02-27 16:30:26 [debug    ] Convert file lv_fixedpoint_unloader/geometry/lv_ellipsoid.msh to dolfin

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
    io4dolfinx.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

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[10], line 13
      2 prestress_problem = pulse.unloading.FixedPointUnloader(
      3     model=model,
      4     geometry=geometry,
   (...)      8     unload_parameters={"ramp_steps": 20, "max_iter": 15, "tol": 1e-4},
      9 )
     11 # Execute unloading
     12 # This returns the inverse displacement u_pre (mapping Target -> Reference)
---> 13 u_pre = prestress_problem.unload()
     15 # Save the result
     16 io4dolfinx.write_function_on_input_mesh(
     17     prestress_fname,
     18     u_pre,
     19     time=0.0,
     20     name="u_pre",
     21 )

File /dolfinx-env/lib/python3.12/site-packages/pulse/unloading.py:149, in FixedPointUnloader.unload(self)
    144 logger.info(f"Unloading Iteration {i + 1}/{max_iter}")
    146 # 1. Solve the Forward Problem on the current reference geometry (X_k)
    147 # We create a new StaticProblem instance to ensure forms are initialized
    148 # with the current mesh coordinates.
--> 149 problem = StaticProblem(
    150     model=self.model,
    151     geometry=self.geometry,
    152     bcs=self.bcs,
    153     parameters=self.problem_parameters,
    154 )
    156 for ramp in np.linspace(0.0, 1.0, self.unload_parameters["ramp_steps"]):
    157     for traction, target_pressure, name in self.targets:

File <string>:8, in __init__(self, model, geometry, parameters, bcs, cavities)

File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:85, in StaticProblem.__post_init__(self)
     83 self.parameters = parameters
     84 self._init_spaces()
---> 85 self._init_forms()
     86 logger.debug("Initialized StaticProblem with parameters:")
     87 for key, value in self.parameters.items():

File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:518, in StaticProblem._init_forms(self)
    515 if self.geometry.markers is None:
    516     raise RuntimeError("Missing markers in geometry")
--> 518 R = self.R
    519 u = self.states
    520 K = self.K(R)

File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:410, in StaticProblem.R(self)
    405 @property
    406 def R(self):
    407     # Order is always (u, cavity pressures, rigid body, p)
    409     R = self._empty_form()
--> 410     R_material = self._material_form(self.u, p=self.p)
    411     R_cavity = self._cavity_pressure_form(self.u, self.cavity_pressures)
    412     R_robin = self._robin_form(self.u)

File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:274, in StaticProblem._material_form(self, u, p)
    272 forms = self._empty_form()
    273 # Integrate over reference configuration dX = J_map * dx
--> 274 forms[0] += ufl.inner(self.model.S(C), 0.5 * var_C) * self.geometry.dx
    276 if self.is_incompressible:
    277     forms[-1] += (J - 1.0) * self.p_test * self.geometry.dx

File /dolfinx-env/lib/python3.12/site-packages/pulse/cardiac_model.py:111, in CardiacModel.S(self, C, C_dot)
    104 def S(
    105     self,
    106     C: ufl.core.expr.Expr,
    107     C_dot: ufl.core.expr.Expr | None = None,
    108 ) -> ufl.core.expr.Expr:
    109     """Cauchy stress for the cardiac model."""
--> 111     S = self.material.S(C, dev=True) + self.active.S(C, dev=True) + self.compressibility.S(C)
    112     if C_dot is not None:
    113         S += self.viscoelasticity.S(C_dot)

File /dolfinx-env/lib/python3.12/site-packages/pulse/active_model.py:133, in ActiveModel.S(self, C, dev)
    131 else:
    132     Cdev = C
--> 133 return 2.0 * ufl.diff(self.strain_energy(Cdev), Cdev)

File /dolfinx-env/lib/python3.12/site-packages/ufl/operators.py:382, in diff(f, v)
    380     return VariableDerivative(f, v)
    381 else:
--> 382     raise ValueError("Expecting a Variable or SpatialCoordinate in diff.")

ValueError: Expecting a Variable or SpatialCoordinate in diff.

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)
io4dolfinx.read_function(
    prestress_fname,
    u_pre,
    time=0.0,
    name="u_pre",
)

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

# 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")

# 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)

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)

vtx.close()

print(
    "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