Problem 1: Deformation of a Beam#
This example implements Problem 1 from the cardiac mechanics benchmark suite [Land et al. 2015].
Problem Description#
Geometry: A cuboid beam with dimensions \(10 \times 1 \times 1\) mm. The domain is defined as \(\Omega = [0, 10] \times [0, 1] \times [0, 1]\).
Material: Transversely isotropic Guccione material.
Fiber direction \(\mathbf{f}_0 = (1, 0, 0)\) (aligned with the long axis).
Constitutive parameters: \(C = 2.0\) kPa, \(b_f = 8.0\), \(b_t = 2.0\), \(b_{fs} = 4.0\).
The material is incompressible.
Boundary Conditions:
Dirichlet: The face at \(X=0\) is fully clamped (\(\mathbf{u} = \mathbf{0}\)).
Neumann: A pressure load \(P\) is applied to the bottom face (\(Z=0\)). The pressure increases linearly from 0 to 0.004 kPa.
Target Quantity: The deflection (vertical displacement) of the point \((10, 0.5, 1.0)\) at the maximum load. The benchmark reference value is approximately 4.0 - 4.2 mm.
1. Geometry and Mesh#
We create a box mesh of dimensions \(10 \times 1 \times 1\).
L = 10.0
W = 1.0
mesh = dolfinx.mesh.create_box(
MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, W, W]], [30, 3, 3], dolfinx.mesh.CellType.hexahedron,
)
We define markers for the boundary conditions.
left(Marker 1): Face at X=0.bottom(Marker 2): Face at Z=0.
left = 1
bottom = 2
boundaries = [
pulse.Marker(name="left", marker=left, dim=2, locator=lambda x: np.isclose(x[0], 0)),
pulse.Marker(name="bottom", marker=bottom, dim=2, locator=lambda x: np.isclose(x[2], 0)),
]
geo = pulse.Geometry(
mesh=mesh,
boundaries=boundaries,
metadata={"quadrature_degree": 4},
)
2. Constitutive Model#
We use the Guccione model as specified in the benchmark.
Since the problem defines a fixed fiber direction along \(x\), we set \(\mathbf{f}_0, \mathbf{s}_0, \mathbf{n}_0\) to align with the coordinate axes.
material_params = {
"C": pulse.Variable(dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(2.0)), "kPa"),
"bf": dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(8.0)),
"bt": dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(2.0)),
"bfs": dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(4.0)),
}
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)))
n0 = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0.0, 0.0, 1.0)))
material = pulse.Guccione(f0=f0, s0=s0, n0=n0, **material_params)
The problem is purely passive (no active contraction) and incompressible.
model = pulse.CardiacModel(
material=material,
active=active_model,
compressibility=comp_model,
)
3. Boundary Conditions#
Dirichlet BC: Fix all displacement components on the ‘left’ face.
def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
facets = geo.facet_tags.find(left)
dofs = dolfinx.fem.locate_dofs_topological(V, geo.facet_dimension, facets)
u_fixed = dolfinx.fem.Function(V)
u_fixed.x.array[:] = 0.0
return [dolfinx.fem.dirichletbc(u_fixed, dofs)]
Neumann BC: Apply a pressure load to the ‘bottom’ face.
Note: In fenicsx-pulse, a NeumannBC with a scalar traction represents a pressure load \(P\) acting normal to the deformed surface: \(\mathbf{t} = P \mathbf{n}\).
Since the pressure is pushing up against the bottom face (normal pointing down), we define a positive pressure.
pressure = pulse.Variable(dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0)), "kPa")
neumann = pulse.NeumannBC(traction=pressure, marker=bottom)
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann,))
4. Solving the Problem#
We initialize the static problem and solve it incrementally up to the target pressure of 0.004 kPa.
problem = pulse.StaticProblem(model=model, geometry=geo, bcs=bcs)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[13], line 1
----> 1 problem = pulse.StaticProblem(model=model, geometry=geo, bcs=bcs)
File <string>:8, in __init__(self, model, geometry, parameters, bcs, cavities)
File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:85, in StaticProblem.__post_init__(self)
83 self.parameters = parameters
84 self._init_spaces()
---> 85 self._init_forms()
86 logger.debug("Initialized StaticProblem with parameters:")
87 for key, value in self.parameters.items():
File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:518, in StaticProblem._init_forms(self)
515 if self.geometry.markers is None:
516 raise RuntimeError("Missing markers in geometry")
--> 518 R = self.R
519 u = self.states
520 K = self.K(R)
File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:410, in StaticProblem.R(self)
405 @property
406 def R(self):
407 # Order is always (u, cavity pressures, rigid body, p)
409 R = self._empty_form()
--> 410 R_material = self._material_form(self.u, p=self.p)
411 R_cavity = self._cavity_pressure_form(self.u, self.cavity_pressures)
412 R_robin = self._robin_form(self.u)
File /dolfinx-env/lib/python3.12/site-packages/pulse/problem.py:274, in StaticProblem._material_form(self, u, p)
272 forms = self._empty_form()
273 # Integrate over reference configuration dX = J_map * dx
--> 274 forms[0] += ufl.inner(self.model.S(C), 0.5 * var_C) * self.geometry.dx
276 if self.is_incompressible:
277 forms[-1] += (J - 1.0) * self.p_test * self.geometry.dx
File /dolfinx-env/lib/python3.12/site-packages/pulse/cardiac_model.py:111, in CardiacModel.S(self, C, C_dot)
104 def S(
105 self,
106 C: ufl.core.expr.Expr,
107 C_dot: ufl.core.expr.Expr | None = None,
108 ) -> ufl.core.expr.Expr:
109 """Cauchy stress for the cardiac model."""
--> 111 S = self.material.S(C, dev=True) + self.active.S(C, dev=True) + self.compressibility.S(C)
112 if C_dot is not None:
113 S += self.viscoelasticity.S(C_dot)
File /dolfinx-env/lib/python3.12/site-packages/pulse/active_model.py:133, in ActiveModel.S(self, C, dev)
131 else:
132 Cdev = C
--> 133 return 2.0 * ufl.diff(self.strain_energy(Cdev), Cdev)
File /dolfinx-env/lib/python3.12/site-packages/ufl/operators.py:382, in diff(f, v)
380 return VariableDerivative(f, v)
381 else:
--> 382 raise ValueError("Expecting a Variable or SpatialCoordinate in diff.")
ValueError: Expecting a Variable or SpatialCoordinate in diff.
target_pressure = 0.004
steps = [0.0005, 0.001, 0.002, 0.003, 0.004]
for p in steps:
print(f"Solving for pressure = {p} kPa")
pressure.assign(p)
problem.solve()
5. Post-processing#
We save the final displacement and evaluate the deflection at the specific target point \((10, 0.5, 1)\).
with dolfinx.io.VTXWriter(mesh.comm, "problem1_displacement.bp", [problem.u], engine="BP4") as vtx:
vtx.write(0.0)
# Evaluate displacement at the target point
point = np.array([10.0, 0.5, 1.0])
# Use bounding box tree for point collision
tree = dolfinx.geometry.bb_tree(mesh, mesh.topology.dim)
cell_candidates = dolfinx.geometry.compute_collisions_points(tree, point)
cell = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates, point)
# Get the z-displacement
if len(cell.array) > 0:
# Evaluate only if the point is found on this process
u_val = problem.u.eval(point, cell.array)[2]
else:
u_val = -np.inf # Placeholder for reduction
# Reduce across all processes (take the max to get the actual value)
uz = mesh.comm.allreduce(u_val, op=MPI.MAX)
result = point[2] + uz
print(f"Target Point: {point}")
print(f"Vertical Deflection (Uz): {uz:.4f} mm")
print(f"Final Z Position: {result:.4f} mm")
# Verify against benchmark tolerance (approx 4.0 - 4.2 mm deflection)
# Note: Result is typically around 4.17 mm for the settings described.
if mesh.comm.rank == 0:
print(f"Benchmark Ref: ~4.17 mm")
# Visualization with PyVista
try:
import pyvista
except ImportError:
print("Pyvista is not installed")
else:
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
uh = dolfinx.fem.Function(V)
uh.interpolate(problem.u)
p = pyvista.Plotter()
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
p.add_mesh(grid, style="wireframe", color="k", opacity=0.5, label="Reference")
warped = grid.warp_by_vector("u", factor=1.0)
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("deflection.png")