Pre-stressing of a Left Ventricle Ellipsoid

Pre-stressing of a Left Ventricle Ellipsoid#

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 Elasticity Problem (IEP) as formulated in [BRR24]. 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 math
import numpy as np
import pulse
import pulse.prestress
import cardiac_geometries
import cardiac_geometries.geometry

We set up logging to monitor the process.

comm = MPI.COMM_WORLD
logging.basicConfig(level=logging.INFO)
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

1. Geometry Generation (Target Configuration)#

We generate an idealized LV ellipsoid which represents our target geometry (e.g., the end-diastolic state). We also generate the fiber architecture.

geodir = Path("lv_ellipsoid-prestress")
if not geodir.exists():
    comm.barrier()
    cardiac_geometries.mesh.lv_ellipsoid(
        outdir=geodir,
        create_fibers=True,
        fiber_space="P_2",
        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-15 07:18:15 [debug    ] Convert file lv_ellipsoid-prestress/lv_ellipsoid.msh to dolfin

INFO:cardiac_geometries.geometry:Reading geometry from lv_ellipsoid-prestress

Info    : Reading 'lv_ellipsoid-prestress/lv_ellipsoid.msh'...
Info    : 54 entities
Info    : 191 nodes
Info    : 990 elements
Info    : Done reading 'lv_ellipsoid-prestress/lv_ellipsoid.msh'
[2025-12-15 07:18:15.097] [info] Extract basic topology: 2828->2828
[2025-12-15 07:18:15.097] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:18:15.097] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:18:15.097] [info] GPS pseudo-diameter:(28) 58-440
[2025-12-15 07:18:15.097] [info] Create topology (generalised)
[2025-12-15 07:18:15.098] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.098] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.098] [info] Compute ghost indices
[2025-12-15 07:18:15.098] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.098] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.098] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.098] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.098] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.098] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.098] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:18:15.098] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.098] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.098] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:18:15.098] [info] Create Geometry (multiple)
[2025-12-15 07:18:15.098] [info] Got 1 dof layouts
[2025-12-15 07:18:15.098] [info] Checking required entities per dimension
[2025-12-15 07:18:15.098] [info] Cell type: 0 dofmap: 707x4
[2025-12-15 07:18:15.098] [info] Global index computation
[2025-12-15 07:18:15.098] [info] Got 1 index_maps
[2025-12-15 07:18:15.098] [info] Get global indices
[2025-12-15 07:18:15.098] [info] Calling compute_local_to_global
[2025-12-15 07:18:15.098] [info] xdofs.size = 2828
[2025-12-15 07:18:15.098] [info] dofmap sizes = 2828 
[2025-12-15 07:18:15.098] [info] all_dofmaps.size = 2828
[2025-12-15 07:18:15.098] [info] nodes.size = 191
[2025-12-15 07:18:15.098] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:18:15.098] [info] XDMF distribute entity data
[2025-12-15 07:18:15.098] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.098] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.098] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.098] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.098] [info] XDMF build map
[2025-12-15 07:18:15.099] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:18:15.099] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.099] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.099] [info] XDMF distribute entity data
[2025-12-15 07:18:15.099] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.099] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.099] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.099] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.099] [info] XDMF build map
[2025-12-15 07:18:15.099] [info] Computing mesh entities of dimension 2
[2025-12-15 07:18:15.099] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.099] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.100] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:15.100] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:18:15.100] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.100] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.100] [info] XDMF distribute entity data
[2025-12-15 07:18:15.100] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.100] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.100] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.100] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.100] [info] XDMF build map
[2025-12-15 07:18:15.100] [info] Computing mesh entities of dimension 1
[2025-12-15 07:18:15.100] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.100] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.101] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:18:15.101] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:18:15.101] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.101] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.101] [info] XDMF distribute entity data
[2025-12-15 07:18:15.101] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.101] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.101] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.101] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.101] [info] XDMF build map
[2025-12-15 07:18:15.101] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:18:15.101] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:18:15.101] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.101] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.102] [info] Opened HDF5 file with id "72057594037927936"
[2025-12-15 07:18:15.102] [info] Adding mesh to node "/Xdmf/Domain"
[2025-12-15 07:18:15.102] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:18:15.103] [info] Adding geometry data to node "/Xdmf/Domain/Grid"
[2025-12-15 07:18:15.103] [info] XDMF: add meshtags (Cell tags)
[2025-12-15 07:18:15.103] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:18:15.103] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:15.103] [info] XDMF: add meshtags (Facet tags)
[2025-12-15 07:18:15.103] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:18:15.103] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:18:15.103] [info] XDMF: add meshtags (Edge tags)
[2025-12-15 07:18:15.103] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:18:15.104] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:18:15.104] [info] XDMF: add meshtags (Vertex tags)
[2025-12-15 07:18:15.104] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:18:15.105] [info] Checking required entities per dimension
[2025-12-15 07:18:15.105] [info] Cell type: 0 dofmap: 707x4
[2025-12-15 07:18:15.105] [info] Global index computation
[2025-12-15 07:18:15.105] [info] Got 1 index_maps
[2025-12-15 07:18:15.105] [info] Get global indices
[2025-12-15 07:18:15.111] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:18:15.113] [info] Checking required entities per dimension
[2025-12-15 07:18:15.113] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.113] [info] Global index computation
[2025-12-15 07:18:15.113] [info] Got 2 index_maps
[2025-12-15 07:18:15.113] [info] Get global indices
[2025-12-15 07:18:15.119] [info] Checking required entities per dimension
[2025-12-15 07:18:15.119] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.119] [info] Global index computation
[2025-12-15 07:18:15.119] [info] Got 2 index_maps
[2025-12-15 07:18:15.119] [info] Get global indices
[2025-12-15 07:18:15.121] [info] Compute face permutations
[2025-12-15 07:18:15.121] [info] Computing permutations for face type 0
[2025-12-15 07:18:15.121] [info] Compute edge permutations
[2025-12-15 07:18:15.127] [info] Opened HDF5 file with id "72057594037927937"
[2025-12-15 07:18:15.127] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-12-15 07:18:15.127] [info] HDF5 Read data rate: 1462.6325316783034 MB/s
[2025-12-15 07:18:15.127] [info] IO permuting cells
[2025-12-15 07:18:15.127] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-12-15 07:18:15.127] [info] HDF5 Read data rate: 574.8683220466517 MB/s
[2025-12-15 07:18:15.128] [info] Using partitioner with cell data (1 cell types)
[2025-12-15 07:18:15.128] [info] Compute partition of cells across ranks
[2025-12-15 07:18:15.128] [info] Building mesh dual graph
[2025-12-15 07:18:15.128] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:18:15.128] [info] Build nonlocal part of mesh dual graph
[2025-12-15 07:18:15.128] [info] Graph edges (local: 2590, non-local: 0)
[2025-12-15 07:18:15.128] [info] Compute graph partition using PT-SCOTCH
[2025-12-15 07:18:15.128] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.128] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.128] [info] Extract basic topology: 2828->2828
[2025-12-15 07:18:15.128] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:18:15.128] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:18:15.129] [info] GPS pseudo-diameter:(28) 706-0
[2025-12-15 07:18:15.129] [info] Create topology (generalised)
[2025-12-15 07:18:15.129] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.129] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.129] [info] Compute ghost indices
[2025-12-15 07:18:15.129] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.129] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.129] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.129] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.129] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.129] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.129] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:18:15.129] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.129] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.129] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:18:15.129] [info] Create Geometry (multiple)
[2025-12-15 07:18:15.129] [info] Got 1 dof layouts
[2025-12-15 07:18:15.129] [info] Checking required entities per dimension
[2025-12-15 07:18:15.129] [info] Cell type: 0 dofmap: 707x4
[2025-12-15 07:18:15.129] [info] Global index computation
[2025-12-15 07:18:15.129] [info] Got 1 index_maps
[2025-12-15 07:18:15.129] [info] Get global indices
[2025-12-15 07:18:15.129] [info] Calling compute_local_to_global
[2025-12-15 07:18:15.129] [info] xdofs.size = 2828
[2025-12-15 07:18:15.129] [info] dofmap sizes = 2828 
[2025-12-15 07:18:15.129] [info] all_dofmaps.size = 2828
[2025-12-15 07:18:15.129] [info] nodes.size = 191
[2025-12-15 07:18:15.129] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:18:15.129] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:18:15.130] [info] XDMF read meshtags (Cell tags)
[2025-12-15 07:18:15.130] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-12-15 07:18:15.130] [info] HDF5 Read data rate: 1731.6494450822809 MB/s
[2025-12-15 07:18:15.130] [info] IO permuting cells
[2025-12-15 07:18:15.130] [info] HDF5 Read data rate: 385.65389335878905 MB/s
[2025-12-15 07:18:15.130] [info] IO permuting cells
[2025-12-15 07:18:15.130] [info] XDMF distribute entity data
[2025-12-15 07:18:15.130] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.130] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.130] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.130] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.130] [info] XDMF build map
[2025-12-15 07:18:15.130] [info] XDMF create meshtags
[2025-12-15 07:18:15.130] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.130] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.130] [info] Computing mesh entities of dimension 2
[2025-12-15 07:18:15.131] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.131] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.131] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:15.131] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:18:15.131] [info] XDMF read meshtags (Facet tags)
[2025-12-15 07:18:15.131] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-12-15 07:18:15.131] [info] HDF5 Read data rate: 577.6114875113764 MB/s
[2025-12-15 07:18:15.131] [info] IO permuting cells
[2025-12-15 07:18:15.131] [info] HDF5 Read data rate: 159.17070723959205 MB/s
[2025-12-15 07:18:15.131] [info] IO permuting cells
[2025-12-15 07:18:15.131] [info] XDMF distribute entity data
[2025-12-15 07:18:15.131] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.131] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.131] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.131] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.131] [info] XDMF build map
[2025-12-15 07:18:15.131] [info] XDMF create meshtags
[2025-12-15 07:18:15.131] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.131] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.131] [info] Computing mesh entities of dimension 1
[2025-12-15 07:18:15.132] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.132] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.132] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:18:15.132] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:18:15.132] [info] XDMF read meshtags (Edge tags)
[2025-12-15 07:18:15.132] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-12-15 07:18:15.132] [info] HDF5 Read data rate: 78.83579695198809 MB/s
[2025-12-15 07:18:15.132] [info] IO permuting cells
[2025-12-15 07:18:15.132] [info] HDF5 Read data rate: 29.44701249785996 MB/s
[2025-12-15 07:18:15.132] [info] IO permuting cells
[2025-12-15 07:18:15.132] [info] XDMF distribute entity data
[2025-12-15 07:18:15.132] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.132] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.132] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.132] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.132] [info] XDMF build map
[2025-12-15 07:18:15.132] [info] XDMF create meshtags
[2025-12-15 07:18:15.132] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.132] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.132] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:18:15.132] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:18:15.132] [info] XDMF read meshtags (Vertex tags)
[2025-12-15 07:18:15.132] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-12-15 07:18:15.133] [info] HDF5 Read data rate: 1.9382192610539069 MB/s
[2025-12-15 07:18:15.133] [info] IO permuting cells
[2025-12-15 07:18:15.133] [info] HDF5 Read data rate: 1.5355086372360844 MB/s
[2025-12-15 07:18:15.133] [info] IO permuting cells
[2025-12-15 07:18:15.133] [info] XDMF distribute entity data
[2025-12-15 07:18:15.133] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.133] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.133] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.133] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.133] [info] XDMF build map
[2025-12-15 07:18:15.133] [info] XDMF create meshtags
[2025-12-15 07:18:15.133] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.133] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.134] [info] Checking required entities per dimension
[2025-12-15 07:18:15.134] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.134] [info] Global index computation
[2025-12-15 07:18:15.134] [info] Got 2 index_maps
[2025-12-15 07:18:15.134] [info] Get global indices
[2025-12-15 07:18:15.135] [info] Compute face permutations
[2025-12-15 07:18:15.135] [info] Computing permutations for face type 0
[2025-12-15 07:18:15.135] [info] Compute edge permutations
[2025-12-15 07:18:15.139] [info] Checking required entities per dimension
[2025-12-15 07:18:15.139] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.139] [info] Global index computation
[2025-12-15 07:18:15.139] [info] Got 2 index_maps
[2025-12-15 07:18:15.140] [info] Get global indices
[2025-12-15 07:18:15.144] [info] Checking required entities per dimension
[2025-12-15 07:18:15.144] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.144] [info] Global index computation
[2025-12-15 07:18:15.144] [info] Got 2 index_maps
[2025-12-15 07:18:15.144] [info] Get global indices

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_ellipsoid-prestress

