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 fenicsx_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-03-05 19:41:41 [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 fenicsx_pulse.Geometry object.
geometry = fenicsx_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 = fenicsx_pulse.Guccione(**material_params)
WARNING:fenicsx_pulse.material_models.guccione:Setting C to c_0 kPa
There are now active contraction, so we choose a pure passive model
active_model = fenicsx_pulse.active_model.Passive()
and the model should be incompressible
comp_model = fenicsx_pulse.Incompressible()
and assembles the CardiacModel
model = fenicsx_pulse.CardiacModel(
material=material,
active=active_model,
compressibility=comp_model,
decouple_deviatoric_volumetric=True,
)
We apply a traction in endocardium
traction = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
neumann = fenicsx_pulse.NeumannBC(
traction=fenicsx_pulse.Variable(traction, "kPa"),
marker=geo.markers["ENDO"][0],
)
and finally combine all the boundary conditions
bcs = fenicsx_pulse.BoundaryConditions(neumann=(neumann,))
and create a Mixed problem
problem = fenicsx_pulse.StaticProblem(
model=model, geometry=geometry, bcs=bcs, parameters={"base_bc": fenicsx_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.750527508181 (tol=1e-06), r (rel) = 0.1695466472667266 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 93.06415070374997 (tol=1e-06), r (rel) = 0.003939015828417427 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.06329548858772499 (tol=1e-06), r (rel) = 2.679033006040392e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 5.777026628947666e-08 (tol=1e-06), r (rel) = 2.4451734809306e-12 (tol=1e-10)
Converged in 5 iterations
Solving problem for traction=2.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 30372.591563072612 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 1825.5635880719399 (tol=1e-06), r (rel) = 0.06010562464783162 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 54.12934758791407 (tol=1e-06), r (rel) = 0.0017821774436174629 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.07896978813910814 (tol=1e-06), r (rel) = 2.6000345731156056e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 3.9048948542960816e-07 (tol=1e-06), r (rel) = 1.2856640323849424e-11 (tol=1e-10)
Converged in 5 iterations
Solving problem for traction=3.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 40306.002325679634 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 35134.91654353103 (tol=1e-06), r (rel) = 0.8717043248207719 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 4916.066118701996 (tol=1e-06), r (rel) = 0.12196858619168707 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 227.628357537395 (tol=1e-06), r (rel) = 0.005647505195333379 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.519410268439401 (tol=1e-06), r (rel) = 3.769687343741751e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 0.00011788422453350133 (tol=1e-06), r (rel) = 2.9247312492312168e-09 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 7: r (abs) = 8.6428538550554e-11 (tol=1e-06), r (rel) = 2.1443093723906454e-15 (tol=1e-10)
Converged in 7 iterations
Solving problem for traction=4.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 20422.60652284542 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 8379.660769812626 (tol=1e-06), r (rel) = 0.4103129911668647 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 673.1241547142881 (tol=1e-06), r (rel) = 0.03295975731409742 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 9.112164288803987 (tol=1e-06), r (rel) = 0.0004461802796137119 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.002458977865296588 (tol=1e-06), r (rel) = 1.2040470262921102e-07 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 6: r (abs) = 2.680276835498316e-10 (tol=1e-06), r (rel) = 1.3124068333295592e-14 (tol=1e-10)
Converged in 6 iterations
Solving problem for traction=5.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 15537.378747403713 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 2195.5410477483642 (tol=1e-06), r (rel) = 0.14130704306318323 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 63.98316768407319 (tol=1e-06), r (rel) = 0.004118015575488545 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.08259278772889403 (tol=1e-06), r (rel) = 5.315747853716653e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.638119472179985e-07 (tol=1e-06), r (rel) = 1.0543087729348902e-11 (tol=1e-10)
Converged in 5 iterations
Solving problem for traction=6.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 12303.37944490715 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 860.3455332163771 (tol=1e-06), r (rel) = 0.06992757860300795 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 9.165681268587791 (tol=1e-06), r (rel) = 0.0007449726564664983 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.001537502383423525 (tol=1e-06), r (rel) = 1.2496585920220135e-07 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 3.4147937093681795e-10 (tol=1e-06), r (rel) = 2.775492477216653e-14 (tol=1e-10)
Converged in 5 iterations
Solving problem for traction=7.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 10087.082182729808 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 416.04880397325 (tol=1e-06), r (rel) = 0.041245703805761716 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 1.8941819691494566 (tol=1e-06), r (rel) = 0.00018778294206747955 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 5.81173510192871e-05 (tol=1e-06), r (rel) = 5.761562161037053e-09 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 4.048984321041039e-10 (tol=1e-06), r (rel) = 4.0140292779346483e-14 (tol=1e-10)
Converged in 5 iterations
Solving problem for traction=8.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 8506.388361606121 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 229.37611181174762 (tol=1e-06), r (rel) = 0.026965158662052704 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.5111504096780219 (tol=1e-06), r (rel) = 6.009018022091691e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 3.736647527207666e-06 (tol=1e-06), r (rel) = 4.392754443323037e-10 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 5: r (abs) = 5.12861135246978e-10 (tol=1e-06), r (rel) = 6.029129090341024e-14 (tol=1e-10)
Converged in 5 iterations
Solving problem for traction=9.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 7330.422209848763 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 138.49085203148363 (tol=1e-06), r (rel) = 0.018892616014042785 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.16757857483509547 (tol=1e-06), r (rel) = 2.286069888443068e-05 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 3.55862492632172e-07 (tol=1e-06), r (rel) = 4.854597490360844e-11 (tol=1e-10)
Converged in 4 iterations
Solving problem for traction=10.0
INFO:scifem.solvers:Newton iteration 1: r (abs) = 6424.6876697761745 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 89.41446624488565 (tol=1e-06), r (rel) = 0.013917324987722041 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.0635887498920717 (tol=1e-06), r (rel) = 9.897562832698296e-06 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 4: r (abs) = 4.574304350199894e-08 (tol=1e-06), r (rel) = 7.119885954486026e-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")
[2025-03-05 19:42:28.590] [info] Checking required entities per dimension
[2025-03-05 19:42:28.590] [info] Cell type: 0 dofmap: 2838x4
[2025-03-05 19:42:28.590] [info] Global index computation
[2025-03-05 19:42:28.590] [info] Got 1 index_maps
[2025-03-05 19:42:28.590] [info] Get global indices
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
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 = fenicsx_pulse.utils.evaluate_at_vertex_tag(U, geo.vfun, geo.markers["ENDOPT"][0])
u_endo_apex = fenicsx_pulse.utils.evaluate_at_vertex_tag(problem.u, geo.vfun, geo.markers["ENDOPT"][0])
endo_apex_pos = fenicsx_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 = fenicsx_pulse.utils.evaluate_at_vertex_tag(U, geo.vfun, geo.markers["EPIPT"][0])
u_epi_apex = fenicsx_pulse.utils.evaluate_at_vertex_tag(problem.u, geo.vfun, geo.markers["EPIPT"][0])
epi_apex_pos = fenicsx_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