Bi-Ventricular Ellipsoid Simulation

Bi-Ventricular Ellipsoid Simulation#

This demo illustrates how to simulate the mechanics of an idealized Bi-Ventricular (BiV) geometry. Unlike the Left Ventricle (LV) only model, this includes both the LV and Right Ventricle (RV), allowing us to study interventricular interactions.

Problem Setup#

Geometry: An idealized BiV geometry is generated using cardiac-geometries. It consists of two truncated ellipsoids joined together.

Microstructure (Fibers): For BiV geometries, analytical fiber definitions are complex. We use the Laplace-Dirichlet Rule-Based (LDRB) algorithm [Bayer et al. 2012] implemented in fenicsx-ldrb to generate realistic fiber (\(\mathbf{f}_0\)), sheet (\(\mathbf{s}_0\)), and sheet-normal (\(\mathbf{n}_0\)) fields.

Physics: We solve the static balance of linear momentum:

\[ \nabla \cdot \mathbf{P} = \mathbf{0} \quad \text{in } \Omega_0 \]

Material Model:

  • Passive: Neo-Hookean (or Holzapfel-Ogden).

  • Active: Active stress along fibers.

  • Compressibility: Either Incompressible or Compressible formulation.

Imports#

from pathlib import Path
from mpi4py import MPI
import numpy as np
import dolfinx
from dolfinx import log
import ldrb
import cardiac_geometries
import cardiac_geometries.geometry
import pulse

We enable info logging to track solver progress.

log.set_log_level(log.LogLevel.INFO)

Geometry and Microstructure#

We generate the mesh and fibers if they don’t already exist.

1. Mesh Generation#

cardiac_geometries.mesh.biv_ellipsoid creates the mesh with markers for:

  • LV_ENDO_FW: LV Endocardium Free Wall

  • LV_SEPTUM: LV Septum

  • RV_ENDO_FW: RV Endocardium Free Wall

  • RV_SEPTUM: RV Septum

  • LV_EPI_FW: LV Epicardium Free Wall

  • RV_EPI_FW: RV Epicardium Free Wall

  • BASE: Base

2. Fiber Generation (LDRB)#

The LDRB algorithm solves Laplace equations to define transmural and apicobasal coordinates. Based on these coordinates, it assigns fiber angles (e.g., +60 to -60 degrees transmurally).

outdir = Path("biv_ellipsoid")
outdir.mkdir(parents=True, exist_ok=True)
geodir = outdir / "geometry"

if not geodir.exists():
    # Generate mesh
    geo = cardiac_geometries.mesh.biv_ellipsoid(outdir=geodir)

    # Map mesh markers to the format expected by LDRB
    markers = cardiac_geometries.mesh.transform_biv_markers(geo.markers)

    # Run LDRB algorithm
    system = ldrb.dolfinx_ldrb(
        mesh=geo.mesh,
        ffun=geo.ffun,
        markers=markers,
        alpha_endo_lv=60,  # Fiber angle at LV endocardium
        alpha_epi_lv=-60,  # Fiber angle at LV epicardium
        beta_endo_lv=0,    # Sheet angle (0 for now)
        beta_epi_lv=0,
        fiber_space="P_2",
    )

    # Save microstructure to XDMF/H5
    cardiac_geometries.fibers.utils.save_microstructure(
        mesh=geo.mesh,
        functions=[system.f0, system.s0, system.n0],
        outdir=geodir,
    )

2025-12-15 07:16:52 [debug    ] Convert file biv_ellipsoid/geometry/biv_ellipsoid.msh to dolfin

