Problem 3: Inflation and Active Contraction of a Ventricle

Problem 3: Inflation and Active Contraction of a Ventricle#

This example implements Problem 3 from the cardiac mechanics benchmark suite [Land et al. 2015].

Problem Description#

Geometry: The same truncated ellipsoid as in Problem 2.

Microstructure: A spatially varying fiber field.

  • \(\alpha_{endo} = +90^\circ\)

  • \(\alpha_{epi} = -90^\circ\)

  • The fiber angle varies linearly across the wall thickness.

Material: Transversely Isotropic Guccione material.

  • Constitutive parameters: \(C = 2.0\) kPa, \(b_f = 8.0\), \(b_t = 2.0\), \(b_{fs} = 4.0\).

  • Incompressible.

Active Contraction: An active stress \(T_a\) is applied along the fiber direction.

Loading Protocol:

  1. Inflation: Increase cavity pressure \(P\) to 15 kPa.

  2. Contraction: While maintaining pressure, increase active tension \(T_a\) to 60 kPa. (Note: The benchmark describes these as varying curves, here we usually perform a coupled ramp or sequential steps).

from pathlib import Path
from mpi4py import MPI
from dolfinx import log
import dolfinx
import numpy as np
import math
import cardiac_geometries
import pulse

log.set_log_level(log.LogLevel.INFO)

1. Geometry and Fibers#

We regenerate the mesh with the specific fiber angles required for Problem 3 (\(+90^\circ\) to \(-90^\circ\)).

geodir = Path("lv_ellipsoid-problem3")
comm = MPI.COMM_WORLD

if not geodir.exists():
    comm.barrier()
    cardiac_geometries.mesh.lv_ellipsoid(
        outdir=geodir,
        r_short_endo=7.0,
        r_short_epi=10.0,
        r_long_endo=17.0,
        r_long_epi=20.0,
        mu_apex_endo=-math.pi,
        mu_base_endo=-math.acos(5 / 17),
        mu_apex_epi=-math.pi,
        mu_base_epi=-math.acos(5 / 20),
        create_fibers=True,
        fiber_angle_endo=90,
        fiber_angle_epi=-90,
        fiber_space="Quadrature_6",  # High order quadrature space for accurate fiber representation
        comm=comm,
    )

2025-12-15 07:15:11 [debug    ] Convert file lv_ellipsoid-problem3/lv_ellipsoid.msh to dolfin