[2025-12-15 07:18:15.153] [info] Opened HDF5 file with id "72057594037927938"
[2025-12-15 07:18:15.153] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-12-15 07:18:15.153] [info] HDF5 Read data rate: 1538.3150880533078 MB/s
[2025-12-15 07:18:15.153] [info] IO permuting cells
[2025-12-15 07:18:15.153] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-12-15 07:18:15.154] [info] HDF5 Read data rate: 622.4877783813146 MB/s
[2025-12-15 07:18:15.154] [info] Using partitioner with cell data (1 cell types)
[2025-12-15 07:18:15.154] [info] Compute partition of cells across ranks
[2025-12-15 07:18:15.154] [info] Building mesh dual graph
[2025-12-15 07:18:15.154] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:18:15.155] [info] Build nonlocal part of mesh dual graph
[2025-12-15 07:18:15.155] [info] Graph edges (local: 2590, non-local: 0)
[2025-12-15 07:18:15.155] [info] Compute graph partition using PT-SCOTCH
[2025-12-15 07:18:15.155] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.155] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.155] [info] Extract basic topology: 2828->2828
[2025-12-15 07:18:15.155] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:18:15.155] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:18:15.155] [info] GPS pseudo-diameter:(28) 706-0
[2025-12-15 07:18:15.155] [info] Create topology (generalised)
[2025-12-15 07:18:15.155] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.155] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.156] [info] Compute ghost indices
[2025-12-15 07:18:15.156] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.156] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.156] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.156] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.156] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.156] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.156] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:18:15.156] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.156] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.156] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:18:15.156] [info] Create Geometry (multiple)
[2025-12-15 07:18:15.156] [info] Got 1 dof layouts
[2025-12-15 07:18:15.156] [info] Checking required entities per dimension
[2025-12-15 07:18:15.156] [info] Cell type: 0 dofmap: 707x4
[2025-12-15 07:18:15.156] [info] Global index computation
[2025-12-15 07:18:15.156] [info] Got 1 index_maps
[2025-12-15 07:18:15.156] [info] Get global indices
[2025-12-15 07:18:15.156] [info] Calling compute_local_to_global
[2025-12-15 07:18:15.156] [info] xdofs.size = 2828
[2025-12-15 07:18:15.156] [info] dofmap sizes = 2828 
[2025-12-15 07:18:15.156] [info] all_dofmaps.size = 2828
[2025-12-15 07:18:15.156] [info] nodes.size = 191
[2025-12-15 07:18:15.156] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:18:15.156] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:18:15.156] [info] XDMF read meshtags (Cell tags)
[2025-12-15 07:18:15.156] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-12-15 07:18:15.156] [info] HDF5 Read data rate: 1604.8804710222032 MB/s
[2025-12-15 07:18:15.156] [info] IO permuting cells
[2025-12-15 07:18:15.156] [info] HDF5 Read data rate: 432.28370528890247 MB/s
[2025-12-15 07:18:15.156] [info] IO permuting cells
[2025-12-15 07:18:15.156] [info] XDMF distribute entity data
[2025-12-15 07:18:15.156] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.156] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.156] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.156] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.156] [info] XDMF build map
[2025-12-15 07:18:15.157] [info] XDMF create meshtags
[2025-12-15 07:18:15.157] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.157] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.157] [info] Computing mesh entities of dimension 2
[2025-12-15 07:18:15.157] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.157] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.157] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:15.157] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:18:15.157] [info] XDMF read meshtags (Facet tags)
[2025-12-15 07:18:15.157] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-12-15 07:18:15.157] [info] HDF5 Read data rate: 608.4363016616958 MB/s
[2025-12-15 07:18:15.157] [info] IO permuting cells
[2025-12-15 07:18:15.157] [info] HDF5 Read data rate: 156.55319848709095 MB/s
[2025-12-15 07:18:15.157] [info] IO permuting cells
[2025-12-15 07:18:15.157] [info] XDMF distribute entity data
[2025-12-15 07:18:15.157] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.157] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.157] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.157] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.157] [info] XDMF build map
[2025-12-15 07:18:15.158] [info] XDMF create meshtags
[2025-12-15 07:18:15.158] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.158] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.158] [info] Computing mesh entities of dimension 1
[2025-12-15 07:18:15.158] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:15.158] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:15.158] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:18:15.158] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:18:15.158] [info] XDMF read meshtags (Edge tags)
[2025-12-15 07:18:15.158] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-12-15 07:18:15.158] [info] HDF5 Read data rate: 82.64264264264263 MB/s
[2025-12-15 07:18:15.158] [info] IO permuting cells
[2025-12-15 07:18:15.158] [info] HDF5 Read data rate: 31.617647058823533 MB/s
[2025-12-15 07:18:15.158] [info] IO permuting cells
[2025-12-15 07:18:15.158] [info] XDMF distribute entity data
[2025-12-15 07:18:15.158] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.158] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.159] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.159] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.159] [info] XDMF build map
[2025-12-15 07:18:15.159] [info] XDMF create meshtags
[2025-12-15 07:18:15.159] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.159] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.159] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:18:15.159] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:18:15.159] [info] XDMF read meshtags (Vertex tags)
[2025-12-15 07:18:15.159] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-12-15 07:18:15.159] [info] HDF5 Read data rate: 2.152562895197094 MB/s
[2025-12-15 07:18:15.159] [info] IO permuting cells
[2025-12-15 07:18:15.159] [info] HDF5 Read data rate: 1.6197610852399271 MB/s
[2025-12-15 07:18:15.159] [info] IO permuting cells
[2025-12-15 07:18:15.159] [info] XDMF distribute entity data
[2025-12-15 07:18:15.159] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.159] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.159] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:15.159] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:15.159] [info] XDMF build map
[2025-12-15 07:18:15.159] [info] XDMF create meshtags
[2025-12-15 07:18:15.159] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:18:15.159] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:18:15.161] [info] Checking required entities per dimension
[2025-12-15 07:18:15.161] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.161] [info] Global index computation
[2025-12-15 07:18:15.161] [info] Got 2 index_maps
[2025-12-15 07:18:15.161] [info] Get global indices
[2025-12-15 07:18:15.162] [info] Compute face permutations
[2025-12-15 07:18:15.162] [info] Computing permutations for face type 0
[2025-12-15 07:18:15.162] [info] Compute edge permutations
[2025-12-15 07:18:15.166] [info] Checking required entities per dimension
[2025-12-15 07:18:15.166] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.166] [info] Global index computation
[2025-12-15 07:18:15.166] [info] Got 2 index_maps
[2025-12-15 07:18:15.166] [info] Get global indices
[2025-12-15 07:18:15.170] [info] Checking required entities per dimension
[2025-12-15 07:18:15.170] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.170] [info] Global index computation
[2025-12-15 07:18:15.170] [info] Got 2 index_maps
[2025-12-15 07:18:15.170] [info] Get global indices

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

