Simple ellipsoid

import dolfin
from fenics_plotly import plot

import pulse

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

gamma_space = "R_0"

geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["simple_ellipsoid"])
# geometry = pulse.geometries.prolate_ellipsoid_geometry(mesh_size_factor=3.0)

2024-09-05 12:04:24,062 [1147] INFO     pulse.geometry_utils: 
Load mesh from h5

if gamma_space == "regional":
    activation = pulse.RegionalParameter(geometry.cfun)
    target_activation = pulse.dolfin_utils.get_constant(0.2, len(activation))
else:
    activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
    target_activation = Constant(0.2)

matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
    activation=activation,
    parameters=matparams,
    f0=geometry.f0,
    s0=geometry.s0,
    n0=geometry.n0,
)

# LV Pressure
lvp = Constant(0.0)
lv_marker = geometry.markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]

# Add spring term at the epicardium of stiffness 1.0 kPa/cm^2 to represent pericardium
base_spring = 1.0
robin_bc = [
    pulse.RobinBC(value=Constant(base_spring), marker=geometry.markers["EPI"][0]),
]

# Fix the basal plane in the longitudinal direction
# 0 in V.sub(0) refers to x-direction, which is the longitudinal direction
def fix_basal_plane(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    bc = DirichletBC(
        V.sub(0),
        Constant(0.0),
        geometry.ffun,
        geometry.markers["BASE"][0],
    )
    return bc

dirichlet_bc = (fix_basal_plane,)

# Collect boundary conditions
bcs = pulse.BoundaryConditions(
    dirichlet=dirichlet_bc,
    neumann=neumann_bc,
    robin=robin_bc,
)

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

# Solve the problem
pulse.iterate.iterate(problem, (lvp, activation), (1.0, target_activation))

2024-09-05 12:04:24,210 [1147] INFO     pulse.iterate: Iterating to target control (f_36)...

2024-09-05 12:04:24,211 [1147] INFO     pulse.iterate: Current control: f_36 = 0.000,0.000

2024-09-05 12:04:24,211 [1147] INFO     pulse.iterate: Target: 1.000,0.200

2024-09-05 12:04:24,336 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,336 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,337 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,337 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,338 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,339 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,340 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,340 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,341 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,341 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:24,342 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:04:25,337 [1147] DEBUG    UFL_LEGACY: Blocks of each mode: 
  99	full

2024-09-05 12:04:25,558 [1147] DEBUG    UFL_LEGACY: Blocks of each mode: 
  18	full

2024-09-05 12:04:25,632 [1147] DEBUG    UFL_LEGACY: Blocks of each mode: 

2024-09-05 12:05:29,640 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,641 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,641 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,642 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,642 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,643 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,643 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,644 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,644 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,645 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,645 [1147] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2024-09-05 12:05:29,819 [1147] DEBUG    UFL_LEGACY: Blocks of each mode: 
  10	full

2024-09-05 12:05:30,008 [1147] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full

2024-09-05 12:05:30,081 [1147] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full

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

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

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

*** Warning: PETSc SNES solver diverged in 3 iterations with divergence reason DIVERGED_DTOL.

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

*** Warning: PETSc SNES solver diverged in 3 iterations with divergence reason DIVERGED_DTOL.

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

*** 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))), 57),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 250),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 315),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 372),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 429),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 478),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 527),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 606),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 655),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 728),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 831),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 918),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1005),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1070),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1112)],
 [(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 53),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 55)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 246),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 248)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 311),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 313)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 368),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 370)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 425),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 427)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 474),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 476)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 523),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 525)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 602),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 604)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 651),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 653)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 724),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 726)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 827),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 829)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 914),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 916)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1001),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1003)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1066),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1068)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1108),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1110))])

# Get the solution
u, p = problem.state.split(deepcopy=True)

# Move mesh accoring to displacement
u_int = interpolate(u, dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1))
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)

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