Info    : Reading 'lv_ellipsoid-problem3/lv_ellipsoid.msh'...
Info    : 53 entities
Info    : 771 nodes
Info    : 4026 elements
Info    : Done reading 'lv_ellipsoid-problem3/lv_ellipsoid.msh'
[2025-12-15 07:15:11.644] [info] Extract basic topology: 11352->11352
[2025-12-15 07:15:11.644] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:15:11.644] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:15:11.646] [info] GPS pseudo-diameter:(52) 22-2041
[2025-12-15 07:15:11.646] [info] Create topology (generalised)
[2025-12-15 07:15:11.646] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.646] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.647] [info] Compute ghost indices
[2025-12-15 07:15:11.647] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:15:11.647] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:15:11.648] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:11.648] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:11.648] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:11.648] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:11.648] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:15:11.648] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:11.648] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:11.648] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:15:11.648] [info] Create Geometry (multiple)
[2025-12-15 07:15:11.648] [info] Got 1 dof layouts
[2025-12-15 07:15:11.648] [info] Checking required entities per dimension
[2025-12-15 07:15:11.648] [info] Cell type: 0 dofmap: 2838x4
[2025-12-15 07:15:11.648] [info] Global index computation
[2025-12-15 07:15:11.648] [info] Got 1 index_maps
[2025-12-15 07:15:11.648] [info] Get global indices
[2025-12-15 07:15:11.648] [info] Calling compute_local_to_global
[2025-12-15 07:15:11.648] [info] xdofs.size = 11352
[2025-12-15 07:15:11.648] [info] dofmap sizes = 11352 
[2025-12-15 07:15:11.648] [info] all_dofmaps.size = 11352
[2025-12-15 07:15:11.648] [info] nodes.size = 771
[2025-12-15 07:15:11.648] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:15:11.648] [info] XDMF distribute entity data
[2025-12-15 07:15:11.648] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.648] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.648] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.648] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.649] [info] XDMF build map
[2025-12-15 07:15:11.650] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:15:11.650] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:11.650] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:11.651] [info] XDMF distribute entity data
[2025-12-15 07:15:11.651] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.651] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.651] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.651] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.651] [info] XDMF build map
[2025-12-15 07:15:11.652] [info] Computing mesh entities of dimension 2
[2025-12-15 07:15:11.653] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:11.653] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:11.654] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:15:11.654] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:15:11.654] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:11.654] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:11.655] [info] XDMF distribute entity data
[2025-12-15 07:15:11.655] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.655] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.655] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.655] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.655] [info] XDMF build map
[2025-12-15 07:15:11.656] [info] Computing mesh entities of dimension 1
[2025-12-15 07:15:11.659] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:11.659] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:11.659] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:15:11.659] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:15:11.659] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:11.659] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:11.660] [info] XDMF distribute entity data
[2025-12-15 07:15:11.660] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.660] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.660] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:11.660] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:11.660] [info] XDMF build map
[2025-12-15 07:15:11.661] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:15:11.661] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:15:11.661] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:11.661] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:11.663] [info] Opened HDF5 file with id "72057594037927936"
[2025-12-15 07:15:11.663] [info] Adding mesh to node "/Xdmf/Domain"
[2025-12-15 07:15:11.663] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:15:11.663] [info] Adding geometry data to node "/Xdmf/Domain/Grid"
[2025-12-15 07:15:11.663] [info] XDMF: add meshtags (Cell tags)
[2025-12-15 07:15:11.663] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:15:11.664] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:15:11.664] [info] XDMF: add meshtags (Facet tags)
[2025-12-15 07:15:11.664] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:15:11.664] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:15:11.664] [info] XDMF: add meshtags (Edge tags)
[2025-12-15 07:15:11.664] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:15:11.665] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:15:11.665] [info] XDMF: add meshtags (Vertex tags)
[2025-12-15 07:15:11.665] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:15:11.668] [info] Checking required entities per dimension
[2025-12-15 07:15:11.668] [info] Cell type: 0 dofmap: 2838x4
[2025-12-15 07:15:11.668] [info] Global index computation
[2025-12-15 07:15:11.668] [info] Got 1 index_maps
[2025-12-15 07:15:11.668] [info] Get global indices

[2025-12-15 07:15:12.374] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:15:12.381] [info] Checking required entities per dimension
[2025-12-15 07:15:12.381] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.381] [info] Global index computation
[2025-12-15 07:15:12.382] [info] Got 1 index_maps
[2025-12-15 07:15:12.382] [info] Get global indices

