LV Ellipsoid

In this demo we will show how to simulate an idealized left ventricular geometry. For this we will use cardiac-geometries to generate an idealized LV ellipsoid.

First we will make the necessary imports.

from pathlib import Path
from mpi4py import MPI
import dolfinx
from dolfinx import log
import fenicsx_pulse
import cardiac_geometries
import cardiac_geometries.geometry

Next we will create the geometry and save it in the folder called lv_ellipsoid. We also make sure to generate fibers which can be done analytically and use a second order Lagrange space for the fibers

outdir = Path("lv_ellipsoid")
outdir.mkdir(parents=True, exist_ok=True)
geodir = outdir / "geometry"
if not geodir.exists():
    cardiac_geometries.mesh.lv_ellipsoid(outdir=geodir, create_fibers=True, fiber_space="P_2")

2025-03-05 19:47:51 [debug    ] Convert file lv_ellipsoid/geometry/lv_ellipsoid.msh to dolfin

Info    : Reading 'lv_ellipsoid/geometry/lv_ellipsoid.msh'...
Info    : 53 entities
Info    : 771 nodes
Info    : 4026 elements
Info    : Done reading 'lv_ellipsoid/geometry/lv_ellipsoid.msh'

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

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

In order to use the geometry with pulse we need to convert it to a fenicsx_pulse.Geometry object. We can do this by using the from_cardiac_geometries method. We also specify that we want to use a quadrature degree of 4

geometry = fenicsx_pulse.Geometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 4})

Next we create the material object, and we will use the transversely isotropic version of the Holzapfel Ogden model

material_params = fenicsx_pulse.HolzapfelOgden.transversely_isotropic_parameters()
material = fenicsx_pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)  # type: ignore

We use an active stress approach with 30% transverse active stress (see fenicsx_pulse.active_stress.transversely_active_stress())

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

We use an incompressible model

comp_model = fenicsx_pulse.Incompressible()

and assembles the CardiacModel

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

We apply a traction in endocardium

traction = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
neumann = fenicsx_pulse.NeumannBC(traction=traction, marker=geometry.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()

8

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

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

i = 1
for plv in [0.1]: #, 0.5, 1.0]:
    print(f"plv: {plv}")
    traction.value = plv
    problem.solve()

    vtx.write(float(i))
    i += 1

plv: 0.1

Now plot with pyvista

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(outdir / "lv_ellipsoid_pressure.png")

[2025-03-05 19:48:03.612] [info] Checking required entities per dimension
[2025-03-05 19:48:03.612] [info] Cell type: 0 dofmap: 2838x4
[2025-03-05 19:48:03.613] [info] Global index computation
[2025-03-05 19:48:03.613] [info] Got 1 index_maps
[2025-03-05 19:48:03.613] [info] Get global indices

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

for ta in [0.1]: #, 0.5, 1.0]:
    print(f"ta: {ta}")
    Ta.value = ta
    problem.solve()
    vtx.write(float(i))
    i += 1

ta: 0.1

log.set_log_level(log.LogLevel.WARNING)
vtx.close()

try:
    import pyvista
except ImportError:
    pass
else:
    # 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(outdir / "lv_ellipsoid_active.png")