Mathematical Background & Implementation Details

This document outlines the mathematical theory of finite hyperelasticity used in fenicsx-pulse and demonstrates how these concepts are mapped to specific functions and classes in the library.

We will assume a standard continuum mechanics framework where a body \(\mathcal{B}\) is identified with a reference configuration \(\Omega_0\). The motion is described by the map \(\mathbf{x} = \chi(\mathbf{X}, t)\), where \(\mathbf{X} \in \Omega_0\) is the reference position and \(\mathbf{x}\) is the current position.

from pathlib import Path
import logging
import dolfinx
import ufl
import pulse
from mpi4py import MPI
import numpy as np

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger("pulse")
for lib in ["trame_server", "wslink"]:
    logging.getLogger(lib).setLevel(logging.WARNING)
logger.setLevel(logging.DEBUG)
# dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

1. Geometry and Mesh

The first step in any Finite Element simulation is defining the domain. Here, we create a simple unit cube mesh to serve as our reference configuration \(\Omega_0\). In a realistic cardiac simulation, this would be replaced by a patient-specific geometry.

comm = MPI.COMM_WORLD
mesh = dolfinx.mesh.create_unit_cube(comm, 3, 3, 3)

Visualization

We can visualize the mesh using pyvista.

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:
    p = pyvista.Plotter()
    topology, cell_types, geometry = dolfinx.plot.vtk_mesh(mesh)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
    p.add_mesh(grid, show_edges=True)
    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        # Save screenshot if running in CI/headless
        p.screenshot("maths_mesh.png")

2026-08-14 16:26:41.907 (   0.856s) [    7F97F13BA140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=:99.0

Boundary Markers

To solve the boundary value problem, we need to identify specific parts of the boundary \(\partial \Omega_0\). We define markers using geometric locators:

• Marker 1 (“X0”): The face where \(X=0\) (to be fixed).

• Marker 2 (“X1”): The face where \(X=1\) (to apply traction).

boundaries = [
    pulse.Marker(name="X0", marker=1, dim=2, locator=lambda x: np.isclose(x[0], 0)),
    pulse.Marker(name="X1", marker=2, dim=2, locator=lambda x: np.isclose(x[0], 1)),
]

We wrap the mesh and markers into a pulse.Geometry object. This object manages the integration measures (dx for volume, ds for surface) and ensures they are set up with the correct quadrature degree.

geo = pulse.Geometry(
    mesh=mesh,
    boundaries=boundaries,
    metadata={"quadrature_degree": 4},
)

DEBUG:pulse.geometry:Created Geometry with 2 boundaries

DEBUG:pulse.geometry:Markers: X0, X1

DEBUG:pulse.geometry:Metadata: {'quadrature_degree': 4}

We can also visualize the boundary markers

geo.mesh.topology.create_connectivity(mesh.topology.dim-1 , mesh.topology.dim)
vtk_bmesh = dolfinx.plot.vtk_mesh(geo.mesh, geo.facet_tags.dim, geo.facet_tags.indices)
bgrid = pyvista.UnstructuredGrid(*vtk_bmesh)
bgrid.cell_data["Facet tags"] = geo.facet_tags.values
bgrid.set_active_scalars("Facet tags")
p = pyvista.Plotter(window_size=[800, 800])
p.add_mesh(bgrid, show_edges=True)
p.add_mesh(grid, show_edges=True, style="wireframe", color="k")
if not pyvista.OFF_SCREEN:
    p.show()
else:
    figure = p.screenshot("facet_tags.png")

2. Constitutive Equations

The material behavior is governed by a Strain Energy Density Function \(\Psi(\mathbf{C})\). In fenicsx-pulse, the total energy is composed of three parts:

\[ \Psi = \Psi_{\text{passive}} + \Psi_{\text{active}} + \Psi_{\text{vol}} \]

The stress tensors are derived from \(\Psi\) via automatic differentiation:

• Second Piola-Kirchhoff stress: \(\mathbf{S} = 2 \frac{\partial \Psi}{\partial \mathbf{C}}\)

• First Piola-Kirchhoff stress: \(\mathbf{P} = \mathbf{F} \mathbf{S}\)

A. Passive Material (\(\Psi_{\text{passive}}\))

We use the Holzapfel-Ogden model (pulse.HolzapfelOgden) [HO09], a standard for ventricular myocardium.

\[ \Psi_{HO} = \frac{a}{2b} (e^{b(I_1-3)} - 1) + \sum_{i=f,s} \frac{a_i}{2b_i} \mathcal{H}(I_{4i}-1) (e^{b_i(I_{4i}-1)^2} - 1) + \frac{a_{fs}}{2b_{fs}} (e^{b_{fs}I_{8fs}^2} - 1) \]

Here \(\mathcal{H}(\cdot)\) is the Heaviside function, ensuring fibers only stiffen in tension.

# Retrieve default parameters (which are wrapped in pulse.units.Variable)
material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
print(f"Parameter 'a': {material_params['a']}")

# Define constant fiber/sheet fields for this simple cube
f0 = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((1.0, 0.0, 0.0)))
s0 = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0.0, 1.0, 0.0)))

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