[2025-12-15 07:15:12.622] [info] Checking required entities per dimension
[2025-12-15 07:15:12.622] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.622] [info] Global index computation
[2025-12-15 07:15:12.622] [info] Got 1 index_maps
[2025-12-15 07:15:12.622] [info] Get global indices
[2025-12-15 07:15:12.635] [info] Compute face permutations
[2025-12-15 07:15:12.635] [info] Computing permutations for face type 0
[2025-12-15 07:15:12.636] [info] Compute edge permutations
[2025-12-15 07:15:12.654] [info] Opened HDF5 file with id "72057594037927938"
[2025-12-15 07:15:12.654] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-12-15 07:15:12.655] [info] HDF5 Read data rate: 1422.8238390674937 MB/s
[2025-12-15 07:15:12.655] [info] IO permuting cells
[2025-12-15 07:15:12.655] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-12-15 07:15:12.655] [info] HDF5 Read data rate: 1407.683529859262 MB/s
[2025-12-15 07:15:12.658] [info] Using partitioner with cell data (1 cell types)
[2025-12-15 07:15:12.658] [info] Compute partition of cells across ranks
[2025-12-15 07:15:12.658] [info] Building mesh dual graph
[2025-12-15 07:15:12.658] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:15:12.659] [info] Build nonlocal part of mesh dual graph
[2025-12-15 07:15:12.659] [info] Graph edges (local: 10252, non-local: 0)
[2025-12-15 07:15:12.659] [info] Compute graph partition using PT-SCOTCH
[2025-12-15 07:15:12.660] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.660] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.660] [info] Extract basic topology: 11352->11352
[2025-12-15 07:15:12.660] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:15:12.660] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:15:12.661] [info] GPS pseudo-diameter:(52) 2837-0
[2025-12-15 07:15:12.662] [info] Create topology (generalised)
[2025-12-15 07:15:12.662] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.662] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.662] [info] Compute ghost indices
[2025-12-15 07:15:12.662] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.662] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.663] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.663] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.663] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.663] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.663] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:15:12.663] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.663] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.663] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:15:12.663] [info] Create Geometry (multiple)
[2025-12-15 07:15:12.663] [info] Got 1 dof layouts
[2025-12-15 07:15:12.663] [info] Checking required entities per dimension
[2025-12-15 07:15:12.663] [info] Cell type: 0 dofmap: 2838x4
[2025-12-15 07:15:12.663] [info] Global index computation
[2025-12-15 07:15:12.663] [info] Got 1 index_maps
[2025-12-15 07:15:12.663] [info] Get global indices
[2025-12-15 07:15:12.663] [info] Calling compute_local_to_global
[2025-12-15 07:15:12.663] [info] xdofs.size = 11352
[2025-12-15 07:15:12.663] [info] dofmap sizes = 11352 
[2025-12-15 07:15:12.663] [info] all_dofmaps.size = 11352
[2025-12-15 07:15:12.663] [info] nodes.size = 771
[2025-12-15 07:15:12.663] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:15:12.664] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:15:12.664] [info] XDMF read meshtags (Cell tags)
[2025-12-15 07:15:12.664] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-12-15 07:15:12.664] [info] HDF5 Read data rate: 3535.4848756180168 MB/s
[2025-12-15 07:15:12.664] [info] IO permuting cells
[2025-12-15 07:15:12.664] [info] HDF5 Read data rate: 1050.0416242715753 MB/s
[2025-12-15 07:15:12.664] [info] IO permuting cells
[2025-12-15 07:15:12.664] [info] XDMF distribute entity data
[2025-12-15 07:15:12.664] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.664] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.664] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.664] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.665] [info] XDMF build map
[2025-12-15 07:15:12.665] [info] XDMF create meshtags
[2025-12-15 07:15:12.665] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.665] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.666] [info] Computing mesh entities of dimension 2
[2025-12-15 07:15:12.668] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.668] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.668] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:15:12.668] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:15:12.668] [info] XDMF read meshtags (Facet tags)
[2025-12-15 07:15:12.668] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-12-15 07:15:12.668] [info] HDF5 Read data rate: 1749.7348886532343 MB/s
[2025-12-15 07:15:12.668] [info] IO permuting cells
[2025-12-15 07:15:12.668] [info] HDF5 Read data rate: 666.4647076643442 MB/s
[2025-12-15 07:15:12.668] [info] IO permuting cells
[2025-12-15 07:15:12.668] [info] XDMF distribute entity data
[2025-12-15 07:15:12.668] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.668] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.668] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.668] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.668] [info] XDMF build map
[2025-12-15 07:15:12.669] [info] XDMF create meshtags
[2025-12-15 07:15:12.669] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.669] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.670] [info] Computing mesh entities of dimension 1
[2025-12-15 07:15:12.672] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.672] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.673] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:15:12.673] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:15:12.673] [info] XDMF read meshtags (Edge tags)
[2025-12-15 07:15:12.673] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-12-15 07:15:12.673] [info] HDF5 Read data rate: 152.78703086831 MB/s
[2025-12-15 07:15:12.673] [info] IO permuting cells
[2025-12-15 07:15:12.673] [info] HDF5 Read data rate: 62.77372262773722 MB/s
[2025-12-15 07:15:12.673] [info] IO permuting cells
[2025-12-15 07:15:12.673] [info] XDMF distribute entity data
[2025-12-15 07:15:12.673] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.673] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.673] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.673] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.673] [info] XDMF build map
[2025-12-15 07:15:12.674] [info] XDMF create meshtags
[2025-12-15 07:15:12.674] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.674] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.674] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:15:12.674] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:15:12.674] [info] XDMF read meshtags (Vertex tags)
[2025-12-15 07:15:12.674] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-12-15 07:15:12.674] [info] HDF5 Read data rate: 0.7167816503897501 MB/s
[2025-12-15 07:15:12.674] [info] IO permuting cells
[2025-12-15 07:15:12.674] [info] HDF5 Read data rate: 1.5240998285387695 MB/s
[2025-12-15 07:15:12.674] [info] IO permuting cells
[2025-12-15 07:15:12.674] [info] XDMF distribute entity data
[2025-12-15 07:15:12.674] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.674] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.674] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.674] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.675] [info] XDMF build map
[2025-12-15 07:15:12.675] [info] XDMF create meshtags
[2025-12-15 07:15:12.675] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.675] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.677] [info] Checking required entities per dimension
[2025-12-15 07:15:12.677] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.677] [info] Global index computation
[2025-12-15 07:15:12.677] [info] Got 1 index_maps
[2025-12-15 07:15:12.678] [info] Get global indices
[2025-12-15 07:15:12.679] [info] Compute face permutations
[2025-12-15 07:15:12.679] [info] Computing permutations for face type 0
[2025-12-15 07:15:12.680] [info] Compute edge permutations
[2025-12-15 07:15:12.709] [info] Checking required entities per dimension
[2025-12-15 07:15:12.709] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.709] [info] Global index computation
[2025-12-15 07:15:12.710] [info] Got 1 index_maps
[2025-12-15 07:15:12.710] [info] Get global indices
[2025-12-15 07:15:12.739] [info] Checking required entities per dimension
[2025-12-15 07:15:12.739] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.740] [info] Global index computation
[2025-12-15 07:15:12.740] [info] Got 1 index_maps
[2025-12-15 07:15:12.740] [info] Get global indices

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

