Problem 3: inflation and active contraction of a ventricle

Problem 3: inflation and active contraction of a ventricle#

In the third problem we will solve the inflation and active contraction of a ventricle. First we import the necessary libraries

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

Next we will create the geometry and save it in the folder called lv_ellipsoid. Now we will also generate fibers and use a sixth order quadrature space for the fibers

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),
        fiber_space="Quadrature_6",
        create_fibers=True,
        fiber_angle_epi=-90,
        fiber_angle_endo=90,
        comm=comm,
    )
    print("Done creating geometry.")

2025-03-05 19:42:31 [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'

Done creating geometry.

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

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

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

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

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

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)),
}
material = fenicsx_pulse.Guccione(f0=geo.f0, s0=geo.s0, n0=geo.n0, **material_params)

Setting C to c_3 kPa

We use an active stress approach with 60% transverse active stress

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

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

1

Now we will solve the problem for a range of active contraction and traction values

target_pressure = 15.0
target_Ta = 60.0

Let us just gradually increase the active contraction and traction with a linear ramp of 40 steps

N = 40

for Ta_value, traction_value in zip(np.linspace(0, target_Ta, N), np.linspace(0, target_pressure, N)):
    print(f"Solving problem for traction={traction_value} and active contraction={Ta_value}")
    Ta.value = Ta_value
    traction.value = traction_value
    problem.solve()

Solving problem for traction=0.0 and active contraction=0.0

Solving problem for traction=0.38461538461538464 and active contraction=1.5384615384615385

Solving problem for traction=0.7692307692307693 and active contraction=3.076923076923077

Solving problem for traction=1.153846153846154 and active contraction=4.615384615384616

Solving problem for traction=1.5384615384615385 and active contraction=6.153846153846154

Solving problem for traction=1.9230769230769231 and active contraction=7.6923076923076925

Solving problem for traction=2.307692307692308 and active contraction=9.230769230769232

Solving problem for traction=2.6923076923076925 and active contraction=10.76923076923077

Solving problem for traction=3.076923076923077 and active contraction=12.307692307692308

Solving problem for traction=3.4615384615384617 and active contraction=13.846153846153847

Solving problem for traction=3.8461538461538463 and active contraction=15.384615384615385

Solving problem for traction=4.230769230769231 and active contraction=16.923076923076923

Solving problem for traction=4.615384615384616 and active contraction=18.461538461538463

Solving problem for traction=5.0 and active contraction=20.0

Solving problem for traction=5.384615384615385 and active contraction=21.53846153846154

Solving problem for traction=5.769230769230769 and active contraction=23.076923076923077

Solving problem for traction=6.153846153846154 and active contraction=24.615384615384617

Solving problem for traction=6.538461538461539 and active contraction=26.153846153846157

Solving problem for traction=6.923076923076923 and active contraction=27.692307692307693

Solving problem for traction=7.307692307692308 and active contraction=29.230769230769234

Solving problem for traction=7.6923076923076925 and active contraction=30.76923076923077

Solving problem for traction=8.076923076923077 and active contraction=32.30769230769231

Solving problem for traction=8.461538461538462 and active contraction=33.84615384615385

Solving problem for traction=8.846153846153847 and active contraction=35.38461538461539

Solving problem for traction=9.230769230769232 and active contraction=36.92307692307693

Solving problem for traction=9.615384615384617 and active contraction=38.46153846153847

Solving problem for traction=10.0 and active contraction=40.0

Solving problem for traction=10.384615384615385 and active contraction=41.53846153846154

Solving problem for traction=10.76923076923077 and active contraction=43.07692307692308

Solving problem for traction=11.153846153846155 and active contraction=44.61538461538462

Solving problem for traction=11.538461538461538 and active contraction=46.15384615384615

Solving problem for traction=11.923076923076923 and active contraction=47.69230769230769

Solving problem for traction=12.307692307692308 and active contraction=49.23076923076923

Solving problem for traction=12.692307692307693 and active contraction=50.769230769230774

Solving problem for traction=13.076923076923078 and active contraction=52.307692307692314

Solving problem for traction=13.461538461538462 and active contraction=53.84615384615385

Solving problem for traction=13.846153846153847 and active contraction=55.38461538461539

Solving problem for traction=14.230769230769232 and active contraction=56.92307692307693

Solving problem for traction=14.615384615384617 and active contraction=58.46153846153847

Solving problem for traction=15.0 and active contraction=60.0

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, "problem3.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:45:49.577] [info] Checking required entities per dimension
[2025-03-05 19:45:49.578] [info] Cell type: 0 dofmap: 2838x4
[2025-03-05 19:45:49.578] [info] Global index computation
[2025-03-05 19:45:49.578] [info] Got 1 index_maps
[2025-03-05 19:45:49.578] [info] Get global indices

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

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: -11.932800 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: -15.140326 mm