Info    : Reading 'biv_ellipsoid/geometry/biv_ellipsoid.msh'...
Info    : 64 entities
Info    : 1472 nodes
Info    : 7421 elements
Info    : Done reading 'biv_ellipsoid/geometry/biv_ellipsoid.msh'
[2025-12-15 07:16:52.504] [info] Extract basic topology: 17924->17924
[2025-12-15 07:16:52.504] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:16:52.504] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:16:52.508] [info] GPS pseudo-diameter:(72) 2399-4155
[2025-12-15 07:16:52.508] [info] Create topology (generalised)
[2025-12-15 07:16:52.509] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.509] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.510] [info] Compute ghost indices
[2025-12-15 07:16:52.510] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.510] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.511] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.511] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.511] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.511] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.511] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:16:52.511] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.511] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.511] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:16:52.511] [info] Create Geometry (multiple)
[2025-12-15 07:16:52.511] [info] Got 1 dof layouts
[2025-12-15 07:16:52.511] [info] Checking required entities per dimension
[2025-12-15 07:16:52.511] [info] Cell type: 0 dofmap: 4481x4
[2025-12-15 07:16:52.511] [info] Global index computation
[2025-12-15 07:16:52.511] [info] Got 1 index_maps
[2025-12-15 07:16:52.511] [info] Get global indices
[2025-12-15 07:16:52.511] [info] Calling compute_local_to_global
[2025-12-15 07:16:52.511] [info] xdofs.size = 17924
[2025-12-15 07:16:52.511] [info] dofmap sizes = 17924 
[2025-12-15 07:16:52.511] [info] all_dofmaps.size = 17924
[2025-12-15 07:16:52.511] [info] nodes.size = 1472
[2025-12-15 07:16:52.511] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:16:52.512] [info] XDMF distribute entity data
[2025-12-15 07:16:52.512] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.512] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.512] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.512] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.512] [info] XDMF build map
[2025-12-15 07:16:52.514] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:16:52.514] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:52.514] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:52.516] [info] XDMF distribute entity data
[2025-12-15 07:16:52.516] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.516] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.516] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.516] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.517] [info] XDMF build map
[2025-12-15 07:16:52.518] [info] Computing mesh entities of dimension 2
[2025-12-15 07:16:52.521] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.521] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.521] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:16:52.521] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:16:52.521] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:52.521] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:52.526] [info] Opened HDF5 file with id "72057594037927936"
[2025-12-15 07:16:52.526] [info] Adding mesh to node "/Xdmf/Domain"
[2025-12-15 07:16:52.526] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:16:52.526] [info] Adding geometry data to node "/Xdmf/Domain/Grid"
[2025-12-15 07:16:52.526] [info] XDMF: add meshtags (Cell tags)
[2025-12-15 07:16:52.526] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:16:52.527] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:16:52.527] [info] XDMF: add meshtags (Facet tags)
[2025-12-15 07:16:52.527] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:16:52.527] [info] Computing mesh entities of dimension 1
[2025-12-15 07:16:52.530] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.530] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.531] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:16:52.531] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:16:52.531] [info] XDMF: add meshtags (Edge tags)
[2025-12-15 07:16:52.531] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:16:52.531] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:16:52.531] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:16:52.531] [info] XDMF: add meshtags (Vertex tags)
[2025-12-15 07:16:52.531] [info] Adding topology data to node /Xdmf/Domain/Grid
[2025-12-15 07:16:52.532] [info] Opened HDF5 file with id "72057594037927937"
[2025-12-15 07:16:52.532] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-12-15 07:16:52.532] [info] HDF5 Read data rate: 4594.719302742887 MB/s
[2025-12-15 07:16:52.532] [info] IO permuting cells
[2025-12-15 07:16:52.532] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-12-15 07:16:52.532] [info] HDF5 Read data rate: 1470.5294705294707 MB/s
[2025-12-15 07:16:52.537] [info] Using partitioner with cell data (1 cell types)
[2025-12-15 07:16:52.537] [info] Compute partition of cells across ranks
[2025-12-15 07:16:52.537] [info] Building mesh dual graph
[2025-12-15 07:16:52.537] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:16:52.539] [info] Build nonlocal part of mesh dual graph
[2025-12-15 07:16:52.539] [info] Graph edges (local: 14984, non-local: 0)
[2025-12-15 07:16:52.539] [info] Compute graph partition using PT-SCOTCH
[2025-12-15 07:16:52.539] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.539] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.539] [info] Extract basic topology: 17924->17924
[2025-12-15 07:16:52.539] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:16:52.539] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:16:52.542] [info] GPS pseudo-diameter:(72) 4478-1
[2025-12-15 07:16:52.543] [info] Create topology (generalised)
[2025-12-15 07:16:52.543] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.543] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.544] [info] Compute ghost indices
[2025-12-15 07:16:52.544] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.544] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.545] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.545] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.545] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.545] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.545] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:16:52.545] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.545] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.545] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:16:52.545] [info] Create Geometry (multiple)
[2025-12-15 07:16:52.545] [info] Got 1 dof layouts
[2025-12-15 07:16:52.545] [info] Checking required entities per dimension
[2025-12-15 07:16:52.545] [info] Cell type: 0 dofmap: 4481x4
[2025-12-15 07:16:52.545] [info] Global index computation
[2025-12-15 07:16:52.545] [info] Got 1 index_maps
[2025-12-15 07:16:52.545] [info] Get global indices
[2025-12-15 07:16:52.545] [info] Calling compute_local_to_global
[2025-12-15 07:16:52.545] [info] xdofs.size = 17924
[2025-12-15 07:16:52.545] [info] dofmap sizes = 17924 
[2025-12-15 07:16:52.545] [info] all_dofmaps.size = 17924
[2025-12-15 07:16:52.545] [info] nodes.size = 1472
[2025-12-15 07:16:52.545] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:16:52.546] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:16:52.546] [info] XDMF read meshtags (Cell tags)
[2025-12-15 07:16:52.546] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-12-15 07:16:52.546] [info] HDF5 Read data rate: 5553.954605314122 MB/s
[2025-12-15 07:16:52.546] [info] IO permuting cells
[2025-12-15 07:16:52.546] [info] HDF5 Read data rate: 1370.9652745907908 MB/s
[2025-12-15 07:16:52.546] [info] IO permuting cells
[2025-12-15 07:16:52.546] [info] XDMF distribute entity data
[2025-12-15 07:16:52.546] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.546] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.546] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.546] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.547] [info] XDMF build map
[2025-12-15 07:16:52.548] [info] XDMF create meshtags
[2025-12-15 07:16:52.548] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:52.548] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:52.550] [info] Computing mesh entities of dimension 2
[2025-12-15 07:16:52.552] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.552] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.552] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:16:52.552] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:16:52.553] [info] XDMF read meshtags (Facet tags)
[2025-12-15 07:16:52.553] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-12-15 07:16:52.553] [info] HDF5 Read data rate: 2880.4702808621814 MB/s
[2025-12-15 07:16:52.553] [info] IO permuting cells
[2025-12-15 07:16:52.553] [info] HDF5 Read data rate: 936.8278499163547 MB/s
[2025-12-15 07:16:52.553] [info] IO permuting cells
[2025-12-15 07:16:52.553] [info] XDMF distribute entity data
[2025-12-15 07:16:52.553] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.553] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.553] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.553] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.553] [info] XDMF build map
[2025-12-15 07:16:52.555] [info] XDMF create meshtags
[2025-12-15 07:16:52.555] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:52.555] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:52.557] [info] Computing mesh entities of dimension 1
[2025-12-15 07:16:52.559] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.559] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.560] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:16:52.560] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:16:52.560] [info] XDMF read meshtags (Edge tags)
[2025-12-15 07:16:52.560] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-12-15 07:16:52.560] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:52.560] [info] IO permuting cells
[2025-12-15 07:16:52.560] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:52.560] [info] IO permuting cells
[2025-12-15 07:16:52.560] [info] XDMF distribute entity data
[2025-12-15 07:16:52.560] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.560] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.560] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.560] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.560] [info] XDMF build map
[2025-12-15 07:16:52.561] [info] XDMF create meshtags
[2025-12-15 07:16:52.561] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:52.561] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:52.563] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:16:52.563] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:16:52.563] [info] XDMF read meshtags (Vertex tags)
[2025-12-15 07:16:52.563] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-12-15 07:16:52.563] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:52.563] [info] IO permuting cells
[2025-12-15 07:16:52.563] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:52.563] [info] IO permuting cells
[2025-12-15 07:16:52.563] [info] XDMF distribute entity data
[2025-12-15 07:16:52.563] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:52.563] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:52.563] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:52.563] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:52.563] [info] XDMF build map
[2025-12-15 07:16:52.564] [info] XDMF create meshtags
[2025-12-15 07:16:52.564] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:52.564] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:52.566] [info] Checking required entities per dimension
[2025-12-15 07:16:52.566] [info] Cell type: 0 dofmap: 4481x4
[2025-12-15 07:16:52.566] [info] Global index computation
[2025-12-15 07:16:52.566] [info] Got 1 index_maps
[2025-12-15 07:16:52.566] [info] Get global indices
[2025-12-15 07:16:52.567] [info] Requesting connectivity (2, 0) - (3, 0)