[2025-12-15 07:15:12.775] [info] Opened HDF5 file with id "72057594037927939"
[2025-12-15 07:15:12.775] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-12-15 07:15:12.775] [info] HDF5 Read data rate: 2387.946675080855 MB/s
[2025-12-15 07:15:12.775] [info] IO permuting cells
[2025-12-15 07:15:12.775] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-12-15 07:15:12.775] [info] HDF5 Read data rate: 1264.193482270957 MB/s
[2025-12-15 07:15:12.778] [info] Using partitioner with cell data (1 cell types)
[2025-12-15 07:15:12.778] [info] Compute partition of cells across ranks
[2025-12-15 07:15:12.778] [info] Building mesh dual graph
[2025-12-15 07:15:12.778] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:15:12.779] [info] Build nonlocal part of mesh dual graph
[2025-12-15 07:15:12.779] [info] Graph edges (local: 10252, non-local: 0)
[2025-12-15 07:15:12.779] [info] Compute graph partition using PT-SCOTCH
[2025-12-15 07:15:12.780] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.780] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.780] [info] Extract basic topology: 11352->11352
[2025-12-15 07:15:12.780] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:15:12.780] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:15:12.781] [info] GPS pseudo-diameter:(52) 2837-0
[2025-12-15 07:15:12.782] [info] Create topology (generalised)
[2025-12-15 07:15:12.782] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.782] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.782] [info] Compute ghost indices
[2025-12-15 07:15:12.782] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.782] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.783] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.783] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.783] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.783] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.783] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:15:12.783] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.783] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.783] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:15:12.783] [info] Create Geometry (multiple)
[2025-12-15 07:15:12.783] [info] Got 1 dof layouts
[2025-12-15 07:15:12.783] [info] Checking required entities per dimension
[2025-12-15 07:15:12.783] [info] Cell type: 0 dofmap: 2838x4
[2025-12-15 07:15:12.783] [info] Global index computation
[2025-12-15 07:15:12.783] [info] Got 1 index_maps
[2025-12-15 07:15:12.783] [info] Get global indices
[2025-12-15 07:15:12.783] [info] Calling compute_local_to_global
[2025-12-15 07:15:12.783] [info] xdofs.size = 11352
[2025-12-15 07:15:12.783] [info] dofmap sizes = 11352 
[2025-12-15 07:15:12.783] [info] all_dofmaps.size = 11352
[2025-12-15 07:15:12.783] [info] nodes.size = 771
[2025-12-15 07:15:12.783] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:15:12.784] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:15:12.784] [info] XDMF read meshtags (Cell tags)
[2025-12-15 07:15:12.784] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-12-15 07:15:12.784] [info] HDF5 Read data rate: 3917.3532329724367 MB/s
[2025-12-15 07:15:12.784] [info] IO permuting cells
[2025-12-15 07:15:12.784] [info] HDF5 Read data rate: 1031.0626702997276 MB/s
[2025-12-15 07:15:12.784] [info] IO permuting cells
[2025-12-15 07:15:12.784] [info] XDMF distribute entity data
[2025-12-15 07:15:12.784] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.784] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.784] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.784] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.784] [info] XDMF build map
[2025-12-15 07:15:12.785] [info] XDMF create meshtags
[2025-12-15 07:15:12.785] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.785] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.786] [info] Computing mesh entities of dimension 2
[2025-12-15 07:15:12.788] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.788] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.788] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:15:12.788] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:15:12.788] [info] XDMF read meshtags (Facet tags)
[2025-12-15 07:15:12.788] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-12-15 07:15:12.788] [info] HDF5 Read data rate: 1690.2490556373646 MB/s
[2025-12-15 07:15:12.788] [info] IO permuting cells
[2025-12-15 07:15:12.788] [info] HDF5 Read data rate: 631.8207926479035 MB/s
[2025-12-15 07:15:12.788] [info] IO permuting cells
[2025-12-15 07:15:12.788] [info] XDMF distribute entity data
[2025-12-15 07:15:12.788] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.788] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.788] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.788] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.788] [info] XDMF build map
[2025-12-15 07:15:12.789] [info] XDMF create meshtags
[2025-12-15 07:15:12.789] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.789] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.790] [info] Computing mesh entities of dimension 1
[2025-12-15 07:15:12.792] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:15:12.792] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:15:12.792] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:15:12.792] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:15:12.792] [info] XDMF read meshtags (Edge tags)
[2025-12-15 07:15:12.792] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-12-15 07:15:12.792] [info] HDF5 Read data rate: 155.19963907060682 MB/s
[2025-12-15 07:15:12.792] [info] IO permuting cells
[2025-12-15 07:15:12.792] [info] HDF5 Read data rate: 62.888482632541134 MB/s
[2025-12-15 07:15:12.792] [info] IO permuting cells
[2025-12-15 07:15:12.792] [info] XDMF distribute entity data
[2025-12-15 07:15:12.792] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.792] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.792] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.792] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.792] [info] XDMF build map
[2025-12-15 07:15:12.793] [info] XDMF create meshtags
[2025-12-15 07:15:12.793] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.793] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.793] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:15:12.793] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:15:12.794] [info] XDMF read meshtags (Vertex tags)
[2025-12-15 07:15:12.794] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-12-15 07:15:12.794] [info] HDF5 Read data rate: 2.04211869814933 MB/s
[2025-12-15 07:15:12.794] [info] IO permuting cells
[2025-12-15 07:15:12.794] [info] HDF5 Read data rate: 1.5717092337917482 MB/s
[2025-12-15 07:15:12.794] [info] IO permuting cells
[2025-12-15 07:15:12.794] [info] XDMF distribute entity data
[2025-12-15 07:15:12.794] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.794] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.794] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:15:12.794] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:15:12.794] [info] XDMF build map
[2025-12-15 07:15:12.794] [info] XDMF create meshtags
[2025-12-15 07:15:12.794] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:15:12.794] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:15:12.796] [info] Checking required entities per dimension
[2025-12-15 07:15:12.796] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.796] [info] Global index computation
[2025-12-15 07:15:12.797] [info] Got 1 index_maps
[2025-12-15 07:15:12.797] [info] Get global indices
[2025-12-15 07:15:12.799] [info] Compute face permutations
[2025-12-15 07:15:12.799] [info] Computing permutations for face type 0
[2025-12-15 07:15:12.799] [info] Compute edge permutations

