Problem 2: Inflation of a ventricle

In the second problem we will solve the inflation of a ventricle. First we import the necessary libraries

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)
# Next we will create the geometry and save it in the folder called `lv_ellipsoid`.

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

2025-07-01 08:29:48 [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'

If the folder already exist, then we just load the geometry

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

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

Now, lets convert the geometry to a pulse.Geometry object.

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

The material model used in this benchmark is the Guccione model.

material_params = {
    "C": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(10.0)),
    "bf": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0)),
    "bt": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0)),
    "bfs": dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0)),
}
material = pulse.Guccione(**material_params)

WARNING:pulse.material_models.guccione:Setting C to c_0 kPa

There are now active contraction, so we choose a pure passive model

active_model = pulse.active_model.Passive()

and the model should be incompressible

comp_model = pulse.Incompressible()

and assembles the CardiacModel

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

We apply a traction in endocardium

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

and finally combine all the boundary conditions

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

and create a Mixed problem

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

# Now we can solve the problem
log.set_log_level(log.LogLevel.INFO)

problem.solve()

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

1

# Now step up the pressure to 10 kPa starting with an increment of 1 kPa
target_value = 10.0
incr = 1.0

Here we use a continuation strategy to speed up the convergence

use_continuation = True

old_u = [problem.u.copy()]
old_p = [problem.p.copy()]
old_tractions = [traction.value.copy()]

while traction.value < target_value:
    value = min(traction.value + incr, target_value)
    print(f"Solving problem for traction={value}")

    if use_continuation and len(old_tractions) > 1:
        # Better initial guess
        d = (value - old_tractions[-2]) / (old_tractions[-1] - old_tractions[-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

    traction.value = value

    try:
        nit = problem.solve()
    except RuntimeError:
        print("Convergence failed, reducing increment")

        # Reset state and half the increment
        traction.value = old_tractions[-1]
        problem.u.x.array[:] = old_u[-1].x.array
        problem.p.x.array[:] = old_p[-1].x.array
        incr *= 0.5
        problem._init_forms()
    else:
        print(f"Converged in {nit} iterations")
        if nit < 3:
            print("Increasing increment")
            # Increase increment
            incr *= 1.5
        old_u.append(problem.u.copy())
        old_p.append(problem.p.copy())
        old_tractions.append(traction.value.copy())

Solving problem for traction=1.0

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

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

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

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

INFO:scifem.solvers:Newton iteration 5: r (abs) = 5.776990451937818e-08 (tol=1e-06), r (rel) = 2.4451581687171067e-12 (tol=1e-10)

Converged in 5 iterations
Solving problem for traction=2.0

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

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

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

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

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

Converged in 5 iterations
Solving problem for traction=3.0

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

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

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

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

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

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

INFO:scifem.solvers:Newton iteration 7: r (abs) = 8.727628928560224e-11 (tol=1e-06), r (rel) = 2.165342238120632e-15 (tol=1e-10)

Converged in 7 iterations
Solving problem for traction=4.0

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

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

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

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

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

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

Converged in 6 iterations
Solving problem for traction=5.0

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

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

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

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

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

Converged in 5 iterations
Solving problem for traction=6.0

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

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

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

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

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

Converged in 5 iterations
Solving problem for traction=7.0

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

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

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

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

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

Converged in 5 iterations
Solving problem for traction=8.0

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

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

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

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

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

Converged in 5 iterations
Solving problem for traction=9.0

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

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

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

INFO:scifem.solvers:Newton iteration 4: r (abs) = 3.558908852702621e-07 (tol=1e-06), r (rel) = 4.8549848164558386e-11 (tol=1e-10)

Converged in 4 iterations
Solving problem for traction=10.0

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

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

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

INFO:scifem.solvers:Newton iteration 4: r (abs) = 4.57375195048159e-08 (tol=1e-06), r (rel) = 7.119026146590678e-12 (tol=1e-10)

Converged in 4 iterations

log.set_log_level(log.LogLevel.INFO)

Now save the displacement to a file that we can view in Paraview

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

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    pyvista.start_xvfb()
    V = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 1, (geometry.mesh.geometry.dim,)))
    uh = dolfinx.fem.Function(V)
    uh.interpolate(problem.u)
    # Create plotter and pyvista grid
    p = pyvista.Plotter()
    topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

    # Attach vector values to grid and warp grid by vector
    grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
    actor_0 = p.add_mesh(grid, style="wireframe", color="k")
    warped = grid.warp_by_vector("u", factor=1.5)
    actor_1 = p.add_mesh(warped, show_edges=True)
    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        figure_as_array = p.screenshot("problem2.png")

INFO:trame_server.utils.namespace:Translator(prefix=None)

[2025-07-01 08:30:24.935] [info] Checking required entities per dimension
[2025-07-01 08:30:24.935] [info] Cell type: 0 dofmap: 2838x4
[2025-07-01 08:30:24.935] [info] Global index computation
[2025-07-01 08:30:24.935] [info] Got 1 index_maps
[2025-07-01 08:30:24.935] [info] Get global indices

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

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

Finally we can extract the longitudinal position of the endocardial and epicardial apex First we create a function to get the coordinates of the apex

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

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