DEBUG:pulse.material_models.holzapfelogden:Created material model: HolzapfelOgden

DEBUG:pulse.material_models.holzapfelogden:Material parameters: {'a': Variable(value=2.28, unit=<Unit('kilogram / meter / second ** 2')>), 'b': Variable(value=9.726, unit=<Unit('dimensionless')>), 'a_f': Variable(value=1.685, unit=<Unit('kilogram / meter / second ** 2')>), 'b_f': Variable(value=15.779, unit=<Unit('dimensionless')>), 'a_s': Variable(value=0.0, unit=<Unit('kilogram / meter / second ** 2')>), 'b_s': Variable(value=0.0, unit=<Unit('dimensionless')>), 'a_fs': Variable(value=0.0, unit=<Unit('kilogram / meter / second ** 2')>), 'b_fs': Variable(value=0.0, unit=<Unit('dimensionless')>)}

Parameter 'a': 2280.0 kilogram / meter / second ** 2 (2.28 kilopascal)

B. Active Contraction (\(\Psi_{\text{active}}\))

Contraction is a change of the constitutive response, not an external load: the myocardium generates stress from within. Two families of models exist for writing that down [AP12]:

• Active strain: decompose the deformation gradient multiplicatively, \(\mathbf{F} = \mathbf{F}_e \mathbf{F}_a\), and evaluate the passive law on the elastic part \(\mathbf{F}_e\) only. In fenicsx-pulse this is what the pulse.ActiveModel.Fe() hook exists for.

• Active stress: add an extra stress (equivalently, an extra strain energy) on top of the passive response. This is what all the active models shipped with fenicsx-pulse do, so Fe is the identity for each of them.

All active models therefore expose the same interface as a material – strain_energy(C), S(C), P(F) – and the total energy is just a sum. Passing pulse.Passive gives a purely passive simulation.

The two active-stress conventions

Given a scalar active tension \(T_a\) and the fiber direction \(\mathbf{f}_0\), there is more than one way to turn \(T_a\) into a stress, and the choice matters numerically. Let

\[ \lambda = \sqrt{I_{4f}} = \sqrt{\mathbf{f}_0 \cdot (\mathbf{C} \mathbf{f}_0)} = \| \mathbf{F} \mathbf{f}_0 \| \]

be the fiber stretch. pulse.ActiveStress takes a formulation argument (pulse.ActiveStressFormulation) selecting between the two conventions:

1. invariant (the default) – the energy is linear in \(I_{4f}\):

\[ \Psi_{\text{active}} = \frac{1}{2} T_a (I_{4f} - 1) \quad \Longrightarrow \quad \mathbf{S}_{\text{active}} = T_a \, \mathbf{f}_0 \otimes \mathbf{f}_0, \quad \mathbf{P}_{\text{active}} = T_a \, \mathbf{F} \mathbf{f}_0 \otimes \mathbf{f}_0 \]