[2025-12-15 07:15:12.828] [info] Checking required entities per dimension
[2025-12-15 07:15:12.828] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.829] [info] Global index computation
[2025-12-15 07:15:12.829] [info] Got 1 index_maps
[2025-12-15 07:15:12.829] [info] Get global indices
[2025-12-15 07:15:12.859] [info] Checking required entities per dimension
[2025-12-15 07:15:12.859] [info] Cell type: 0 dofmap: 2838x24
[2025-12-15 07:15:12.859] [info] Global index computation
[2025-12-15 07:15:12.859] [info] Got 1 index_maps
[2025-12-15 07:15:12.860] [info] Get global indices

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

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

2. Constitutive Model#

Material: Transversely Isotropic Guccione. Active: Active Stress model \(S_{active} = T_a (\mathbf{f}_0 \otimes \mathbf{f}_0)\).

material_params = {
    "C": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(2.0)),
    "bf": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(8.0)),
    "bt": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(2.0)),
    "bfs": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(4.0)),
}

# Note: f0, s0, n0 come from the geometry object where they were loaded from files
material = pulse.Guccione(f0=geo.f0, s0=geo.s0, n0=geo.n0, **material_params)

Setting C to c_3 kPa

Ta = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
active_model = pulse.ActiveStress(geo.f0, activation=pulse.Variable(Ta, "kPa"))

