Pre-stressing of a left ventricle ellipsoid

In this demo we will show how to simulate an idealized left ventricular geometry with pre-stressing using the approach described in [BRR24].

from pathlib import Path
from mpi4py import MPI
import dolfinx
import logging

import math
import numpy as np
import pulse.prestress
import pulse
import cardiac_geometries
import cardiac_geometries.geometry

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

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-09-15 20:55:11 [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-09-15 20:55:11.703] [info] Extract basic topology: 2828->2828
[2025-09-15 20:55:11.703] [info] Build local dual graph
[2025-09-15 20:55:11.703] [info] Build local part of mesh dual graph (mixed)
[2025-09-15 20:55:11.704] [info] GPS pseudo-diameter:(28) 58-440
[2025-09-15 20:55:11.704] [info] Create topology (single cell type)
[2025-09-15 20:55:11.704] [info] Create topology (generalised)
[2025-09-15 20:55:11.704] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.704] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.704] [info] Compute ghost indices
[2025-09-15 20:55:11.704] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.704] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.704] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.704] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.704] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.704] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.704] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-09-15 20:55:11.704] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.704] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.704] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-09-15 20:55:11.704] [info] Checking required entities per dimension
[2025-09-15 20:55:11.704] [info] Cell type: 0 dofmap: 707x4
[2025-09-15 20:55:11.705] [info] Global index computation
[2025-09-15 20:55:11.705] [info] Got 1 index_maps
[2025-09-15 20:55:11.705] [info] Get global indices
[2025-09-15 20:55:11.705] [info] XDMF distribute entity data
[2025-09-15 20:55:11.705] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.705] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.705] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.705] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.705] [info] XDMF build map
[2025-09-15 20:55:11.705] [info] Requesting connectivity (3, 0) - (0, 0)
[2025-09-15 20:55:11.705] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.705] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.706] [info] XDMF distribute entity data
[2025-09-15 20:55:11.706] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.706] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.706] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.706] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.706] [info] XDMF build map
[2025-09-15 20:55:11.706] [info] Computing mesh entities of dimension 2
[2025-09-15 20:55:11.706] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.706] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.706] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:11.706] [info] Computing mesh connectivity 2-3 from transpose.
[2025-09-15 20:55:11.706] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.706] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.707] [info] XDMF distribute entity data
[2025-09-15 20:55:11.707] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.707] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.707] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.707] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.707] [info] XDMF build map
[2025-09-15 20:55:11.707] [info] Computing mesh entities of dimension 1
[2025-09-15 20:55:11.708] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.708] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.708] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-09-15 20:55:11.708] [info] Computing mesh connectivity 1-3 from transpose.
[2025-09-15 20:55:11.708] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.708] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.708] [info] XDMF distribute entity data
[2025-09-15 20:55:11.708] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.708] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.709] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.709] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.709] [info] XDMF build map
[2025-09-15 20:55:11.709] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-09-15 20:55:11.709] [info] Computing mesh connectivity 0-3 from transpose.
[2025-09-15 20:55:11.709] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.709] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.710] [info] Opened HDF5 file with id "72057594037927936"
[2025-09-15 20:55:11.710] [info] Adding mesh to node "/Xdmf/Domain"
[2025-09-15 20:55:11.710] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-09-15 20:55:11.711] [info] Adding geometry data to node "/Xdmf/Domain/Grid"
[2025-09-15 20:55:11.711] [info] XDMF: add meshtags (Cell tags)
[2025-09-15 20:55:11.711] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-09-15 20:55:11.711] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:11.711] [info] XDMF: add meshtags (Facet tags)
[2025-09-15 20:55:11.711] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-09-15 20:55:11.711] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-09-15 20:55:11.711] [info] XDMF: add meshtags (Edge tags)
[2025-09-15 20:55:11.711] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-09-15 20:55:11.711] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-09-15 20:55:11.711] [info] XDMF: add meshtags (Vertex tags)
[2025-09-15 20:55:11.711] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-09-15 20:55:11.714] [info] Checking required entities per dimension
[2025-09-15 20:55:11.714] [info] Cell type: 0 dofmap: 707x4
[2025-09-15 20:55:11.714] [info] Global index computation
[2025-09-15 20:55:11.714] [info] Got 1 index_maps
[2025-09-15 20:55:11.714] [info] Get global indices
[2025-09-15 20:55:11.770] [info] Column ghost size increased from 0 to 0
[2025-09-15 20:55:11.774] [info] Checking required entities per dimension
[2025-09-15 20:55:11.774] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.774] [info] Global index computation
[2025-09-15 20:55:11.774] [info] Got 2 index_maps
[2025-09-15 20:55:11.774] [info] Get global indices
[2025-09-15 20:55:11.780] [info] Checking required entities per dimension
[2025-09-15 20:55:11.780] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.780] [info] Global index computation
[2025-09-15 20:55:11.780] [info] Got 2 index_maps
[2025-09-15 20:55:11.780] [info] Get global indices
[2025-09-15 20:55:11.782] [info] Compute face permutations
[2025-09-15 20:55:11.782] [info] Computing permutations for face type 0
[2025-09-15 20:55:11.782] [info] Compute edge permutations
[2025-09-15 20:55:11.789] [info] Opened HDF5 file with id "72057594037927937"
[2025-09-15 20:55:11.789] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-09-15 20:55:11.789] [info] HDF5 Read data rate: 945.2266555253813 MB/s
[2025-09-15 20:55:11.789] [info] IO permuting cells
[2025-09-15 20:55:11.789] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-09-15 20:55:11.789] [info] HDF5 Read data rate: 576.2413576367065 MB/s
[2025-09-15 20:55:11.790] [info] Using partitioner with 2828 cell data
[2025-09-15 20:55:11.790] [info] Compute partition of cells across ranks
[2025-09-15 20:55:11.790] [info] Building mesh dual graph
[2025-09-15 20:55:11.790] [info] Build local part of mesh dual graph (mixed)
[2025-09-15 20:55:11.790] [info] Build nonlocal part of mesh dual graph
[2025-09-15 20:55:11.790] [info] Graph edges (local: 2590, non-local: 0)
[2025-09-15 20:55:11.790] [info] Compute graph partition using PT-SCOTCH
[2025-09-15 20:55:11.790] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.790] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.790] [info] Extract basic topology: 2828->2828
[2025-09-15 20:55:11.790] [info] Build local dual graph
[2025-09-15 20:55:11.790] [info] Build local part of mesh dual graph (mixed)
[2025-09-15 20:55:11.791] [info] GPS pseudo-diameter:(28) 706-0
[2025-09-15 20:55:11.791] [info] Create topology (single cell type)
[2025-09-15 20:55:11.791] [info] Create topology (generalised)
[2025-09-15 20:55:11.791] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.791] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.791] [info] Compute ghost indices
[2025-09-15 20:55:11.791] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.791] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.791] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.791] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.791] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.791] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.791] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-09-15 20:55:11.791] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.791] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.791] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-09-15 20:55:11.791] [info] Checking required entities per dimension
[2025-09-15 20:55:11.791] [info] Cell type: 0 dofmap: 707x4
[2025-09-15 20:55:11.791] [info] Global index computation
[2025-09-15 20:55:11.791] [info] Got 1 index_maps
[2025-09-15 20:55:11.791] [info] Get global indices
[2025-09-15 20:55:11.791] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-09-15 20:55:11.792] [info] XDMF read meshtags (Cell tags)
[2025-09-15 20:55:11.792] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-09-15 20:55:11.792] [info] HDF5 Read data rate: 1754.614549402823 MB/s
[2025-09-15 20:55:11.792] [info] IO permuting cells
[2025-09-15 20:55:11.792] [info] HDF5 Read data rate: 391.5270663159352 MB/s
[2025-09-15 20:55:11.792] [info] IO permuting cells
[2025-09-15 20:55:11.792] [info] XDMF distribute entity data
[2025-09-15 20:55:11.792] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.792] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.792] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.792] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.792] [info] XDMF build map
[2025-09-15 20:55:11.792] [info] XDMF create meshtags
[2025-09-15 20:55:11.792] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.792] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.792] [info] Computing mesh entities of dimension 2
[2025-09-15 20:55:11.793] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.793] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.793] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:11.793] [info] Computing mesh connectivity 2-3 from transpose.
[2025-09-15 20:55:11.793] [info] XDMF read meshtags (Facet tags)
[2025-09-15 20:55:11.793] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-09-15 20:55:11.793] [info] HDF5 Read data rate: 553.0060993319779 MB/s
[2025-09-15 20:55:11.793] [info] IO permuting cells
[2025-09-15 20:55:11.793] [info] HDF5 Read data rate: 155.27646387212525 MB/s
[2025-09-15 20:55:11.793] [info] IO permuting cells
[2025-09-15 20:55:11.793] [info] XDMF distribute entity data
[2025-09-15 20:55:11.793] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.793] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.793] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.793] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.793] [info] XDMF build map
[2025-09-15 20:55:11.793] [info] XDMF create meshtags
[2025-09-15 20:55:11.793] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.793] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.793] [info] Computing mesh entities of dimension 1
[2025-09-15 20:55:11.794] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.794] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.794] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-09-15 20:55:11.794] [info] Computing mesh connectivity 1-3 from transpose.
[2025-09-15 20:55:11.794] [info] XDMF read meshtags (Edge tags)
[2025-09-15 20:55:11.794] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-09-15 20:55:11.794] [info] HDF5 Read data rate: 81.84630026171781 MB/s
[2025-09-15 20:55:11.794] [info] IO permuting cells
[2025-09-15 20:55:11.794] [info] HDF5 Read data rate: 30.228471001757473 MB/s
[2025-09-15 20:55:11.794] [info] IO permuting cells
[2025-09-15 20:55:11.794] [info] XDMF distribute entity data
[2025-09-15 20:55:11.794] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.794] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.794] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.794] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.794] [info] XDMF build map
[2025-09-15 20:55:11.794] [info] XDMF create meshtags
[2025-09-15 20:55:11.794] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.794] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.794] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-09-15 20:55:11.794] [info] Computing mesh connectivity 0-3 from transpose.
[2025-09-15 20:55:11.794] [info] XDMF read meshtags (Vertex tags)
[2025-09-15 20:55:11.794] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-09-15 20:55:11.794] [info] HDF5 Read data rate: 2.071465561885034 MB/s
[2025-09-15 20:55:11.794] [info] IO permuting cells
[2025-09-15 20:55:11.794] [info] HDF5 Read data rate: 1.6464293064416546 MB/s
[2025-09-15 20:55:11.794] [info] IO permuting cells
[2025-09-15 20:55:11.794] [info] XDMF distribute entity data
[2025-09-15 20:55:11.794] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.794] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.794] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.794] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.794] [info] XDMF build map
[2025-09-15 20:55:11.795] [info] XDMF create meshtags
[2025-09-15 20:55:11.795] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.795] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.796] [info] Checking required entities per dimension
[2025-09-15 20:55:11.796] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.797] [info] Global index computation
[2025-09-15 20:55:11.797] [info] Got 2 index_maps
[2025-09-15 20:55:11.797] [info] Get global indices
[2025-09-15 20:55:11.797] [info] Compute face permutations
[2025-09-15 20:55:11.797] [info] Computing permutations for face type 0
[2025-09-15 20:55:11.797] [info] Compute edge permutations
[2025-09-15 20:55:11.803] [info] Checking required entities per dimension
[2025-09-15 20:55:11.803] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.803] [info] Global index computation
[2025-09-15 20:55:11.803] [info] Got 2 index_maps
[2025-09-15 20:55:11.803] [info] Get global indices
[2025-09-15 20:55:11.808] [info] Checking required entities per dimension
[2025-09-15 20:55:11.808] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.808] [info] Global index computation
[2025-09-15 20:55:11.808] [info] Got 2 index_maps
[2025-09-15 20:55:11.808] [info] Get global indices

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

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