[2025-12-15 07:16:53.395] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:16:53.406] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:16:53.407] [info] Requesting connectivity (2, 0) - (0, 0)
[2025-12-15 07:16:53.407] [info] Requesting connectivity (2, 0) - (0, 0)
[2025-12-15 07:16:53.407] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:16:53.407] [info] Requesting connectivity (3, 0) - (0, 0)
[2025-12-15 07:16:53.411] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:16:53.426] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:16:53.436] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:16:53.447] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:16:53.457] [info] Column ghost size increased from 0 to 0
[2025-12-15 07:16:53.466] [info] Checking required entities per dimension
[2025-12-15 07:16:53.466] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:53.467] [info] Global index computation
[2025-12-15 07:16:53.467] [info] Got 2 index_maps
[2025-12-15 07:16:53.467] [info] Get global indices
[2025-12-15 07:16:53.468] [info] Checking required entities per dimension
[2025-12-15 07:16:53.468] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:53.469] [info] Global index computation
[2025-12-15 07:16:53.469] [info] Got 2 index_maps
[2025-12-15 07:16:53.469] [info] Get global indices

[2025-12-15 07:16:53.929] [info] Checking required entities per dimension
[2025-12-15 07:16:53.929] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:53.929] [info] Global index computation
[2025-12-15 07:16:53.929] [info] Got 2 index_maps
[2025-12-15 07:16:53.930] [info] Get global indices
[2025-12-15 07:16:53.931] [info] Checking required entities per dimension
[2025-12-15 07:16:53.931] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:53.931] [info] Global index computation
[2025-12-15 07:16:53.931] [info] Got 2 index_maps
[2025-12-15 07:16:53.931] [info] Get global indices

