LV ellipsoid

In this example we will create a simple ellipsoid model of the left ventricle (LV) of the heart.

from pathlib import Path
from mpi4py import MPI
import pyvista
import dolfinx
import numpy as np
import cardiac_geometries

geodir = Path("lv_ellipsoid")
if not geodir.exists():
    cardiac_geometries.mesh.lv_ellipsoid(outdir=geodir, create_fibers=True, fiber_space="P_1")

2025-10-06 08:42:35 [debug    ] Convert file lv_ellipsoid/lv_ellipsoid.msh to dolfin

Info    : Reading 'lv_ellipsoid/lv_ellipsoid.msh'...
Info    : 53 entities
Info    : 771 nodes
Info    : 4026 elements
Info    : Done reading 'lv_ellipsoid/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,
)

Next we will use pyvista to plot the mesh

vtk_mesh = dolfinx.plot.vtk_mesh(geo.mesh, geo.mesh.topology.dim)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("lv_mesh.png")

2025-10-06 08:42:35.899 (   1.607s) [    7F0C291FE140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

The facets of the mesh are also marked with a MeshTags object. We can access the facet tags with geo.ffun.values and the markers with geo.markers.

print(geo.markers)

{'ENDOPT': [1, 0], 'EPIPT': [2, 0], 'EPIRING': [3, 1], 'ENDORING': [4, 1], 'BASE': [5, 2], 'ENDO': [6, 2], 'EPI': [7, 2], 'MYOCARDIUM': [8, 3]}

Next lets plot the facet tags with pyvista.

assert geo.ffun is not None
vtk_bmesh = dolfinx.plot.vtk_mesh(geo.mesh, geo.ffun.dim, geo.ffun.indices)
bgrid = pyvista.UnstructuredGrid(*vtk_bmesh)
bgrid.cell_data["Facet tags"] = geo.ffun.values
bgrid.set_active_scalars("Facet tags")
p = pyvista.Plotter(window_size=[800, 800])
p.add_mesh(bgrid, show_edges=True)
if not pyvista.OFF_SCREEN:
    p.show()
else:
    figure = p.screenshot("facet_tags.png")

# Now let us look at the fibers of the LV ellipsoid. The fibers are stored in the `f0` attribute of the `Geometry` object.
assert geo.f0 is not None
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(geo.f0.function_space)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, : len(geo.f0)] = geo.f0.x.array.real.reshape((geometry.shape[0], len(geo.f0)))
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=2.0)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe", color="r")
plotter.add_mesh(glyphs)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    fig_as_array = plotter.screenshot("fiber.png")

Now we could try to redo the exercise but with create a thicker ellipsoid with fibers oriented more longitudinally at the endocardium and more circumferentially at the epicardium.

geodir2 = Path("lv_ellipsoid2")
if not geodir2.exists():
    cardiac_geometries.mesh.lv_ellipsoid(
        outdir=geodir2,
        create_fibers=True,
        fiber_space="P_1",
        fiber_angle_endo=80,
        fiber_angle_epi=-30,
        r_short_endo=5.0,
        r_short_epi=10.0,
    )

2025-10-06 08:42:37 [debug    ] Convert file lv_ellipsoid2/lv_ellipsoid.msh to dolfin

Info    : Reading 'lv_ellipsoid2/lv_ellipsoid.msh'...
Info    : 53 entities
Info    : 871 nodes
Info    : 4583 elements
Info    : Done reading 'lv_ellipsoid2/lv_ellipsoid.msh'

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

geo2 = cardiac_geometries.geometry.Geometry.from_folder(
    comm=MPI.COMM_WORLD,
    folder=geodir2,
)

Next we will use pyvista to plot the mesh

vtk_mesh = dolfinx.plot.vtk_mesh(geo2.mesh, geo2.mesh.topology.dim)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("lv_mesh2.png")

assert geo2.ffun is not None
vtk_bmesh = dolfinx.plot.vtk_mesh(geo2.mesh, geo2.ffun.dim, geo2.ffun.indices)
bgrid = pyvista.UnstructuredGrid(*vtk_bmesh)
bgrid.cell_data["Facet tags"] = geo2.ffun.values
bgrid.set_active_scalars("Facet tags")
p = pyvista.Plotter(window_size=[800, 800])
p.add_mesh(bgrid, show_edges=True)
if not pyvista.OFF_SCREEN:
    p.show()
else:
    figure = p.screenshot("facet_tags2.png")

assert geo2.f0 is not None
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(geo2.f0.function_space)
values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
values[:, : len(geo2.f0)] = geo2.f0.x.array.real.reshape((geometry.shape[0], len(geo2.f0)))
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=2.0)
grid = pyvista.UnstructuredGrid(*vtk_mesh)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe", color="r")
plotter.add_mesh(glyphs)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    fig_as_array = plotter.screenshot("fiber2.png")