Problem 2: Passive Inflation of a Ventricle#
This example implements Problem 2 from the cardiac mechanics benchmark suite [Land et al. 2015].
Problem Description#
Geometry: A truncated thick-walled ellipsoid representing an idealized Left Ventricle (LV).
Material: Isotropic Guccione material.
Constitutive parameters: \(C = 10.0\) kPa, \(b_f = b_t = b_{fs} = 1.0\).
The material is incompressible.
Boundary Conditions:
Base: The basal plane is fixed (\(u_x = u_y = u_z = 0\)). Note: The original benchmark specifies a sliding boundary condition in the plane, but fixing it is a common variation for stability in simple tests. We use a fixed base here.
Neumann: A pressure load \(P\) is applied to the endocardial surface. The pressure increases to 10 kPa.
Target Quantity: The inflation of the ventricle and the corresponding volume change.
from pathlib import Path
from mpi4py import MPI
from dolfinx import log
import logging
import dolfinx
import numpy as np
import math
import cardiac_geometries
import pulse
logging.basicConfig(level=logging.INFO)
1. Geometry#
We generate the truncated ellipsoid using cardiac_geometries.
\(r_{short}^{endo} = 7\) mm, \(r_{short}^{epi} = 10\) mm
\(r_{long}^{endo} = 17\) mm, \(r_{long}^{epi} = 20\) mm
Base cut at specific coordinate.
comm = MPI.COMM_WORLD
geodir = Path("lv_ellipsoid-problem2")
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),
comm=comm,
create_fibers=False, # Isotropic problem, fibers not strictly needed
)
2025-12-15 07:14:35 [debug ] Convert file lv_ellipsoid-problem2/lv_ellipsoid.msh to dolfin
INFO:cardiac_geometries.geometry:Reading geometry from lv_ellipsoid-problem2
Info : Reading 'lv_ellipsoid-problem2/lv_ellipsoid.msh'...
Info : 53 entities
Info : 771 nodes
Info : 4026 elements
Info : Done reading 'lv_ellipsoid-problem2/lv_ellipsoid.msh'
geo = cardiac_geometries.geometry.Geometry.from_folder(
comm=comm,
folder=geodir,
)
INFO:cardiac_geometries.geometry:Reading geometry from lv_ellipsoid-problem2
# We convert to `pulse.Geometry`. We assume the mesh units are mm (consistent with parameters).
geometry = pulse.Geometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 4})
# ## 2. Constitutive Model
#
# **Isotropic Guccione Model**:
# By setting $b_f = b_t = b_{fs} = 1.0$, the exponent $Q$ becomes $Q = (E_{11}^2 + E_{22}^2 + E_{33}^2 + 2E_{12}^2 + \dots) = \text{tr}(\mathbf{E}^2)$, making the model isotropic.
# For an isotropic material, the fiber direction vectors don't affect the energy (as b parameters are equal),
# but the class requires them. We can use dummy fields or the ones from the mesh.
material = pulse.Guccione(
C=pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(10.0)), "kPa"),
bf=pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0)), "dimensionless"),
bt=pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0)), "dimensionless"),
bfs=pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0)), "dimensionless"),
)
model = pulse.CardiacModel(
material=material,
active=active_model,
compressibility=comp_model,
)
3. Boundary Conditions#
Pressure: Applied to the
ENDOsurface.Base: Fixed (Dirichlet).
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. Solving#
We specify base_bc=pulse.BaseBC.fixed to clamp the base.
problem = pulse.StaticProblem(
model=model, geometry=geometry, bcs=bcs, parameters={"base_bc": pulse.BaseBC.fixed},
)
log.set_log_level(log.LogLevel.INFO)
Continuation Solver#
We ramp the pressure up to 10 kPa. We use a custom continuation loop to handle the nonlinearity.
target_pressure = 10.0
d_pressure = 1.0
current_pressure = 0.0
# Initial solve
problem.solve()
INFO:scifem.solvers:Newton iteration 1: r (abs) = 0.0 (tol=1e-06), r (rel) = 0.0 (tol=1e-10)
1
# Ramping
use_continuation = True
old_u = [problem.u.copy()]
old_p = [problem.p.copy()]
old_pressures = [pressure.value.copy()]
while pressure.value < target_pressure:
value = min(pressure.value + d_pressure, target_pressure)
print(f"Solving problem for pressure={value}")
if use_continuation and len(old_pressures) > 1:
# Better initial guess
d = (value - old_pressures[-2]) / (old_pressures[-1] - old_pressures[-2])
problem.u.x.array[:] = (1 - d) * old_u[-2].x.array + d * old_u[-1].x.array
problem.p.x.array[:] = (1 - d) * old_p[-2].x.array + d * old_p[-1].x.array
pressure.value = value
try:
nit = problem.solve()
except RuntimeError:
print("Convergence failed, reducing increment")
# Reset state and half the increment
pressure.value = old_pressures[-1]
problem.u.x.array[:] = old_u[-1].x.array
problem.p.x.array[:] = old_p[-1].x.array
d_pressure *= 0.5
problem._init_forms()
else:
print(f"Converged in {nit} iterations")
if nit < 3:
print("Increasing increment")
# Increase increment
d_pressure *= 1.5
old_u.append(problem.u.copy())
old_p.append(problem.p.copy())
old_pressures.append(pressure.value.copy())
# ## 5. Post-processing
Solving problem for pressure=1.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 23626.244412503063 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 4005.7505278976996 (tol=1e-06), r (rel) = 0.16954664727745927 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 93.06415071636705 (tol=1e-06), r (rel) = 0.003939015828817774 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.06329548859332128 (tol=1e-06), r (rel) = 2.679033006186339e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 5.776942190397818e-08 (tol=1e-06), r (rel) = 2.4451377415450115e-12 (tol=1e-10)
Converged in 5 iterations
Solving problem for pressure=2.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 30372.591565380582 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 1825.563588178098 (tol=1e-06), r (rel) = 0.06010562464675947 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 54.129347622456805 (tol=1e-06), r (rel) = 0.0017821774446193372 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.07896978832512105 (tol=1e-06), r (rel) = 2.6000345790424e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 3.904844964779365e-07 (tol=1e-06), r (rel) = 1.2856476064525894e-11 (tol=1e-10)
Converged in 5 iterations
Solving problem for pressure=3.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 40306.00232702902 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 35134.91655815246 (tol=1e-06), r (rel) = 0.8717043251543491 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 4916.066124290104 (tol=1e-06), r (rel) = 0.12196858632624583 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 227.62835748728503 (tol=1e-06), r (rel) = 0.005647505193901071 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.5194102701744425 (tol=1e-06), r (rel) = 3.76968734792022e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 0.00011788423243749854 (tol=1e-06), r (rel) = 2.924731445233057e-09 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 7: r (abs) = 9.642815094900373e-11 (tol=1e-06), r (rel) = 2.3924017610731755e-15 (tol=1e-10)
Converged in 7 iterations
Solving problem for pressure=4.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 20422.606522419763 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 8379.660769280172 (tol=1e-06), r (rel) = 0.41031299114934483 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 673.1241546056639 (tol=1e-06), r (rel) = 0.032959757309465564 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 9.112164288674096 (tol=1e-06), r (rel) = 0.0004461802796166513 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.0024589778670100677 (tol=1e-06), r (rel) = 1.2040470271562166e-07 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 2.696073307260546e-10 (tol=1e-06), r (rel) = 1.3201416304529099e-14 (tol=1e-10)
Converged in 6 iterations
Solving problem for pressure=5.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 15537.37874696135 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 2195.54104763528 (tol=1e-06), r (rel) = 0.14130704305992814 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 63.98316767813709 (tol=1e-06), r (rel) = 0.004118015575223736 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.08259278772724973 (tol=1e-06), r (rel) = 5.315747853762168e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.637900016314637e-07 (tol=1e-06), r (rel) = 1.0541675291496395e-11 (tol=1e-10)
Converged in 5 iterations
Solving problem for pressure=6.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 12303.379444282542 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 860.3455332151356 (tol=1e-06), r (rel) = 0.06992757860645707 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 9.165681268256671 (tol=1e-06), r (rel) = 0.0007449726564774055 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.0015375023974988067 (tol=1e-06), r (rel) = 1.2496586035256304e-07 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 3.3260094490482836e-10 (tol=1e-06), r (rel) = 2.703329978653874e-14 (tol=1e-10)
Converged in 5 iterations
Solving problem for pressure=7.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 10087.08218212225 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 416.0488039893757 (tol=1e-06), r (rel) = 0.04124570380984464 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 1.894181969156594 (tol=1e-06), r (rel) = 0.0001877829420794975 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 5.8117381663056444e-05 (tol=1e-06), r (rel) = 5.761565199306125e-09 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 4.486586497245761e-10 (tol=1e-06), r (rel) = 4.4478536173696716e-14 (tol=1e-10)
Converged in 5 iterations
Solving problem for pressure=8.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 8506.388361082778 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 229.37611182657773 (tol=1e-06), r (rel) = 0.026965158665455107 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.5111504097408017 (tol=1e-06), r (rel) = 6.009018023199418e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 3.7366582563431016e-06 (tol=1e-06), r (rel) = 4.3927670566259717e-10 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 5.464419107178475e-10 (tol=1e-06), r (rel) = 6.423900338454461e-14 (tol=1e-10)
Converged in 5 iterations
Solving problem for pressure=9.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 7330.422209421166 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 138.49085204257378 (tol=1e-06), r (rel) = 0.018892616016657718 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.1675785748573879 (tol=1e-06), r (rel) = 2.2860698888805265e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 3.558994465346871e-07 (tol=1e-06), r (rel) = 4.855101607616542e-11 (tol=1e-10)
Converged in 4 iterations
Solving problem for pressure=10.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 6424.687669435144 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 89.41446625275093 (tol=1e-06), r (rel) = 0.013917324989685018 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.06358874996926557 (tol=1e-06), r (rel) = 9.897562845238866e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 4.556663042549151e-08 (tol=1e-06), r (rel) = 7.092427331879576e-12 (tol=1e-10)
Converged in 4 iterations
with dolfinx.io.VTXWriter(geometry.mesh.comm, "problem2.bp", [problem.u], engine="BP4") as vtx:
vtx.write(0.0)
# Visualization
try:
import pyvista
except ImportError:
pass
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, label="Inflated")
p.show_axes()
if not pyvista.OFF_SCREEN:
p.show()
else:
p.screenshot("problem2.png")
[2025-12-15 07:15:08.082] [info] Checking required entities per dimension
[2025-12-15 07:15:08.082] [info] Cell type: 0 dofmap: 2838x4
[2025-12-15 07:15:08.082] [info] Global index computation
[2025-12-15 07:15:08.082] [info] Got 1 index_maps
[2025-12-15 07:15:08.083] [info] Get global indices
2025-12-15 07:15:08.694 ( 0.909s) [ 7F648066E140]vtkXOpenGLRenderWindow.:1458 WARN| bad X server connection. DISPLAY=:99.0
INFO:trame_server.utils.namespace:Translator(prefix=None)
INFO:wslink.backends.aiohttp:awaiting runner setup
INFO:wslink.backends.aiohttp:awaiting site startup
INFO:wslink.backends.aiohttp:Print WSLINK_READY_MSG
INFO:wslink.backends.aiohttp:Schedule auto shutdown with timout 0
INFO:wslink.backends.aiohttp:awaiting running future
# Check apex displacement
U = dolfinx.fem.Function(problem.u.function_space)
U.interpolate(lambda x: (x[0], x[1], x[2]))
# Locate apex (assumed to be the point with min z? or max z depending on orientation)
# In this mesh generation, the apex is typically at the bottom.
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: -26.523660 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: -28.172577 mm