[2025-09-15 20:55:11.817] [info] Opened HDF5 file with id "72057594037927938"
[2025-09-15 20:55:11.817] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-09-15 20:55:11.817] [info] HDF5 Read data rate: 1217.391304347826 MB/s
[2025-09-15 20:55:11.817] [info] IO permuting cells
[2025-09-15 20:55:11.817] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-09-15 20:55:11.817] [info] HDF5 Read data rate: 441.61849710982654 MB/s
[2025-09-15 20:55:11.818] [info] Using partitioner with 2828 cell data
[2025-09-15 20:55:11.818] [info] Compute partition of cells across ranks
[2025-09-15 20:55:11.818] [info] Building mesh dual graph
[2025-09-15 20:55:11.818] [info] Build local part of mesh dual graph (mixed)
[2025-09-15 20:55:11.818] [info] Build nonlocal part of mesh dual graph
[2025-09-15 20:55:11.818] [info] Graph edges (local: 2590, non-local: 0)
[2025-09-15 20:55:11.818] [info] Compute graph partition using PT-SCOTCH
[2025-09-15 20:55:11.818] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.818] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.818] [info] Extract basic topology: 2828->2828
[2025-09-15 20:55:11.818] [info] Build local dual graph
[2025-09-15 20:55:11.818] [info] Build local part of mesh dual graph (mixed)
[2025-09-15 20:55:11.819] [info] GPS pseudo-diameter:(28) 706-0
[2025-09-15 20:55:11.819] [info] Create topology (single cell type)
[2025-09-15 20:55:11.819] [info] Create topology (generalised)
[2025-09-15 20:55:11.819] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.819] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.819] [info] Compute ghost indices
[2025-09-15 20:55:11.819] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.819] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.819] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.819] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.819] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.819] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.819] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-09-15 20:55:11.819] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.819] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.819] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-09-15 20:55:11.819] [info] Checking required entities per dimension
[2025-09-15 20:55:11.819] [info] Cell type: 0 dofmap: 707x4
[2025-09-15 20:55:11.819] [info] Global index computation
[2025-09-15 20:55:11.819] [info] Got 1 index_maps
[2025-09-15 20:55:11.819] [info] Get global indices
[2025-09-15 20:55:11.820] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-09-15 20:55:11.820] [info] XDMF read meshtags (Cell tags)
[2025-09-15 20:55:11.820] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-09-15 20:55:11.820] [info] HDF5 Read data rate: 1764.1921397379913 MB/s
[2025-09-15 20:55:11.820] [info] IO permuting cells
[2025-09-15 20:55:11.820] [info] HDF5 Read data rate: 440.36125817502335 MB/s
[2025-09-15 20:55:11.820] [info] IO permuting cells
[2025-09-15 20:55:11.820] [info] XDMF distribute entity data
[2025-09-15 20:55:11.820] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.820] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.820] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.820] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.820] [info] XDMF build map
[2025-09-15 20:55:11.820] [info] XDMF create meshtags
[2025-09-15 20:55:11.820] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.820] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.820] [info] Computing mesh entities of dimension 2
[2025-09-15 20:55:11.821] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.821] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.821] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:11.821] [info] Computing mesh connectivity 2-3 from transpose.
[2025-09-15 20:55:11.821] [info] XDMF read meshtags (Facet tags)
[2025-09-15 20:55:11.821] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-09-15 20:55:11.821] [info] HDF5 Read data rate: 609.8003629764065 MB/s
[2025-09-15 20:55:11.821] [info] IO permuting cells
[2025-09-15 20:55:11.821] [info] HDF5 Read data rate: 167.57613096285866 MB/s
[2025-09-15 20:55:11.821] [info] IO permuting cells
[2025-09-15 20:55:11.821] [info] XDMF distribute entity data
[2025-09-15 20:55:11.821] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.821] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.821] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.821] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.821] [info] XDMF build map
[2025-09-15 20:55:11.821] [info] XDMF create meshtags
[2025-09-15 20:55:11.821] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.821] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.822] [info] Computing mesh entities of dimension 1
[2025-09-15 20:55:11.822] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-09-15 20:55:11.822] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-09-15 20:55:11.822] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-09-15 20:55:11.822] [info] Computing mesh connectivity 1-3 from transpose.
[2025-09-15 20:55:11.822] [info] XDMF read meshtags (Edge tags)
[2025-09-15 20:55:11.822] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-09-15 20:55:11.822] [info] HDF5 Read data rate: 83.03161959932416 MB/s
[2025-09-15 20:55:11.822] [info] IO permuting cells
[2025-09-15 20:55:11.822] [info] HDF5 Read data rate: 31.38113482941069 MB/s
[2025-09-15 20:55:11.822] [info] IO permuting cells
[2025-09-15 20:55:11.822] [info] XDMF distribute entity data
[2025-09-15 20:55:11.822] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.822] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.822] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.822] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.822] [info] XDMF build map
[2025-09-15 20:55:11.822] [info] XDMF create meshtags
[2025-09-15 20:55:11.822] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.822] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.823] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-09-15 20:55:11.823] [info] Computing mesh connectivity 0-3 from transpose.
[2025-09-15 20:55:11.823] [info] XDMF read meshtags (Vertex tags)
[2025-09-15 20:55:11.823] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-09-15 20:55:11.823] [info] HDF5 Read data rate: 1.2389654638376955 MB/s
[2025-09-15 20:55:11.823] [info] IO permuting cells
[2025-09-15 20:55:11.823] [info] HDF5 Read data rate: 0.7072761029086729 MB/s
[2025-09-15 20:55:11.823] [info] IO permuting cells
[2025-09-15 20:55:11.823] [info] XDMF distribute entity data
[2025-09-15 20:55:11.823] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.823] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.823] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-09-15 20:55:11.823] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-09-15 20:55:11.823] [info] XDMF build map
[2025-09-15 20:55:11.823] [info] XDMF create meshtags
[2025-09-15 20:55:11.823] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-09-15 20:55:11.823] [info] Build list of mesh entity indices from the entity vertices.
[2025-09-15 20:55:11.825] [info] Checking required entities per dimension
[2025-09-15 20:55:11.825] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.825] [info] Global index computation
[2025-09-15 20:55:11.825] [info] Got 2 index_maps
[2025-09-15 20:55:11.825] [info] Get global indices
[2025-09-15 20:55:11.826] [info] Compute face permutations
[2025-09-15 20:55:11.826] [info] Computing permutations for face type 0
[2025-09-15 20:55:11.826] [info] Compute edge permutations
[2025-09-15 20:55:11.830] [info] Checking required entities per dimension
[2025-09-15 20:55:11.830] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.830] [info] Global index computation
[2025-09-15 20:55:11.830] [info] Got 2 index_maps
[2025-09-15 20:55:11.830] [info] Get global indices
[2025-09-15 20:55:11.835] [info] Checking required entities per dimension
[2025-09-15 20:55:11.835] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.835] [info] Global index computation
[2025-09-15 20:55:11.835] [info] Got 2 index_maps
[2025-09-15 20:55:11.835] [info] Get global indices

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

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

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

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

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

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

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