Adding \(T_a \mathbf{f}_0 \otimes \mathbf{f}_0\) to the second Piola-Kirchhoff stress is by far the most widespread form of the active stress approach in cardiac mechanics, which is why it is the default here. Note that the fiber traction it actually delivers scales with the stretch, \(\| \mathbf{P}_{\text{active}} \mathbf{f}_0 \| = T_a \lambda\), so \(T_a\) is only equal to the generated tension in the reference configuration.

2. stretch – the energy is linear in \(\lambda\) instead:

\[ \Psi_{\text{active}} = T_a (\lambda - 1) \quad \Longrightarrow \quad \mathbf{S}_{\text{active}} = \frac{T_a}{\lambda} \, \mathbf{f}_0 \otimes \mathbf{f}_0, \quad \mathbf{P}_{\text{active}} = T_a \, \frac{\mathbf{F} \mathbf{f}_0 \otimes \mathbf{f}_0}{\| \mathbf{F} \mathbf{f}_0 \|} \]

Here \(T_a\) is the first Piola-Kirchhoff fiber traction, independent of \(\lambda\). This is the convention used by Regazzoni and Quarteroni [RQ21].

The two differ by exactly a factor of \(\lambda\), so switching between them changes the answer – it is a modelling decision, not a refactoring. Pick stretch when \(T_a\) comes from a force-generation model whose tension and stiffness are defined against \(\lambda\) (as in [RDedeQ20] or [LPHS+17]); keep invariant otherwise.

Transverse activation

Experimentally, contraction also generates some stress across the fibers. The invariant formulation accepts a parameter \(\eta \in [0, 1]\) that blends the fiber-directed energy towards an isotropic one:

\[ \Psi_{\text{active}} = \frac{1}{2} T_a \left[ (I_{4f} - 1) + \eta \left( (I_1 - 3) - (I_{4f} - 1) \right) \right] \]

so \(\eta = 0\) (the default) puts all active stress along \(\mathbf{f}_0\) and \(\eta = 1\) puts it all in the transverse plane. There is no accepted way to write the same blend for an energy expressed in \(\lambda\), so the stretch formulation requires \(\eta = 0\) and raises NotImplementedError otherwise.

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

DEBUG:pulse.active_stress:Created ActiveStress model with Isotropy: ActiveStressModels.transversely

The same model in the stretch convention would be constructed as

active_model = pulse.ActiveStress(
    f0,
    activation=Ta,
    formulation=pulse.ActiveStressFormulation.stretch,
)

Length-dependent activation (Frank-Starling)

In real myocardium the tension a sarcomere develops depends on how far it has been stretched. pulse.FrankStarlingActiveStress captures this cheaply, by scaling the supplied activation with a piecewise-linear ascending limb \(g(\lambda)\):

\[\begin{split} T_a \mapsto g(\lambda) \, T_a, \qquad g(\lambda) = \begin{cases} a_{\min} & \lambda \le \lambda_{\text{thresh}} \\ a_{\min} + m (\lambda - \lambda_{\text{thresh}}) & \lambda_{\text{thresh}} < \lambda \le \lambda_{\text{opt}} \\ a_{\max} & \lambda > \lambda_{\text{opt}} \end{cases} \end{split}\]

with slope \(m = (a_{\max} - a_{\min}) / (\lambda_{\text{opt}} - \lambda_{\text{thresh}})\). Because \(g\) depends on the unknown displacement, the model must be told which field to read via register(u); the Problem classes do this for you. This is a curve fit, not a mechanism – biophysical cross-bridge models such as [LPHS+17] or [LMCN24] produce length-dependent activation from their own kinetics instead. See the isometric twitch demos for a comparison.

Stabilized active stress

When \(T_a\) is computed by an external force-generation solver – a cell model advanced once per time step, with mechanics then solved at fixed \(T_a\) – the resulting staggered scheme has a failure mode that is easy to hit. Regazzoni and Quarteroni [RQ21] showed that whenever the active stiffness \(K_a\) exceeds the passive stiffness \(K_p\) (routine in contracting myocardium) the scheme is not merely inaccurate but non-convergent: its amplification factor tends to \(-K_a/K_p < -1\) as \(\Delta t \to 0\), so refining the time step makes the oscillations worse rather than better.

