Mathematical Background & Implementation Details

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")
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)

[2025-12-15 07:18:00.829] [info] Extract basic topology: 648->648
[2025-12-15 07:18:00.829] [info] Build local dual graphs, re-order cells, and compute process boundary vertices.
[2025-12-15 07:18:00.829] [info] Build local part of mesh dual graph (mixed)
[2025-12-15 07:18:00.829] [info] GPS pseudo-diameter:(18) 148-12
[2025-12-15 07:18:00.830] [info] Create topology (generalised)
[2025-12-15 07:18:00.830] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 1
[2025-12-15 07:18:00.830] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
[2025-12-15 07:18:00.830] [info] Compute ghost indices
[2025-12-15 07:18:00.830] [info] Computing communication graph edges (using PCX algorithm). Number of input edges: 0
[2025-12-15 07:18:00.830] [info] Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
[2025-12-15 07:18:00.830] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:00.830] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:00.830] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:00.830] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:00.830] [info] Number of neighbourhood source ranks in distribute_to_postoffice: 0
[2025-12-15 07:18:00.830] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:00.830] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:00.830] [info] Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
[2025-12-15 07:18:00.830] [info] Create Geometry (multiple)
[2025-12-15 07:18:00.830] [info] Got 1 dof layouts
[2025-12-15 07:18:00.830] [info] Checking required entities per dimension
[2025-12-15 07:18:00.830] [info] Cell type:

 0 dofmap: 162x4
[2025-12-15 07:18:00.830] [info] Global index computation
[2025-12-15 07:18:00.830] [info] Got 1 index_maps
[2025-12-15 07:18:00.830] [info] Get global indices
[2025-12-15 07:18:00.830] [info] Calling compute_local_to_global
[2025-12-15 07:18:00.830] [info] xdofs.size = 648
[2025-12-15 07:18:00.830] [info] dofmap sizes = 648 
[2025-12-15 07:18:00.830] [info] all_dofmaps.size = 648
[2025-12-15 07:18:00.830] [info] nodes.size = 64
[2025-12-15 07:18:00.830] [info] Creating geometry with 1 dofmaps

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")

2025-12-15 07:18:01.710 (   0.859s) [    7F1D6380A140]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=:99.0
INFO:trame_server.utils.namespace:Translator(prefix=None)

INFO:wslink.backends.aiohttp:awaiting runner setup

INFO:wslink.backends.aiohttp:awaiting site startup

INFO:wslink.backends.aiohttp:Print WSLINK_READY_MSG

INFO:wslink.backends.aiohttp:Schedule auto shutdown with timout 0

INFO:wslink.backends.aiohttp:awaiting running future

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}

[2025-12-15 07:18:02.171] [info] Computing mesh entities of dimension 2
[2025-12-15 07:18:02.171] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:02.171] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:02.171] [info] Requesting connectivity (2, 0) - (0, 0)
[2025-12-15 07:18:02.171] [info] Requesting connectivity (2, 0) - (0, 0)
[2025-12-15 07:18:02.171] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:02.171] [info] Computing mesh connectivity 2-3 from transpose.

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")

[2025-12-15 07:18:02.179] [info] Requesting connectivity (2, 0) - (3, 0)

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), 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}}$)#

We model contraction using the Active Stress approach (pulse.ActiveStress). An active component is added to the energy:

\[ \Psi_{\text{active}} = \frac{1}{2} T_a (I_{4f} - 1) \]

Differentiating this yields the active stress contribution $\mathbf{S}_{\text{active}} = T_a \mathbf{f}_0 \otimes \mathbf{f}_0$.

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

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:

Dirichlet BC: Fix displacement at $X=0$.
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...

[2025-12-15 07:18:02.771] [info] Computing mesh entities of dimension 1
[2025-12-15 07:18:02.771] [info] Computing communication graph edges (using NBX algorithm). Number of input edges: 0
[2025-12-15 07:18:02.772] [info] Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
[2025-12-15 07:18:02.772] [info] Checking required entities per dimension
[2025-12-15 07:18:02.772] [info] Cell type: 0 dofmap: 162x10
[2025-12-15 07:18:02.772] [info] Global index computation
[2025-12-15 07:18:02.772] [info] Got 2 index_maps
[2025-12-15 07:18:02.772] [info] Get global indices
[2025-12-15 07:18:02.774] [info] Checking required entities per dimension
[2025-12-15 07:18:02.774] [info] Cell type: 0 dofmap: 162x4
[2025-12-15 07:18:02.775] [info] Global index computation
[2025-12-15 07:18:02.775] [info] Got 1 index_maps
[2025-12-15 07:18:02.775] [info] Get global indices
[2025-12-15 07:18:02.781] [info] Requesting connectivity (2, 0) - (3, 0)

[2025-12-15 07:18:03.777] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:03.777] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:03.777] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:03.777] [info] Requesting connectivity (3, 0) - (2, 0)

[2025-12-15 07:18:09.507] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:09.507] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:09.507] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:09.507] [info] Requesting connectivity (3, 0) - (2, 0)

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'}

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 0x7f1ced4cede0>,), 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...

INFO:scifem.solvers:Newton iteration 1: r (abs) = 22132.799231152836 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 2: r (abs) = 7489.241631882697 (tol=1e-06), r (rel) = 0.33837751626740814 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 3: r (abs) = 259.6378106630803 (tol=1e-06), r (rel) = 0.011730907055698101 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 4: r (abs) = 5.281488194414575 (tol=1e-06), r (rel) = 0.00023862721290945704 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 5: r (abs) = 0.002699330231142725 (tol=1e-06), r (rel) = 1.219606342131051e-07 (tol=1e-10)

INFO:scifem.solvers:Newton iteration 6: r (abs) = 9.37185661658365e-10 (tol=1e-06), r (rel) = 4.234374747949808e-14 (tol=1e-10)

[2025-12-15 07:18:10.021] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:10.021] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-15 07:18:10.024] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-15 07:18:10.028] [info] Column ghost size increased from 0 to 0

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

0.9999999999999979

Here we first compiles form

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

then we assemble the form

dolfinx.fem.assemble_scalar(dolfinx.fem.form(J * ufl.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 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:

Right Cauchy-Green tensor: $\mathbf{C} = \mathbf{F}^T \mathbf{F}$
Left Cauchy-Green tensor: $\mathbf{B} = \mathbf{F} \mathbf{F}^T$
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 we have $\mathbf{u}$ being $\mathbb{P}_2$, i.e second order polynomial but only continous and not continously differentiable that dofs, the graient $\nabla \mathbf{u}$ belongs to a discountinous first order space. Since $\mathbf{E}$ is a function of the $\nabla \mathbf{u}^T \nabla \mathbf{u}$ a reasonable space would be a second order discontinous space

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

[2025-12-15 07:18:10.992] [info] Checking required entities per dimension
[2025-12-15 07:18:10.992] [info] Cell type: 0 dofmap: 162x10
[2025-12-15 07:18:10.992] [info] Global index computation
[2025-12-15 07:18:10.992] [info] Got 1 index_maps
[2025-12-15 07:18:10.992] [info] Get global indices

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 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 [Hol02]

References#

[Hol02]

Gerhard A Holzapfel. Nonlinear solid mechanics: a continuum approach for engineering science. 2002.