Hyperelasticity#
Author: Jørgen S. Dokken and Garth N. Wells
This section shows how to solve the hyperelasticity problem for deformation of a beam.
We will also show how to create a constant boundary condition for a vector function space.
We start by importing DOLFINx and some additional dependencies.
Then, we create a slender cantilever consisting of hexahedral elements and create the function space V for our unknown.
from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
import pyvista
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh, plot
L = 20.0
domain = mesh.create_box(
MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron
)
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim,)))
We create two python functions for determining the facets to apply boundary conditions to
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], L)
fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
Next, we create a marker based on these two functions
Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets]
)
We then create a function for supplying the boundary condition on the left side, which is fixed.
u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
To apply the boundary condition, we identity the dofs located on the facets marked by the MeshTag.
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]
Next, we define the body force on the reference configuration (B), and nominal (first Piola-Kirchhoff) traction (T).
B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
Define the test and solution functions on the space \(V\)
v = ufl.TestFunction(V)
u = fem.Function(V)
Define kinematic quantities used in the problem
# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
Define the elasticity model via a stored strain energy density function \(\psi\), and create the expression for the first Piola-Kirchhoff stress:
Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
Stored strain energy density (compressible neo-Hookean model)
Hyper-elasticity
P = ufl.diff(psi, F)
Comparison to linear elasticity
To illustrate the difference between linear and hyperelasticity, the following lines can be uncommented to solve the linear elasticity problem.
# P = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
Define the variational form with traction integral over all facets with value 2. We set the quadrature degree for the integrals to 4.
metadata = {"quadrature_degree": 4}
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
Define the residual of the equation (we want to find u such that residual(u) = 0)
As the varitional form is non-linear and written on residual form, we use the non-linear problem class from DOLFINx to set up required structures to use a Newton solver.
petsc_options = {
"snes_type": "newtonls",
"snes_linesearch_type": "none",
"snes_monitor": None,
"snes_atol": 1e-8,
"snes_rtol": 1e-8,
"snes_stol": 1e-8,
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
}
problem = NonlinearProblem(
residual,
u,
bcs=bcs,
petsc_options=petsc_options,
petsc_options_prefix="hyperelasticity",
)
We create a function to plot the solution at each time step.
Each process builds the grid of the cells it owns, the pieces are gathered on one process and merged into a single grid, which is then drawn, see Plotting in parallel.
The mesh does not change during the simulation, so the pieces are gathered once, and only the values are sent to the process that writes the GIF at each frame.
topology, cells, geometry = plot.vtk_mesh(u.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, : len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
function_grid["u"] = values
function_grid.set_active_vectors("u")
# Compute magnitude of displacement to visualize in GIF
Vs = fem.functionspace(domain, ("Lagrange", 2))
magnitude = fem.Function(Vs)
us = fem.Expression(
ufl.sqrt(sum([u[i] ** 2 for i in range(len(u))])), Vs.element.interpolation_points
)
magnitude.interpolate(us)
function_grid["mag"] = magnitude.x.array
root = 0
assert root < domain.comm.size, f"Cannot gather on process {root} of {domain.comm.size}"
pieces = domain.comm.gather(function_grid, root=root)
plotter = None
if pieces is not None:
plotter = pyvista.Plotter()
plotter.open_gif("deformation.gif", fps=3)
def write_frame(displacements: list[np.ndarray], magnitudes: list[np.ndarray]):
"""Update the values of each piece, and draw the merged grid as a single frame."""
assert plotter is not None and pieces is not None
for piece, displacement, mag in zip(pieces, displacements, magnitudes):
piece["u"][:, : displacement.shape[1]] = displacement
piece["mag"][:] = mag
warped = pyvista.merge(pieces).warp_by_vector("u", factor=1)
warped.set_active_scalars("mag")
plotter.clear()
plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10])
plotter.update_scalar_bar_range([0, 10])
plotter.write_frame()
2026-10-19 15:06:05.576 ( 1.887s) [ 7FE00BC85140]vtkXOpenGLRenderWindow.:1460 WARN| bad X server connection. DISPLAY=
Finally, we solve the problem over several time steps, updating the z-component of the traction
log.set_log_level(log.LogLevel.INFO)
tval0 = -1.5
for n in range(1, 10):
T.value[2] = n * tval0
problem.solve()
converged = problem.solver.getConvergedReason()
num_its = problem.solver.getIterationNumber()
assert converged > 0, f"Solver did not converge with reason {converged}."
print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")
magnitude.interpolate(us)
# The gathers are collective, so every process takes part in them, while only the process
# that received the pieces draws the frame
local_u = u.x.array.reshape(geometry.shape[0], len(u))
displacements = domain.comm.gather(local_u.copy(), root=root)
magnitudes = domain.comm.gather(magnitude.x.array.copy(), root=root)
if displacements is not None and magnitudes is not None:
write_frame(displacements, magnitudes)
if plotter is not None:
plotter.close()
0 SNES Function norm 1.633333333333e-01
1 SNES Function norm 1.330753372046e+03
2 SNES Function norm 7.390366887036e+01
3 SNES Function norm 4.714331050662e-01
4 SNES Function norm 9.045604307033e-01
5 SNES Function norm 1.011600114129e-03
6 SNES Function norm 2.328091526068e-05
7 SNES Function norm 1.801158348430e-10
Time step 1, Number of iterations 7, Load [ 0. 0. -1.5]
0 SNES Function norm 1.633333333294e-01
1 SNES Function norm 1.053155719429e+03
2 SNES Function norm 4.497738445445e+01
3 SNES Function norm 1.597191951891e+00
4 SNES Function norm 1.008829277909e+01
5 SNES Function norm 3.363230749376e-02
6 SNES Function norm 1.161567348044e-01
7 SNES Function norm 4.338363890860e-06
8 SNES Function norm 2.488038214109e-09
Time step 2, Number of iterations 8, Load [ 0. 0. -3.]
0 SNES Function norm 1.633333331821e-01
1 SNES Function norm 6.157586973378e+02
2 SNES Function norm 1.460711281490e+01
3 SNES Function norm 3.318947791600e+00
4 SNES Function norm 1.009685403497e+01
5 SNES Function norm 4.185327730234e-01
6 SNES Function norm 1.058465191042e+00
7 SNES Function norm 3.889936028319e-03
8 SNES Function norm 1.252948658089e-04
9 SNES Function norm 3.763962823680e-10
Time step 3, Number of iterations 9, Load [ 0. 0. -4.5]
0 SNES Function norm 1.633333333282e-01
1 SNES Function norm 3.351285834259e+02
2 SNES Function norm 4.135214428429e+00
3 SNES Function norm 4.493073543128e+01
4 SNES Function norm 3.405823793819e-01
5 SNES Function norm 1.363506117312e+01
6 SNES Function norm 7.670312200733e-03
7 SNES Function norm 5.774813103132e-03
8 SNES Function norm 7.146043024778e-09
Time step 4, Number of iterations 8, Load [ 0. 0. -6.]
0 SNES Function norm 1.633333325082e-01
1 SNES Function norm 1.880704733240e+02
2 SNES Function norm 1.252286182893e+00
3 SNES Function norm 3.609070569483e+00
4 SNES Function norm 4.953401642992e-02
5 SNES Function norm 9.623194816626e-03
6 SNES Function norm 3.341941014384e-07
7 SNES Function norm 4.215321827480e-10
Time step 5, Number of iterations 7, Load [ 0. 0. -7.5]
0 SNES Function norm 1.633333333370e-01
1 SNES Function norm 1.116094096259e+02
2 SNES Function norm 4.247396796581e-01
3 SNES Function norm 1.262748205164e+00
4 SNES Function norm 2.317735947565e-03
5 SNES Function norm 4.305376600761e-05
6 SNES Function norm 4.326984373333e-10
Time step 6, Number of iterations 6, Load [ 0. 0. -9.]
0 SNES Function norm 1.633333333329e-01
1 SNES Function norm 7.005111549876e+01
2 SNES Function norm 1.606542965747e-01
3 SNES Function norm 5.406329906822e-01
4 SNES Function norm 1.096042492949e-04
5 SNES Function norm 3.000172474539e-07
6 SNES Function norm 4.580860336378e-10
Time step 7, Number of iterations 6, Load [ 0. 0. -10.5]
0 SNES Function norm 1.633333333307e-01
1 SNES Function norm 4.622497138501e+01
2 SNES Function norm 6.753406519374e-02
3 SNES Function norm 2.415083558299e-01
4 SNES Function norm 4.656928264076e-06
5 SNES Function norm 2.074258950994e-09
Time step 8, Number of iterations 5, Load [ 0. 0. -12.]
0 SNES Function norm 1.633333331308e-01
1 SNES Function norm 3.185930868983e+01
2 SNES Function norm 3.414437561952e-02
3 SNES Function norm 1.076837932702e-01
4 SNES Function norm 4.272869537342e-07
5 SNES Function norm 4.934057166211e-10
Time step 9, Number of iterations 5, Load [ 0. 0. -13.5]