[2025-09-15 20:55:11.869] [info] Checking required entities per dimension
[2025-09-15 20:55:11.869] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:11.869] [info] Global index computation
[2025-09-15 20:55:11.869] [info] Got 2 index_maps
[2025-09-15 20:55:11.869] [info] Get global indices

[2025-09-15 20:55:13.028] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:13.028] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-15 20:55:13.028] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:13.028] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-09-15 20:55:23.007] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:23.007] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-15 20:55:23.007] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:23.007] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-15 20:55:23.007] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:23.023] [info] Column ghost size increased from 0 to 0

ramp_steps = 5
for ramp in np.linspace(0.0, 1.0, ramp_steps):
    pressure.assign(target_pressure * ramp)
    prestress_problem.solve()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

geometry.deform(prestress_problem.u)

[2025-09-15 20:55:26.754] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-09-15 20:55:26.755] [info] Computed bounding box tree with 1413 nodes for 707 entities

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

[2025-09-15 20:55:26.762] [info] Checking required entities per dimension
[2025-09-15 20:55:26.762] [info] Cell type: 0 dofmap: 707x10
[2025-09-15 20:55:26.762] [info] Global index computation
[2025-09-15 20:55:26.762] [info] Got 2 index_maps
[2025-09-15 20:55:26.762] [info] Get global indices

[2025-09-15 20:55:27.806] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:27.806] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-15 20:55:27.806] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:27.806] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-09-15 20:55:36.131] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:36.131] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-09-15 20:55:36.131] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:36.131] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-15 20:55:36.131] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-15 20:55:36.148] [info] Column ghost size increased from 0 to 0

for ramp in np.linspace(0.0, 1.0, ramp_steps):
    pressure.assign(target_pressure * ramp)
    forward_problem.solve()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.