The cure is to stop treating the active tension as a dead load over the mechanics solve and restore the fact that the cross-bridges behave as springs. pulse.StabilizedActiveStress implements the resulting consistent stabilization,

\[ \Psi_{\text{active}} = T_a \Delta\lambda + \frac{1}{2} K_a \Delta\lambda^2, \qquad \Delta\lambda = \lambda - \lambda_{\text{prev}} \]

\[ \mathbf{P}_{\text{active}} = \left[ T_a + K_a \Delta\lambda \right] \frac{\mathbf{F} \mathbf{f}_0 \otimes \mathbf{f}_0}{\| \mathbf{F} \mathbf{f}_0 \|} \]

where \(\lambda_{\text{prev}}\) is the fiber stretch at the previous time step. The added term is \(\mathcal{O}(\Delta t)\) and vanishes in the limit, so the scheme remains consistent with the same continuous problem – it is a numerical device, not a change of model – while becoming unconditionally stable. Note that it is built on the stretch convention above, so \(T_a\) and \(K_a = \partial \dot{T_a} / \partial \dot{\lambda}\) must refer to the same kinematic variable.

Using it requires one extra call per time step, in this order: advance the force-generation model using \(\lambda_{\text{prev}}\), assign \(T_a\) and \(K_a\), solve mechanics, then call update_prev(u) so that \(\lambda_{\text{prev}}\) tracks the same stretch that drove the cell model.

C. Compressibility (\(\Psi_{\text{vol}}\))

Myocardium is nearly incompressible (\(J \approx 1\)). We enforce this using the Incompressible model (pulse.Incompressible), which uses a Lagrange multiplier \(p\) (hydrostatic pressure).

\[ \Psi_{\text{vol}} = p (J - 1) \]

comp_model = pulse.Incompressible()

DEBUG:pulse.compressibility:Created Incompressible compressibility model

Assembly: The Cardiac Model

The pulse.CardiacModel class aggregates these components into a single object that provides the total \(\mathbf{S}\) and \(\mathbf{P}\) tensors.

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

DEBUG:pulse.cardiac_model:Created CardiacModel with components:

DEBUG:pulse.cardiac_model:  Material: HolzapfelOgden

DEBUG:pulse.cardiac_model:  Active Model: ActiveStress

DEBUG:pulse.cardiac_model:  Compressibility: Incompressible

DEBUG:pulse.cardiac_model:  Viscoelasticity: NoneViscoElasticity

3. Balance Laws & Boundary Value Problem

We solve the balance of linear momentum in the reference configuration:

\[ \nabla \cdot \mathbf{P} + \rho_0 \mathbf{B} = \mathbf{0} \quad \text{in } \Omega_0 \]

The weak (variational) form used in pulse.StaticProblem is derived by multiplying by a test function \(\delta \mathbf{u}\) and integrating by parts:

\[ \int_{\Omega_0} \mathbf{P} : \nabla \delta \mathbf{u} \, \text{d}X - \int_{\partial \Omega_N} \mathbf{t} \cdot \delta \mathbf{u} \, \text{d}S = 0 \]

For the incompressible case, we also add the constraint equation: $\( \int_{\Omega_0} (J - 1) \delta p \, \text{d}X = 0 \)$

Boundary Conditions

We define the specific conditions for our cube:

1. Dirichlet BC: Fix displacement at \(X=0\).

2. Neumann BC: Apply traction at \(X=1\).

# 1. Dirichlet: Fix X0
def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
    facets = geo.facet_tags.find(1)
    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
    dofs = dolfinx.fem.locate_dofs_topological(V, 2, facets)
    u_fixed = dolfinx.fem.Function(V)
    u_fixed.x.array[:] = 0.0
    return [dolfinx.fem.dirichletbc(u_fixed, dofs)]

# 2. Neumann: Traction on X1 (Marker 2)
traction = pulse.Variable(dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(1.0)), "kPa")
neumann_bc = pulse.NeumannBC(traction=traction, marker=2)

# Collect BCs
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))

