From XML data
From XML dataΒΆ
Example on how to load your data if it is currently in xml format
import json
import dolfin
from fenics_plotly import plot
import pulse
try:
from dolfin_adjoint import Constant, DirichletBC, Function, Mesh, interpolate
except ImportError:
from dolfin import Mesh, Function, Constant, DirichletBC, interpolate
# Mesh
mesh = Mesh(pulse.utils.mpi_comm_world(), "data/mesh.xml")
# Marker functions
facet_function = dolfin.MeshFunction("size_t", mesh, "data/facet_function.xml")
marker_functions = pulse.MarkerFunctions(ffun=facet_function)
# Markers
with open("data/markers.json", "r") as f:
markers = json.load(f)
# Fiber
fiber_element = dolfin.VectorElement(
family="Quadrature",
cell=mesh.ufl_cell(),
degree=4,
quad_scheme="default",
)
fiber_space = dolfin.FunctionSpace(mesh, fiber_element)
fiber = Function(fiber_space, "data/fiber.xml")
microstructure = pulse.Microstructure(f0=fiber)
# Create the geometry
geometry = pulse.HeartGeometry(
mesh=mesh,
markers=markers,
marker_functions=marker_functions,
microstructure=microstructure,
)
activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
activation.assign(Constant(0.0))
matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
activation=activation,
parameters=matparams,
f0=geometry.f0,
)
# LV Pressure
lvp = Constant(0.0)
lv_marker = markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]
# Add spring term at the base with stiffness 1.0 kPa/cm^2
base_spring = 1.0
robin_bc = [pulse.RobinBC(value=Constant(base_spring), marker=markers["BASE"][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, markers["BASE"][0])
return bc
dirichlet_bc = [fix_basal_plane]
# You can also use a built in function for this
# from functools import partial
# dirichlet_bc = partial(pulse.mechanicsproblem.dirichlet_fix_base_directional,
# ffun=geometry.ffun,
# marker=geometry.markers["BASE"][0])
# Collect boundary conditions
bcs = pulse.BoundaryConditions(
dirichlet=dirichlet_bc,
neumann=neumann_bc,
robin=robin_bc,
)
# Create the problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
problem.solve()
(2, True)
# Solve the problem
pulse.iterate.iterate(problem, (lvp, activation), (1.0, 0.2))
2022-06-05 11:38:35,760 [565] INFO pulse.iterate: Iterating to target control (f_26)...
2022-06-05 11:38:35,763 [565] INFO pulse.iterate: Current control: f_26 = 0.000,0.000
2022-06-05 11:38:35,764 [565] INFO pulse.iterate: Target: 1.000,0.200
*** Warning: PETSc SNES solver diverged in 1 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.
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 76),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 299),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 366),
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))), 484),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 535),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 586),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 645),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 728),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 796)],
[(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 72),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 74)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 295),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 297)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 362),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 364)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 421),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 423)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 480),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 482)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 531),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 533)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 582),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 584)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 641),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 643)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 724),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 726)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 792),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 794))])
# Get the solution
u, p = problem.state.split(deepcopy=True)
volume = geometry.cavity_volume(u=u)
print(f"Cavity volume = {volume}")
Cavity volume = 1.118996302419786
# 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, color="red", show=False)
fig.add_plot(plot(mesh, show=False))
fig.show()