comp_model = pulse.Incompressible()

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

3. Boundary Conditions#

  • Base: Fixed.

  • Pressure: Endocardial load.

pressure = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
neumann = pulse.NeumannBC(
    traction=pulse.Variable(pressure, "kPa"),
    marker=geo.markers["ENDO"][0],
)

bcs = pulse.BoundaryConditions(neumann=(neumann,))

4. Solution Process#

We simultaneously ramp up Pressure (\(0 \to 15\) kPa) and Active Tension (\(0 \to 60\) kPa). This represents a path in the P-Ta space.

problem = pulse.StaticProblem(
    model=model, geometry=geometry, bcs=bcs, parameters={"base_bc": pulse.BaseBC.fixed},
)

[2025-12-15 07:15:12.929] [info] Checking required entities per dimension
[2025-12-15 07:15:12.929] [info] Cell type: 0 dofmap: 2838x10
[2025-12-15 07:15:12.930] [info] Global index computation
[2025-12-15 07:15:12.930] [info] Got 2 index_maps
[2025-12-15 07:15:12.930] [info] Get global indices
[2025-12-15 07:15:12.931] [info] Checking required entities per dimension
[2025-12-15 07:15:12.931] [info] Cell type: 0 dofmap: 2838x4
[2025-12-15 07:15:12.931] [info] Global index computation
[2025-12-15 07:15:12.931] [info] Got 1 index_maps
[2025-12-15 07:15:12.931] [info] Get global indices

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

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

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

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

Initial solve

problem.solve()

1

Ramping

target_pressure = 15.0
target_Ta = 60.0
steps = 20

pressures = np.linspace(0, target_pressure, steps)
tensions = np.linspace(0, target_Ta, steps)

vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, "problem3.bp", [problem.u], engine="BP4")
vtx.write(0.0)

for i, (p, ta) in enumerate(zip(pressures, tensions)):
    if i == 0:
        continue
    print(f"Step {i}/{steps}: Pressure={p:.2f} kPa, Ta={ta:.2f} kPa")

    pressure.value = p
    Ta.value = ta

    problem.solve()
    vtx.write(float(i))