DEBUG:pulse.boundary_conditions:Created NeumannBC on marker 2 with traction 1000.0 * c_5 kilogram / meter / second ** 2 (c_5 kilopascal)

4. Solving the Problem

The StaticProblem class handles the assembly of the mixed function space (for \(\mathbf{u}\) and \(p\)), the construction of the variational forms, and the Newton solver configuration.

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

# Apply active tension to simulate contraction
Ta.value = 2.0

# Solve the system
problem.solve()

DEBUG:pulse.problem:Initializing function spaces...

DEBUG:pulse.problem:Initializing displacement function space...

DEBUG:pulse.problem:Displacement space: family=P, degree=2

DEBUG:pulse.problem:Initializing pressure function space...

DEBUG:pulse.problem:Model is incompressible, initializing pressure space with Lagrange multiplier

DEBUG:pulse.problem:Initializing cavity pressure function spaces...

DEBUG:pulse.problem:No cavity pressure states needed

DEBUG:pulse.problem:Initializing rigid body function space...

DEBUG:pulse.problem:No rigid body constraint needed

DEBUG:pulse.problem:Initializing ufl forms...

DEBUG:pulse.problem:Creating material form...

DEBUG:pulse.problem:Creating cavity pressure form...

DEBUG:pulse.problem:Creating Neumann boundary condition form...

DEBUG:pulse.problem:Creating Newton solver...

DEBUG:pulse.problem:Initialized StaticProblem with parameters:

DEBUG:pulse.problem:  u_space: P_2

DEBUG:pulse.problem:  p_space: P_1

DEBUG:pulse.problem:  base_bc: BaseBC.free

DEBUG:pulse.problem:  rigid_body_constraint: False

DEBUG:pulse.problem:  mesh_unit: m

DEBUG:pulse.problem:  base_marker: BASE

DEBUG:pulse.problem:  petsc_options: {'ksp_type': 'preonly', 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps', 'snes_error_if_not_converged': True, 'ksp_error_if_not_converged': True, 'snes_type': 'newtonls', 'snes_atol': 1e-06, 'snes_rtol': 1e-10, 'snes_stol': 1e-08, 'snes_max_it': 50, 'snes_linesearch_type': 'l2'}

DEBUG:pulse.problem:Number of cavities: 0

DEBUG:pulse.problem:Boundary conditions: BoundaryConditions(neumann=(NeumannBC(traction=Variable(value=Constant(Mesh(blocked element (Basix element (P, tetrahedron, 1, gll_warped, unset, False, float64, []), (3,)), 0), (), 5), unit=<Unit('kilogram / meter / second ** 2')>), marker=2),), dirichlet=(<function dirichlet_bc at 0x7f9774ad2700>,), robin=(), body_force=())

DEBUG:pulse.problem:Solving the system...

DEBUG:pulse.problem:Updating old states to current values...

DEBUG:pulse.problem:Updating old displacement state to current value...

DEBUG:pulse.problem:Updating old pressure state to current value...

DEBUG:pulse.problem:Solved in 4 iterations, converged: True

True

Visualization of Result

Finally, we can visualize the deformed configuration.

try:
    import pyvista
except ImportError:
    pass
else:
    # Interpolate solution to a standard space for plotting
    p = pyvista.Plotter()
    topology, cell_types, geometry = dolfinx.plot.vtk_mesh(problem.u_space)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

    # Add reference mesh (wireframe)
    p.add_mesh(grid, style="wireframe", color="black", opacity=0.3, label="Reference")

    # Warp by displacement
    grid["u"] = problem.u.x.array.reshape((-1, 3))
    warped = grid.warp_by_vector("u", factor=1.0)

    # Add deformed mesh
    p.add_mesh(warped, show_edges=True, label="Deformed")
    p.add_legend()
    p.show_axes()

    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot("maths_deformed.png")

5. Kinematics definitions

The primary unknown in our problem is the Displacement field \(\mathbf{u}(\mathbf{X})\). We define a standard Lagrange function space and the function \(\mathbf{u}\).

The Deformation Gradient \(\mathbf{F}\) is defined as:

\[ \mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \mathbf{I} + \nabla \mathbf{u} \]

In fenicsx-pulse, this is computed via pulse.kinematics.DeformationGradient().

u = problem.u
F = pulse.kinematics.DeformationGradient(u)
print(f"Shape of F: {F.ufl_shape}")

Shape of F: (3, 3)

The volume change is measured by the Jacobian \(J = \det \mathbf{F}\).

J = pulse.kinematics.Jacobian(F)

Since this is an incompressible problem we expect \(J\) to be equal to 1.0

mesh_volume = comm.allreduce(dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(mesh, 1.0) * geo.dx)), op=MPI.SUM)
comm.allreduce(dolfinx.fem.assemble_scalar(dolfinx.fem.form(J * geo.dx)), op=MPI.SUM) / mesh_volume