[2025-12-15 07:16:55.044] [info] Checking required entities per dimension
[2025-12-15 07:16:55.044] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:55.044] [info] Global index computation
[2025-12-15 07:16:55.044] [info] Got 2 index_maps
[2025-12-15 07:16:55.044] [info] Get global indices
[2025-12-15 07:16:55.046] [info] Checking required entities per dimension
[2025-12-15 07:16:55.046] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:55.047] [info] Global index computation
[2025-12-15 07:16:55.047] [info] Got 2 index_maps
[2025-12-15 07:16:55.047] [info] Get global indices
[2025-12-15 07:16:55.052] [info] Compute face permutations
[2025-12-15 07:16:55.052] [info] Computing permutations for face type 0
[2025-12-15 07:16:55.053] [info] Compute edge permutations

# Load the geometry with fibers
geo = cardiac_geometries.geometry.Geometry.from_folder(
    comm=MPI.COMM_WORLD,
    folder=geodir,
)

[2025-12-15 07:16:55.069] [info] Opened HDF5 file with id "72057594037927938"
[2025-12-15 07:16:55.069] [info] Read topology data "Mesh" at /Xdmf/Domain
[2025-12-15 07:16:55.069] [info] HDF5 Read data rate: 3754.601869550417 MB/s
[2025-12-15 07:16:55.069] [info] IO permuting cells
[2025-12-15 07:16:55.069] [info] Read geometry data "Mesh" at /Xdmf/Domain
[2025-12-15 07:16:55.069] [info] HDF5 Read data rate: 1254.9019607843136 MB/s
[2025-12-15 07:16:55.074] [info] Using partitioner with cell data (1 cell types)
[2025-12-15 07:16:55.074] [info] Compute partition of cells across ranks
[2025-12-15 07:16:55.074] [info] Building mesh dual graph
[2025-12-15 07:16:55.074] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:16:55.076] [info] Build nonlocal part of mesh dual graph
[2025-12-15 07:16:55.076] [info] Graph edges (local: 14984, non-local: 0)
[2025-12-15 07:16:55.076] [info] Compute graph partition using PT-SCOTCH
[2025-12-15 07:16:55.076] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.076] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.076] [info] Extract basic topology: 17924->17924
[2025-12-15 07:16:55.076] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:16:55.076] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:16:55.079] [info] GPS pseudo-diameter:(72) 4478-1
[2025-12-15 07:16:55.080] [info] Create topology (generalised)
[2025-12-15 07:16:55.080] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.080] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.081] [info] Compute ghost indices
[2025-12-15 07:16:55.081] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.081] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.082] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.082] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.082] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.082] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.082] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:16:55.082] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.082] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.082] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:16:55.082] [info] Create Geometry (multiple)
[2025-12-15 07:16:55.082] [info] Got 1 dof layouts
[2025-12-15 07:16:55.082] [info] Checking required entities per dimension
[2025-12-15 07:16:55.082] [info] Cell type: 0 dofmap: 4481x4
[2025-12-15 07:16:55.082] [info] Global index computation
[2025-12-15 07:16:55.082] [info] Got 1 index_maps
[2025-12-15 07:16:55.082] [info] Get global indices
[2025-12-15 07:16:55.082] [info] Calling compute_local_to_global
[2025-12-15 07:16:55.082] [info] xdofs.size = 17924
[2025-12-15 07:16:55.082] [info] dofmap sizes = 17924 
[2025-12-15 07:16:55.082] [info] all_dofmaps.size = 17924
[2025-12-15 07:16:55.082] [info] nodes.size = 1472
[2025-12-15 07:16:55.082] [info] Creating geometry with 1 dofmaps
[2025-12-15 07:16:55.083] [info] Requesting connectivity (3, 0) - (3, 0)
[2025-12-15 07:16:55.083] [info] XDMF read meshtags (Cell tags)
[2025-12-15 07:16:55.083] [info] Read topology data "Cell tags" at /Xdmf/Domain
[2025-12-15 07:16:55.083] [info] HDF5 Read data rate: 5553.954605314122 MB/s
[2025-12-15 07:16:55.083] [info] IO permuting cells
[2025-12-15 07:16:55.083] [info] HDF5 Read data rate: 1556.984016678249 MB/s
[2025-12-15 07:16:55.083] [info] IO permuting cells
[2025-12-15 07:16:55.083] [info] XDMF distribute entity data
[2025-12-15 07:16:55.083] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.083] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.083] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.083] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.084] [info] XDMF build map
[2025-12-15 07:16:55.085] [info] XDMF create meshtags
[2025-12-15 07:16:55.085] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:55.085] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:55.087] [info] Computing mesh entities of dimension 2
[2025-12-15 07:16:55.089] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.089] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.089] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:16:55.089] [info] Computing mesh connectivity 2-3 from transpose.
[2025-12-15 07:16:55.090] [info] XDMF read meshtags (Facet tags)
[2025-12-15 07:16:55.091] [info] Read topology data "Facet tags" at /Xdmf/Domain
[2025-12-15 07:16:55.091] [info] HDF5 Read data rate: 3961.151967664066 MB/s
[2025-12-15 07:16:55.091] [info] IO permuting cells
[2025-12-15 07:16:55.091] [info] HDF5 Read data rate: 1381.0921902524956 MB/s
[2025-12-15 07:16:55.091] [info] IO permuting cells
[2025-12-15 07:16:55.091] [info] XDMF distribute entity data
[2025-12-15 07:16:55.091] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.091] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.091] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.091] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.091] [info] XDMF build map
[2025-12-15 07:16:55.092] [info] XDMF create meshtags
[2025-12-15 07:16:55.092] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:55.092] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:55.095] [info] Computing mesh entities of dimension 1
[2025-12-15 07:16:55.097] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.097] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.098] [info] Requesting connectivity (1, 0) - (3, 0)
[2025-12-15 07:16:55.098] [info] Computing mesh connectivity 1-3 from transpose.
[2025-12-15 07:16:55.098] [info] XDMF read meshtags (Edge tags)
[2025-12-15 07:16:55.098] [info] Read topology data "Edge tags" at /Xdmf/Domain
[2025-12-15 07:16:55.098] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:55.098] [info] IO permuting cells
[2025-12-15 07:16:55.098] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:55.099] [info] IO permuting cells
[2025-12-15 07:16:55.099] [info] XDMF distribute entity data
[2025-12-15 07:16:55.099] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.099] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.099] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.099] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.099] [info] XDMF build map
[2025-12-15 07:16:55.100] [info] XDMF create meshtags
[2025-12-15 07:16:55.100] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:55.100] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:55.101] [info] Requesting connectivity (0, 0) - (3, 0)
[2025-12-15 07:16:55.101] [info] Computing mesh connectivity 0-3 from transpose.
[2025-12-15 07:16:55.101] [info] XDMF read meshtags (Vertex tags)
[2025-12-15 07:16:55.101] [info] Read topology data "Vertex tags" at /Xdmf/Domain
[2025-12-15 07:16:55.101] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:55.101] [info] IO permuting cells
[2025-12-15 07:16:55.101] [info] HDF5 Read data rate: 0 MB/s
[2025-12-15 07:16:55.101] [info] IO permuting cells
[2025-12-15 07:16:55.101] [info] XDMF distribute entity data
[2025-12-15 07:16:55.101] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:16:55.101] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:16:55.101] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:16:55.101] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:16:55.101] [info] XDMF build map
[2025-12-15 07:16:55.102] [info] XDMF create meshtags
[2025-12-15 07:16:55.102] [info] Building MeshTags object from tagged entities (defined by vertices).
[2025-12-15 07:16:55.102] [info] Build list of mesh entity indices from the entity vertices.
[2025-12-15 07:16:55.106] [info] Checking required entities per dimension
[2025-12-15 07:16:55.106] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:55.106] [info] Global index computation
[2025-12-15 07:16:55.106] [info] Got 2 index_maps
[2025-12-15 07:16:55.107] [info] Get global indices
[2025-12-15 07:16:55.108] [info] Compute face permutations
[2025-12-15 07:16:55.108] [info] Computing permutations for face type 0
[2025-12-15 07:16:55.108] [info] Compute edge permutations
[2025-12-15 07:16:55.120] [info] Checking required entities per dimension
[2025-12-15 07:16:55.120] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:55.120] [info] Global index computation
[2025-12-15 07:16:55.120] [info] Got 2 index_maps
[2025-12-15 07:16:55.120] [info] Get global indices
[2025-12-15 07:16:55.132] [info] Checking required entities per dimension
[2025-12-15 07:16:55.132] [info] Cell type: 0 dofmap: 4481x10
[2025-12-15 07:16:55.133] [info] Global index computation
[2025-12-15 07:16:55.133] [info] Got 2 index_maps
[2025-12-15 07:16:55.133] [info] Get global indices