[2025-12-15 07:18:15.178] [info] Requesting connectivity (2, 0) - (3, 0)

2. Constitutive Model#

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

Note that we use a compressible formulation with a high bulk modulus.

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

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

3. Boundary Conditions#

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

  • Neumann (Pressure): The target end-diastolic pressure \(P_{target} = 2000\) Pa applied to the endocardium.

  • Robin (Springs): To prevent rigid body motion and mimic the pericardium.

target_pressure = 2000.0
pressure = pulse.Variable(dolfinx.fem.Constant(geo.mesh, 0.0), "Pa") # Variable to ramp up pressure

# Epicardial 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_epi_perp = pulse.Variable(
    dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(2e5 / 10)),
    "Pa / m",
)
robin_epi_perp = pulse.RobinBC(
    value=alpha_epi_perp, marker=geometry.markers["EPI"][0], perpendicular=True,
)

# Endocardial Pressure
neumann = pulse.NeumannBC(traction=pressure, marker=geometry.markers["ENDO"][0])

bcs = pulse.BoundaryConditions(neumann=(neumann,), robin=(robin_epi, robin_epi_perp))

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 pressure from 0 to the target pressure.

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 displacement required to map from this fixed target back to the corresponding reference state for the current pressure.

prestress_problem = pulse.prestress.PrestressProblem(
    geometry=geometry,
    model=model,
    bcs=bcs,
    parameters={"u_space": "P_2"},
)