1.0000000000004865

Here we first compile the form

dolfinx.fem.form(J * geo.dx)

then we assemble the form

dolfinx.fem.assemble_scalar(dolfinx.fem.form(J * geo.dx))

which will assemble the form locally on each process, and finally we perform an allreduce using the MPI communicator to sum up the contributions from all the processors. Note that we also divide by the volume (which is 1.0 in the case of the Unit Cube) which is computed in the same fashion.

Strain Tensors

To define material laws that are independent of rigid body rotations, we use strain tensors derived from \(\mathbf{F}\). pulse.kinematics provides standard tensors:

1. Right Cauchy-Green tensor: \(\mathbf{C} = \mathbf{F}^T \mathbf{F}\)

2. Left Cauchy-Green tensor: \(\mathbf{B} = \mathbf{F} \mathbf{F}^T\)

3. Green-Lagrange strain: \(\mathbf{E} = \frac{1}{2}(\mathbf{C} - \mathbf{I})\)

C = pulse.kinematics.RightCauchyGreen(F)
B = pulse.kinematics.LeftCauchyGreen(F)
E = pulse.kinematics.GreenLagrangeStrain(F)

Now, say you are interested in the strain in a given direction, e.g the fiber strain (\(E_{ff}\)). Then one can get the by computing the inner product

Eff_ufl_expr = ufl.inner(E * f0, f0)

If we now would like to visualize this in Pyvista or Paraview then we need to first interpolate this into a function space. Since \(\mathbf{u}\) is \(\mathbb{P}_2\), i.e. a second order polynomial that is continuous but not continuously differentiable across elements, the gradient \(\nabla \mathbf{u}\) belongs to a discontinuous first order space. Since \(\mathbf{E}\) is a function of \(\nabla \mathbf{u}^T \nabla \mathbf{u}\), a reasonable space would be a second order discontinuous space.

V_strain = dolfinx.fem.functionspace(mesh, ("DG", 2))
Eff = dolfinx.fem.Function(V_strain)

We can now interpolate the fiber strain into this space

Eff_expr = dolfinx.fem.Expression(Eff_ufl_expr, V_strain.element.interpolation_points)
Eff.interpolate(Eff_expr)

We can now visualize the strain

try:
    import pyvista
except ImportError:
    pass
else:
    # Interpolate solution to a standard space for plotting
    p = pyvista.Plotter()
    topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V_strain)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

    grid["Eff"] = Eff.x.array
    # Add reference mesh (wireframe)
    p.add_mesh(grid, show_edges=True, cmap="inferno")
    p.show_axes()

    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot("Eff.png")

6. Material Invariants

Hyperelastic constitutive laws are typically expressed in terms of the invariants of \(\mathbf{C}\). pulse.invariants provides helper functions for these.

Isotropic Invariants

For isotropic materials, the strain energy \(\Psi\) depends on:

\[ I_1 = \text{tr}(\mathbf{C}), \quad I_2 = \frac{1}{2}(I_1^2 - \text{tr}(\mathbf{C}^2)), \quad I_3 = \det \mathbf{C} = J^2 \]

I1 = pulse.invariants.I1(C)
I2 = pulse.invariants.I2(C)
I3 = pulse.invariants.I3(C)

Anisotropic Invariants

Cardiac tissue is orthotropic. Its behavior depends on the local microstructure defined by fiber (\(\mathbf{f}_0\)), sheet (\(\mathbf{s}_0\)), and normal (\(\mathbf{n}_0\)) directions.