Marker Consolidation#

We group the detailed surface markers into broader categories for applying boundary conditions.

  • ENDO_LV: All LV endocardial surfaces (Free Wall + Septum).

  • ENDO_RV: All RV endocardial surfaces (Free Wall + Septum).

  • EPI: All epicardial surfaces.

  • BASE: The base.

markers = {"ENDO_LV": [1, 2], "ENDO_RV": [2, 2], "BASE": [3, 2], "EPI": [4, 2]}

# Create a new marker array based on the old one
marker_values = geo.ffun.values.copy()
marker_indices = geo.ffun.indices

# Helper to update markers
def update_marker(old_name, new_id):
    old_id = geo.markers[old_name][0]
    marker_values[np.isin(marker_indices, geo.ffun.find(old_id))] = new_id

# Update LV Endo
update_marker("LV_ENDO_FW", markers["ENDO_LV"][0])
update_marker("LV_SEPTUM", markers["ENDO_LV"][0])

# Update RV Endo
update_marker("RV_ENDO_FW", markers["ENDO_RV"][0])
update_marker("RV_SEPTUM", markers["ENDO_RV"][0])

# Update Base
update_marker("BASE", markers["BASE"][0])

# Update Epi
update_marker("LV_EPI_FW", markers["EPI"][0])
update_marker("RV_EPI_FW", markers["EPI"][0])