[2025-12-15 07:18:15.215] [info] Checking required entities per dimension
[2025-12-15 07:18:15.215] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:15.215] [info] Global index computation
[2025-12-15 07:18:15.215] [info] Got 2 index_maps
[2025-12-15 07:18:15.215] [info] Get global indices

[2025-12-15 07:18:16.343] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:16.343] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:16.343] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:16.343] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-12-15 07:18:26.092] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:26.092] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:26.092] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:26.092] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:26.092] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:26.109] [info] Column ghost size increased from 0 to 0

ramp_steps = 5
print(f"Starting prestress algorithm with target pressure {target_pressure} Pa")

Starting prestress algorithm with target pressure 2000.0 Pa

for ramp in np.linspace(0.0, 1.0, ramp_steps):
    current_p = target_pressure * ramp
    pressure.assign(current_p)
    print(f"Solving for pressure fraction: {ramp:.2f} (P = {current_p:.2f} Pa)")
    prestress_problem.solve()

Solving for pressure fraction: 0.00 (P = 0.00 Pa)

INFO:scifem.solvers:Newton iteration 1: r (abs) = 0.0 (tol=1e-06), r (rel) = 0.0 (tol=1e-08)

Solving for pressure fraction: 0.25 (P = 500.00 Pa)

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

INFO:scifem.solvers:Newton iteration 2: r (abs) = 0.005180289597146483 (tol=1e-06), r (rel) = 0.08292473546666401 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.00025071093723477944 (tol=1e-06), r (rel) = 0.00401331581158041 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 1.2573235282551247e-06 (tol=1e-06), r (rel) = 2.0126909706747157e-05 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.0028895859572815e-10 (tol=1e-06), r (rel) = 1.6053996993448036e-09 (tol=1e-08)

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