vtx.close()

Step 1/20: Pressure=0.79 kPa, Ta=3.16 kPa

Step 2/20: Pressure=1.58 kPa, Ta=6.32 kPa

Step 3/20: Pressure=2.37 kPa, Ta=9.47 kPa

Step 4/20: Pressure=3.16 kPa, Ta=12.63 kPa

Step 5/20: Pressure=3.95 kPa, Ta=15.79 kPa

Step 6/20: Pressure=4.74 kPa, Ta=18.95 kPa

Step 7/20: Pressure=5.53 kPa, Ta=22.11 kPa

Step 8/20: Pressure=6.32 kPa, Ta=25.26 kPa

Step 9/20: Pressure=7.11 kPa, Ta=28.42 kPa

Step 10/20: Pressure=7.89 kPa, Ta=31.58 kPa

Step 11/20: Pressure=8.68 kPa, Ta=34.74 kPa

Step 12/20: Pressure=9.47 kPa, Ta=37.89 kPa

Step 13/20: Pressure=10.26 kPa, Ta=41.05 kPa

Step 14/20: Pressure=11.05 kPa, Ta=44.21 kPa

Step 15/20: Pressure=11.84 kPa, Ta=47.37 kPa

Step 16/20: Pressure=12.63 kPa, Ta=50.53 kPa

Step 17/20: Pressure=13.42 kPa, Ta=53.68 kPa

Step 18/20: Pressure=14.21 kPa, Ta=56.84 kPa

Step 19/20: Pressure=15.00 kPa, Ta=60.00 kPa

5. Post-processing#

Compute apex position

U = dolfinx.fem.Function(problem.u.function_space)
U.interpolate(lambda x: (x[0], x[1], x[2]))

Use utility to evaluate at the specific vertex tag for the apex (ENDOPT/EPIPT generated by cardiac-geometries)

endo_apex_coord = pulse.utils.evaluate_at_vertex_tag(U, geo.vfun, geo.markers["ENDOPT"][0])
u_endo_apex = pulse.utils.evaluate_at_vertex_tag(problem.u, geo.vfun, geo.markers["ENDOPT"][0])
endo_apex_pos = pulse.utils.gather_broadcast_array(geo.mesh.comm, endo_apex_coord + u_endo_apex)
print(f"\nGet longitudinal position of endocardial apex: {endo_apex_pos[0, 0]:4f} mm")

Get longitudinal position of endocardial apex: -11.992595 mm

epi_apex_coord = pulse.utils.evaluate_at_vertex_tag(U, geo.vfun, geo.markers["EPIPT"][0])
u_epi_apex = pulse.utils.evaluate_at_vertex_tag(problem.u, geo.vfun, geo.markers["EPIPT"][0])
epi_apex_pos = pulse.utils.gather_broadcast_array(geo.mesh.comm, epi_apex_coord + u_epi_apex)
print(f"\nGet longitudinal position of epicardial apex: {epi_apex_pos[0, 0]:4f} mm")

Get longitudinal position of epicardial apex: -15.228496 mm

Visualization

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    V = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 1, (geometry.mesh.geometry.dim,)))
    uh = dolfinx.fem.Function(V)
    uh.interpolate(problem.u)

    p = pyvista.Plotter()
    topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))

    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")
    warped = grid.warp_by_vector("u", factor=1.0)
    p.add_mesh(warped, show_edges=True, color="blue", label="Contracted")

    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot("problem3.png")

[2025-12-15 07:16:41.561] [info] Checking required entities per dimension
[2025-12-15 07:16:41.561] [info] Cell type: 0 dofmap: 2838x4
[2025-12-15 07:16:41.562] [info] Global index computation
[2025-12-15 07:16:41.562] [info] Got 1 index_maps
[2025-12-15 07:16:41.562] [info] Get global indices

2025-12-15 07:16:42.167 (   0.843s) [    7F9247E23140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=:99.0
Contents