# Assign back to geometry object
geo.markers = markers
geo.ffun = dolfinx.mesh.meshtags(
    geo.mesh,
    geo.ffun.dim,
    geo.ffun.indices,
    marker_values,
)

# Scale the geometry from mm to meters (approximate scale factor)
geo.mesh.geometry.x[:] *= 1.4e-2

# Convert to `pulse.Geometry`
geometry = pulse.Geometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 4})

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

Constitutive Model#

1. Passive Material#

We use a Neo-Hookean model for simplicity in this demo, though Holzapfel-Ogden is also available.

\[ \Psi_{NH} = \frac{\mu}{2} (I_1 - 3) \]

deviatoric=True (default) means we use the isochoric invariant \(\bar{I}_1 = J^{-2/3} I_1\).

material = pulse.NeoHookean(mu=dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(15.0)))

Setting mu to c_6 kPa

2. Active Contraction#

We use the Active Stress model.

\[ \mathbf{S}_{active} = T_a (\mathbf{f}_0 \otimes \mathbf{f}_0) \]

Here \(T_a\) acts only along the fiber direction (\(\eta=0\)).

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

Activation is not a Variable, defaulting to kPa

3. Compressibility#

We can choose between an Incompressible formulation (using a Lagrange multiplier \(p\)) or a Compressible formulation (using a penalty function).