Solving for pressure fraction: 0.50 (P = 1000.00 Pa)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 0.008857807674380434 (tol=1e-06), r (rel) = 0.16227231169962789 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.0006696870553721266 (tol=1e-06), r (rel) = 0.01226846084103455 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 1.350697757398982e-05 (tol=1e-06), r (rel) = 0.0002474436740533167 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 9.697271897975652e-09 (tol=1e-06), r (rel) = 1.776510380345793e-07 (tol=1e-08)

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

Solving for pressure fraction: 0.75 (P = 1500.00 Pa)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 0.0063143775275905124 (tol=1e-06), r (rel) = 0.14424911694124193 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.0003100549487638711 (tol=1e-06), r (rel) = 0.007083065966047325 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 1.4813095563304532e-06 (tol=1e-06), r (rel) = 3.38398511149566e-05 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 8.386012594363524e-11 (tol=1e-06), r (rel) = 1.9157468905041357e-09 (tol=1e-08)

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

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

INFO:scifem.solvers:Newton iteration 2: r (abs) = 0.0037014977912994415 (tol=1e-06), r (rel) = 0.09422444124965854 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.00012592721490264593 (tol=1e-06), r (rel) = 0.0032055730224175094 (tol=1e-08)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 2.0624746712989918e-07 (tol=1e-06), r (rel) = 5.250186126046494e-06 (tol=1e-08)