We define pseudo-invariants to capture stretch along these directions and shear between them:

\[ I_{4f} = \mathbf{f}_0 \cdot (\mathbf{C} \mathbf{f}_0) \quad (\text{Fiber stretch squared}) \]

\[ I_{8fs} = \mathbf{f}_0 \cdot (\mathbf{C} \mathbf{s}_0) \quad (\text{Fiber-sheet shear coupling}) \]

I4f = pulse.invariants.I4(C, f0)
I8fs = pulse.invariants.I8(C, f0, s0)

7. Stress tensors

We can also compute the stress tensors directly from the cardiac model. Note that we need to wrap the deformation gradient and right Cauchy-Green tensor into a ufl.variable in order to be able to differentiate the strain energy function.

S = model.S(ufl.variable(C))
P = model.P(ufl.variable(F))
T = model.sigma(ufl.variable(F))

Similar to the strain case, one might also be interested in visualizing the fiber stress

\[\sigma_{ff} = \mathbf{f} \cdot (\mathbf{\sigma} \mathbf{f}),\]

where

\[\mathbf{f} = \frac{\mathbf{F} \mathbf{f}_0}{\| \mathbf{F} \mathbf{f}_0 \|} \]

f = F * f0
f /= ufl.sqrt(f**2)
Tff_ufl_expr = ufl.inner(T * f, f)
Tff = dolfinx.fem.Function(V_strain)
Tff_expr = dolfinx.fem.Expression(Tff_ufl_expr, V_strain.element.interpolation_points)
Tff.interpolate(Tff_expr)
try:
    import pyvista
except ImportError:
    pass
else:
    # Interpolate solution to a standard space for plotting
    p = pyvista.Plotter()
    topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V_strain)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

    grid["Tff"] = Tff.x.array
    # Add reference mesh (wireframe)
    p.add_mesh(grid, show_edges=True, cmap="inferno", clim=(0, 1000))
    p.show_axes()

    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        p.screenshot("Eff.png")

Learn more

To learn more you could check out Holzapfel’s Nonlinear continuum mechanics book [Hol00]

References

[AP12]

Davide Ambrosi and Simone Pezzuto. Active stress vs. active strain in mechanobiology: constitutive issues. Journal of Elasticity, 107(2):199–212, 2012. doi:10.1007/s10659-011-9351-4.

[Hol00]

Gerhard A Holzapfel. Nonlinear solid mechanics: a continuum approach for engineering. John Wiley & Sons, Chichester, 2000. ISBN 9780471823193.

[HO09]

Gerhard A Holzapfel and Ray W Ogden. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 367(1902):3445–3475, 2009. doi:10.1098/rsta.2009.0091.

[LPHS+17] (1,2)

Sander Land, So-Jin Park-Holohan, Nicolas P Smith, Cristobal G Dos Remedios, Jonathan C Kentish, and Steven A Niederer. A model of cardiac contraction based on novel measurements of tension development in human cardiomyocytes. Journal of molecular and cellular cardiology, 106:68–83, 2017. doi:10.1016/j.yjmcc.2017.03.008.

[LMCN24]

Alexandre Lewalle, Gregory Milburn, Kenneth S Campbell, and Steven A Niederer. Cardiac length-dependent activation driven by force-dependent thick-filament dynamics. Biophysical Journal, 123(18):2996–3009, 2024. doi:10.1016/j.bpj.2024.05.025.

[RDedeQ20]

Francesco Regazzoni, Luca Dedè, and Alfio Quarteroni. Biophysically detailed mathematical models of multiscale cardiac active mechanics. PLOS Computational Biology, 16(10):e1008294, 2020. doi:10.1371/journal.pcbi.1008294.

[RQ21] (1,2)

Francesco Regazzoni and Alfio Quarteroni. An oscillation-free fully staggered algorithm for velocity-dependent active models of cardiac mechanics. Computer Methods in Applied Mechanics and Engineering, 373:113506, 2021. doi:10.1016/j.cma.2020.113506.