incompressible = False

if incompressible:
    comp_model: pulse.Compressibility = pulse.Incompressible()
else:
    # Compressible model with bulk modulus kappa (default in class)
    comp_model = pulse.Compressible()

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

Boundary Conditions#

Neumann BCs: Ventricular Pressures#

We apply separate pressures to the LV and RV endocardiums.

  • \(P_{LV}\) on \(\Gamma_{endo, LV}\)

  • \(P_{RV}\) on \(\Gamma_{endo, RV}\)

lvp = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
neumann_lv = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO_LV"][0])

Traction is not a Variable, defaulting to kPa

rvp = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
neumann_rv = pulse.NeumannBC(traction=rvp, marker=geometry.markers["ENDO_RV"][0])

Traction is not a Variable, defaulting to kPa

Robin BC: Pericardial Constraint#

We model the pericardium as a spring-like boundary condition on the epicardium.

\[ \mathbf{P}\mathbf{N} + k_{per} \mathbf{u} \cdot \mathbf{N} = 0 \quad \text{on } \Gamma_{epi} \]

This prevents rigid body motion and mimics the surrounding tissue support.

pericardium_stiffness = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0))
robin_per = pulse.RobinBC(value=pericardium_stiffness, marker=geometry.markers["EPI"][0])

Value is not a Variable, defaulting to Pa / m

We collect all BCs.

bcs = pulse.BoundaryConditions(neumann=(neumann_lv, neumann_rv), robin=(robin_per,))

Solving the Problem#

We initialize the StaticProblem. Note: base_bc=pulse.BaseBC.fixed will fix the base, preventing movement in all directions. Combined with the pericardial spring, this fully constrains the model.

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

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

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

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

[2025-12-15 07:16:59.847] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:16:59.847] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:16:59.848] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:16:59.996] [info] Column ghost size increased from 0 to 0

# Initial solve (zero load) to initialize system.
problem.solve()

1

Phase 1: Passive Inflation#

We inflate both ventricles. Typically \(P_{LV} > P_{RV}\).

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

# Pressures in kPa
lv_pressures = [0.1, 0.5, 1.0]
for i, plv in enumerate(lv_pressures, start=1):
    prv = plv * 0.2  # RV pressure is typically lower
    print(f"Solving: P_LV = {plv:.2f} kPa, P_RV = {prv:.2f} kPa")

    lvp.value = plv
    rvp.value = prv
    problem.solve()
    vtx.write(float(i))

Solving: P_LV = 0.10 kPa, P_RV = 0.02 kPa

Solving: P_LV = 0.50 kPa, P_RV = 0.10 kPa

Solving: P_LV = 1.00 kPa, P_RV = 0.20 kPa

Phase 2: Active Contraction#

We keep the pressure constant and increase active tension \(T_a\).

active_tensions = [1.0, 5.0, 10.0] # kPa
start_step = len(lv_pressures) + 1

for i, ta in enumerate(active_tensions, start=start_step):
    print(f"Solving: Ta = {ta:.2f} kPa")
    Ta.value = ta
    problem.solve()
    vtx.write(float(i))

Solving: Ta = 1.00 kPa

Solving: Ta = 5.00 kPa

Solving: Ta = 10.00 kPa

vtx.close()

Visualization#

We visualize the final state using PyVista.

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

    # Setup plotter
    p = pyvista.Plotter()
    topology, cell_types, geometry_data = dolfinx.plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_data)

    # Add reference mesh
    grid["u"] = uh.x.array.reshape((geometry_data.shape[0], 3))
    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")

    # Add deformed mesh
    warped = grid.warp_by_vector("u", factor=1.0)
    p.add_mesh(warped, show_edges=False, color="blue", opacity=1.0, label="Deformed")

    p.add_legend()
    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot(outdir / "biv_ellipsoid_pressure.png")

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

2025-12-15 07:17:20.169 (   0.996s) [    7FB379259140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=:99.0
Contents