We save the computed inverse displacement field \(\mathbf{u}\). This field maps: Target Geometry (\(\mathbf{x}\)) \(\to\) Reference Geometry (\(\mathbf{X}\)).

with dolfinx.io.VTXWriter(
    comm, "prestress_backward.bp", [prestress_problem.u], engine="BP4",
) as vtx:
    vtx.write(0.0)

We can also visualize comparison between the original target geometry and the recovered unloaded geometry.

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(prestress_problem.u_space)
    grid = pyvista.UnstructuredGrid(topology, cell_types, vtk_geometry)

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

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

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

2025-12-15 07:18:30.799 (   0.893s) [    7FBA6EEB2140]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

5. Verification (Forward Problem)#

To verify the result, we now explicitly deform the geometry to the recovered reference configuration.

\[ \mathbf{X}_{ref} = \mathbf{X}_{target} + \mathbf{u}_{inverse} \]

Note: geometry.deform(u) adds u to the nodal coordinates. Since we defined \(\mathbf{X} = \mathbf{x} + \mathbf{u}\), this correctly moves the mesh to \(\Omega_0\).

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

Deforming mesh to recovered reference configuration...
[2025-12-15 07:18:31.369] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:18:31.370] [info] Computed bounding box tree with 1413 nodes for 707 entities

