<no title>

# # 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.
#
from pathlib import Path
import dolfin

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

import pulse
import cardiac_geometries as cg
from fenics_plotly import plot

geo_path = Path("geometry")
if not geo_path.is_dir():
    cg.create_benchmark_geometry_land15(outdir=geo_path)
geo = cg.geometry.Geometry.from_folder(geo_path)
geometry = pulse.HeartGeometry(
    mesh=geo.mesh,
    markers=geo.markers,
    marker_functions=pulse.MarkerFunctions(vfun=geo.vfun, ffun=geo.ffun),
    microstructure=pulse.Microstructure(f0=geo.f0, s0=geo.s0, n0=geo.n0),
)

# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["C"] = 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)

2024-09-05 11:59:02,228 [941] INFO     pulse.iterate: Iterating to target control (f_51)...

2024-09-05 11:59:02,229 [941] INFO     pulse.iterate: Current control: f_51 = 0.000

2024-09-05 11:59:02,230 [941] INFO     pulse.iterate: Target: 15.000

2024-09-05 11:59:02,346 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,347 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,348 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,348 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,349 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,350 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,351 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,351 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:02,352 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 11:59:03,336 [941] DEBUG    UFL_LEGACY: Blocks of each mode: 
  99	full

2024-09-05 11:59:03,545 [941] DEBUG    UFL_LEGACY: Blocks of each mode: 
  18	full

2024-09-05 12:00:15,154 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,155 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,155 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,156 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,156 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,157 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,157 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,158 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,158 [941] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:00:15,339 [941] DEBUG    UFL_LEGACY: Blocks of each mode: 
  10	full

2024-09-05 12:00:15,531 [941] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full

2024-09-05 12:00:37,536 [941] INFO     pulse.iterate: Iterating to target control (f_44)...

2024-09-05 12:00:37,536 [941] INFO     pulse.iterate: Current control: f_44 = 0.000

2024-09-05 12:00:37,537 [941] 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), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 425),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 458),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 491),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 524),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 557),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 598),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 639),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 688),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 774),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 839)],
 [Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 423),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 456),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 489),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 522),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 555),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 596),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 637),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 686),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 772),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 837)])

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

endo_apex_marker = geometry.markers["ENDOPT"][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: -11.9845 mm

epi_apex_marker = geometry.markers["EPIPT"][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: -15.2197 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()

Contents