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:
Apply Robin Boundary Conditions (springs) to the epicardium and base to mimic surrounding tissue support.
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.
Generate Mesh
if not geodir.exists():
cardiac_geometries.mesh.lv_ellipsoid(
outdir=geodir,
create_fibers=True,
fiber_space="P_2",
)
2026-02-27 16:26:27 [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'
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})
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.
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,
)
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
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[15], line 5
3 print(f"Solving for pressure: {plv} kPa")
4 traction.assign(plv)
----> 5 problem.solve()
6 vtx.write(float(i))
File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:573, in StaticProblem.solve(self, rtol, atol, beta, update_old_states)
571 if update_old_states:
572 self.update_old_states()
--> 573 ret = self._solver.solve(rtol=rtol, atol=atol, beta=beta)
574 self.update_fields()
576 return ret
File /dolfinx-env/lib/python3.12/site-packages/scifem/solvers.py:253, in NewtonSolver.solve(self, atol, rtol, beta)
251 if self._error_on_convergence:
252 if (status := self._solver.getConvergedReason()) <= 0:
--> 253 raise RuntimeError(f"Linear solver did not converge, got reason: {status}")
255 # Update solution
256 offset_start = 0
RuntimeError: Linear solver did not converge, got reason: -11
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")
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))
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")