We then solve the Forward mechanics problem starting from this new reference configuration. If we apply the target pressure to this reference geometry, the resulting deformed state should match our original target geometry.

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

[2025-12-15 07:18:31.384] [info] Checking required entities per dimension
[2025-12-15 07:18:31.384] [info] Cell type: 0 dofmap: 707x10
[2025-12-15 07:18:31.384] [info] Global index computation
[2025-12-15 07:18:31.384] [info] Got 2 index_maps
[2025-12-15 07:18:31.384] [info] Get global indices

[2025-12-15 07:18:32.499] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:32.499] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:32.500] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:32.500] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-12-15 07:18:40.607] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:40.607] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-12-15 07:18:40.607] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:40.607] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:40.608] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:40.624] [info] Column ghost size increased from 0 to 0

print("Solving forward problem to recover target geometry...")
# Ramp pressure again for robust forward solution
for ramp in np.linspace(0.0, 1.0, ramp_steps):
    current_p = target_pressure * ramp
    pressure.assign(current_p)
    forward_problem.solve()

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

Solving forward problem to recover target geometry...

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# Save the forward displacement. This maps **Reference Geometry** $\to$ **Recovered Target**.
with dolfinx.io.VTXWriter(comm, "prestress_forward.bp", [forward_problem.u], engine="BP4") as vtx:
    vtx.write(0.0)

print("Done. You can now compare the original geometry with the recovered geometry in Paraview.")

Done. You can now compare the original geometry with the recovered geometry in Paraview.

We can also visualize comparison between the original target geometry and the recovered unloaded geometry with the target pressure applied.

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(prestress_problem.u_space)
    grid = pyvista.UnstructuredGrid(topology, cell_types, vtk_geometry)

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

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

    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.

Contents