Cardiac mechanics Benchmark - Problem 3¶

This code implements problem 3 in the Cardiac Mechanic Benchmark paper

Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.

import dolfin

try:
    from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
    from dolfin import DirichletBC, Constant, interpolate, Mesh

import pulse
from fenics_plotly import plot

geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["benchmark"])
# geometry = pulse.geometries.benchmark_ellipsoid_geometry()

2022-06-05 11:55:55,404 [908] INFO     pulse.geometry_utils: 
Load mesh from h5

# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["CC"] = 2.0
material_parameters["bf"] = 8.0
material_parameters["bfs"] = 4.0
material_parameters["bt"] = 2.0

activation = Constant(0.0)
material = pulse.Guccione(
    parameters=material_parameters,
    active_model=pulse.ActiveModels.active_stress,
    activation=activation,
)

# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    return DirichletBC(
        V,
        Constant((0.0, 0.0, 0.0)),
        geometry.ffun,
        geometry.markers["BASE"][0],
    )

# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])

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

# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)

# Solve problem
pulse.iterate.iterate(problem, lvp, 15.0, initial_number_of_steps=50)
pulse.iterate.iterate(problem, activation, 60.0, initial_number_of_steps=50)

2022-06-05 11:55:55,653 [908] INFO     pulse.iterate: Iterating to target control (f_28)...

2022-06-05 11:55:55,655 [908] INFO     pulse.iterate: Current control: f_28 = 0.000

2022-06-05 11:55:55,656 [908] INFO     pulse.iterate: Target: 15.000

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

2022-06-05 11:59:21,692 [908] INFO     pulse.iterate: Iterating to target control (f_20)...

2022-06-05 11:59:21,693 [908] INFO     pulse.iterate: Current control: f_20 = 0.000

2022-06-05 11:59:21,694 [908] INFO     pulse.iterate: Target: 60.000

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 496),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 529),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 554),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 587),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 620),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 653),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 694),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 743),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 813),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 862)],
 [Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 494),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 527),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 552),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 585),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 618),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 651),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 692),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 741),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 811),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 860)])

# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)

endo_apex_marker = geometry.markers["APEX_ENDO"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)

print(
    ("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
        endo_apex_pos[0],
    ),
)

Get longitudinal position of endocardial apex: 19.3891 mm

epi_apex_marker = geometry.markers["APEX_EPI"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)

print(
    ("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
        epi_apex_pos[0],
    ),
)

Get longitudinal position of epicardial apex: 20.7619 mm

V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)

fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()