Nonlinear elasticity with blocked Newton solver

Nonlinear elasticity with blocked Newton solver#

Author: Henrik N. T. Finsberg

SPDX-License-Identifier: MIT

In this example we will solve a nonlinear elasticity problem with a non-linear Lagrange multiplier.

We consider a unit cube domain \(\Omega = [0, 1]^3\) with Dirichlet boundary conditions on the left face and traction force on the right face, and we seek a displacement field \(\mathbf{u}: \Omega \to \mathbb{R}^3\) that solves the momentum balance equation

\[\begin{split} \begin{align} \nabla \cdot \mathbf{P} = 0, \quad \mathbf{X} \in \Omega, \\ \end{align} \end{split}\]

where \(\mathbf{P}\) is the first Piola-Kirchhoff stress tensor, and \(\mathbf{X}\) is the reference configuration. We consider a Neo-Hookean material model, where the strain energy density is given by

\[ \begin{align} \psi = \frac{\mu}{2}(\text{tr}(\mathbf{C}) - 3), \end{align} \]

and the first Piola-Kirchhoff stress tensor is given by

\[ \begin{align} \mathbf{P} = \frac{\partial \psi}{\partial \mathbf{F}} \end{align} \]

where \(\mathbf{F} = \nabla \mathbf{u} + \mathbf{I}\) is the deformation gradient, \(\mathbf{C} = \mathbf{F}^T \mathbf{F}\) is the right Cauchy-Green tensor, \(\mu\) is the shear modulus, and \(p\) is the pressure. We also enforce the incompressibility constraint

\[ \begin{align} J = \det(\mathbf{F}) = 1, \end{align} \]

so that the total Lagrangian is given by

\[ \begin{align} \mathcal{L}(\mathbf{u}, p) = \int_{\Omega} \psi \, dx - \int_{\partial \Omega} t \cdot \mathbf{u} \, ds + \int_{\Omega} p(J - 1) \, dx. \end{align} \]

Here \(t\) is the traction force which is set to \(10\) on the right face of the cube and \(0\) elsewhere. The Euler-Lagrange equations for this problem are given by: Find \(\mathbf{u} \in V\) and \(p \in Q\) such that

\[\begin{split} \begin{align} D_{\delta \mathbf{u} } \mathcal{L}(\mathbf{u}, p) = 0, \quad \forall \delta \mathbf{u} \in V, \\ D_{\delta p} \mathcal{L}(\mathbf{u}, p) = 0, \quad \forall \delta p \in Q, \end{align} \end{split}\]

where \(V\) is the displacement space and \(Q\) is the pressure space. For this we select \(Q_2/P_1\) elements i.e second order Lagrange elements for \(\mathbf{u}\) and first order discontinuous polynomial cubical elements for \(p\), which is a stable element for incompressible elasticity [ABeiraodVL+13]. Note also that the Euler-Lagrange equations can be derived automatically using ufl.

from pathlib import Path
from mpi4py import MPI
import numpy as np
import ufl
import dolfinx
import dolfinx.fem.petsc
import scifem

We create the mesh and the function spaces

mesh = dolfinx.mesh.create_unit_cube(
    MPI.COMM_WORLD, 3, 3, 3, dolfinx.mesh.CellType.hexahedron, dtype=np.float64
)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2, (3,)))
Q = dolfinx.fem.functionspace(mesh, ("DPC", 1,))
W = ufl.MixedFunctionSpace(V, Q)

And the test and trial functions

v, q = ufl.TestFunctions(W)
du, dp = ufl.TrialFunctions(W)
u = dolfinx.fem.Function(V)
p = dolfinx.fem.Function(Q)

Next we create the facet tags for the left and right faces

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], 1)

facet_tags = scifem.create_entity_markers(
    mesh, mesh.topology.dim - 1, [(1, left, True), (2, right, True)]
)

We create the Dirichlet boundary conditions on the left face

facets_left = facet_tags.find(1)
dofs_left = dolfinx.fem.locate_dofs_topological(V, 2, facets_left)
u_bc_left = dolfinx.fem.Function(V)
u_bc_left.x.array[:] = 0
bc = dolfinx.fem.dirichletbc(u_bc_left, dofs_left)

Define the deformation gradient, right Cauchy-Green tensor, and invariants of the deformation tensors

d = len(u)
I = ufl.Identity(d)             # Identity tensor
F = I + ufl.grad(u)             # Deformation gradient
C = F.T*F                       # Right Cauchy-Green tensor
I1 = ufl.tr(C)                  # First invariant of C
J  = ufl.det(F)                 # Jacobian of F

Traction for to be applied on the right face

t = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(10.0))

N = ufl.FacetNormal(mesh)

# Material parameters and strain energy density
mu = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(10.0))
psi = (mu / 2)*(I1 - 3)

We for the total Lagrangian

