Left Ventricular Ellipsoid with Custom Boundary Conditions

Left Ventricular Ellipsoid with Custom Boundary Conditions#

This demo builds upon the Left Ventricular Ellipsoid Simulation.

While the previous example utilized the convenience parameter base_bc=pulse.BaseBC.fixed to fully clamp the base, this example demonstrates how to:

  1. Apply Robin Boundary Conditions (springs) to the epicardium and base to mimic surrounding tissue support.

  2. Manually define Dirichlet Boundary Conditions to constrain specific degrees of freedom.

The geometry generation, constitutive modeling, and active stress implementation remain identical to the base example.

Imports#

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

Geometry and Materials#

We generate the same truncated ellipsoid geometry and define the material model (Holzapfel-Ogden with Active Stress) as in the base example.

  1. Generate Mesh

outdir = Path("lv_ellipsoid_custom_bcs")
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-12-31 10:14:03 [debug    ] Convert file lv_ellipsoid_custom_bcs/geometry/lv_ellipsoid.msh to dolfin

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

  1. Load Geometry

geo = cardiac_geometries.geometry.Geometry.from_folder(
    comm=MPI.COMM_WORLD,
    folder=geodir,
)
geometry = pulse.Geometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 4})

  1. Define Constitutive Model

material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
material = pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)

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

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

Custom Boundary Conditions#

Instead of using the preset base_bc parameter, we will explicitly define our boundary conditions.

1. Neumann BC (Cavity Pressure)#

Standard pressure load on the endocardium.

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

2. Robin BCs (Elastic Support)#

We apply a spring-like penalty on the Epicardium and the Base to represent the resistance provided by the pericardium and surrounding tissue.

\[ \mathbf{P}\mathbf{N} + k \mathbf{u} = 0 \]
robin_epi = pulse.RobinBC(
    value=pulse.Variable(
        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
        "Pa / m",
    ),
    marker=geometry.markers["EPI"][0],
)

robin_base = pulse.RobinBC(
    value=pulse.Variable(
        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
        "Pa / m",
    ),
    marker=geometry.markers["BASE"][0],
)

3. Manual Dirichlet BC#

We manually define a Dirichlet condition. In this specific case, we constrain the displacement in the x-direction (\(u_x = 0\)) on the Base. Ideally, the Base plane is perpendicular to the X-axis in this geometry, so this prevents the base from moving longitudinally while allowing expansion/sliding in the Y-Z plane (resisted only by the Robin spring).

def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
    # Find facets for the BASE marker
    facets = geo.ffun.find(geo.markers["BASE"][0])
    # Locate degrees of freedom for the x-component (sub(0)) on these facets
    dofs = dolfinx.fem.locate_dofs_topological(V.sub(0), 2, facets)
    # Return the Dirichlet BC object
    return [dolfinx.fem.dirichletbc(0.0, dofs, V.sub(0))]

# We collect all conditions into the `BoundaryConditions` container.
bcs = pulse.BoundaryConditions(
    neumann=(neumann,),
    dirichlet=(dirichlet_bc,),
    robin=(robin_epi, robin_base),
)

Solving the Problem#

We initialize the StaticProblem. Note: We do not pass parameters={"base_bc": ...} here, as we are fully controlling the boundaries via the bcs argument.

problem = pulse.StaticProblem(
    model=model,
    geometry=geometry,
    bcs=bcs,
)

log.set_log_level(log.LogLevel.INFO)

Phase 1: Passive Inflation#

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

pressures = [1.0] # kPa
for i, plv in enumerate(pressures, start=1):
    print(f"Solving for pressure: {plv} kPa")
    traction.assign(plv)
    problem.solve()
    vtx.write(float(i))

Solving for pressure: 1.0 kPa

Visualization#

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    V = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 1, (geometry.mesh.geometry.dim,)))
    uh = dolfinx.fem.Function(V)
    uh.interpolate(problem.u)

    p = pyvista.Plotter()
    topology, cell_types, geometry_data = dolfinx.plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_data)
    grid["u"] = uh.x.array.reshape((geometry_data.shape[0], 3))

    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")
    warped = grid.warp_by_vector("u", factor=1.0)
    p.add_mesh(warped, show_edges=True, color="firebrick", label="Inflated")

    p.add_legend()
    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot(outdir / "lv_ellipsoid_pressure.png")

[2025-12-31 10:14:20.122] [info] Checking required entities per dimension
[2025-12-31 10:14:20.122] [info] Cell type: 0 dofmap: 2838x4
[2025-12-31 10:14:20.122] [info] Global index computation
[2025-12-31 10:14:20.122] [info] Got 1 index_maps
[2025-12-31 10:14:20.122] [info] Get global indices

2025-12-31 10:14:20.826 (   0.944s) [    7FD8D92F5140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=:99.0

Phase 2: Active Contraction#

active_tensions = [3.0] # kPa
for i, ta in enumerate(active_tensions, start=len(pressures) + 1):
    print(f"Solving for active tension: {ta} kPa")
    Ta.assign(ta)
    problem.solve()
    vtx.write(float(i))

Solving for active tension: 3.0 kPa

vtx.close()

try:
    import pyvista
except ImportError:
    pass
else:
    uh.interpolate(problem.u)
    grid["u"] = uh.x.array.reshape((geometry_data.shape[0], 3))

    p = pyvista.Plotter()
    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")
    warped = grid.warp_by_vector("u", factor=1.0)
    p.add_mesh(warped, show_edges=False, color="red", label="Contracted")

    p.add_legend()
    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot(outdir / "lv_ellipsoid_active.png")
Contents