L = psi*ufl.dx - ufl.inner(u, t * N)*ufl.ds(subdomain_data=facet_tags, subdomain_id=2)  + p * (J - 1) * ufl.dx

and take the first variation of the total Lagrangian to obtain the residual

R = ufl.derivative(L, [u, p], [v, q])
if np.issubdtype(dolfinx.default_scalar_type, np.complexfloating):
    # The derivative direction must be an Argument, so conjugate the test functions afterwards
    R = ufl.replace(R, {v: ufl.conj(v), q: ufl.conj(q)})
R_blocked = ufl.extract_blocks(R)

We do the same for the second variation to obtain the Jacobian

K = ufl.derivative(R, [u, p], [du, dp])
K_blocked = ufl.extract_blocks(K)

Now we can create the Newton solver and solve the problem

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True,
                 "snes_error_if_not_converged": True, "snes_monitor": None}
problem = dolfinx.fem.petsc.NonlinearProblem(R_blocked, [u, p], bcs=[bc], J=K_blocked, petsc_options=petsc_options,
                                             petsc_options_prefix="elasticity_")

The underlying PETSc SNES solver is accessible through problem.solver, which we can inspect before and after the solve

def pre_solve(problem: dolfinx.fem.petsc.NonlinearProblem):
    max_iterations = problem.solver.getTolerances()[3]
    print(f"Starting solve with at most {max_iterations} iterations")

def post_solve(problem: dolfinx.fem.petsc.NonlinearProblem):
    num_iterations = problem.solver.getIterationNumber()
    dx_norm = problem.solver.getSolutionUpdate().norm()
    print(f"Solve completed in {num_iterations} iterations with correction norm {dx_norm}")

pre_solve(problem)
problem.solve()
post_solve(problem)

Starting solve with at most 50 iterations
SNES Function norm 3.549248521349e+00
SNES Function norm 8.491655649000e-01
SNES Function norm 1.046210159005e-01
SNES Function norm 8.058185269344e-04
SNES Function norm 5.106834087519e-08
SNES Function norm 7.304484941955e-15
Solve completed in 5 iterations with correction norm 1.4961407110982043e-06

Finally, we can visualize the solution using pyvista We first generate the grids locally on each MPI rank, before we gather the results on rank root for visualization.

import pyvista
pyvista.global_theme.allow_empty_mesh= True
root = 0
tdim = mesh.topology.dim
cell_map = mesh.topology.index_map(tdim)
num_owned_cells = cell_map.size_local
owned_cells = np.arange(num_owned_cells, dtype=np.int32)
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V, owned_cells)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
linear_topology, linear_cell_types, linear_geometry = dolfinx.plot.vtk_mesh(mesh, tdim, owned_cells)

if mesh.geometry.cmaps[0].degree == 1:
    _dx_to_pv_lin = {dolfinx.mesh.CellType.hexahedron: pyvista.CellType.HEXAHEDRON,
                     dolfinx.mesh.CellType.tetrahedron: pyvista.CellType.TETRA,
                     dolfinx.mesh.CellType.triangle: pyvista.CellType.TRIANGLE,
                     dolfinx.mesh.CellType.quadrilateral: pyvista.CellType.QUAD}
    pv_ct = _dx_to_pv_lin[mesh.topology.cell_type]
    linear_cell_types = np.full_like(linear_cell_types, pv_ct)
    linear_grid = pyvista.UnstructuredGrid(linear_topology, linear_cell_types, linear_geometry)
grid["u"] = u.x.array.real.reshape((geometry.shape[0], 3))
gathered_grid = mesh.comm.gather(grid, root=root)
gathered_linear_grid = mesh.comm.gather(linear_grid, root=root)
if gathered_grid is not None:
    p = pyvista.Plotter()
    merged_grid = pyvista.merge(gathered_grid)
    merged_linear_grid = pyvista.merge(gathered_linear_grid)
    p.add_mesh(merged_linear_grid, style="wireframe", color="k")
    warped = merged_grid.warp_by_vector("u", factor=1.5)
    p.add_mesh(warped, show_edges=False)
    p.show_axes()
    figure = Path("blocked_solver.png")
    p.export_vtksz(figure.with_suffix(".vtksz"))
    p.show(screenshot=figure)

/tmp/ipykernel_487/1404549386.py:21: PyVistaDeprecationWarning: Plotter.export_vtksz is deprecated. Use plotter.trame.export_vtksz(...) instead. Install via `pip install trame-pyvista`.
  p.export_vtksz(figure.with_suffix(".vtksz"))

References#

[ABeiraodVL+13]

Ferdinando Auricchio, Lourenço Beirão da Veiga, Carlo Lovadina, Alessandro Reali, Robert L Taylor, and Peter Wriggers. Approximation of incompressible large deformation elastic problems: some unresolved issues. Computational Mechanics, 52:1153–1167, 2013.

Nonlinear elasticity with blocked Newton solver

Contents

Nonlinear elasticity with blocked Newton solver#

References#