Linear elasticity¶
This work is adapted from the FEniCS tutorial: https://fenicsproject.org/pub/tutorial/html/._ftut1008.html#ftut:elast
Having studied a few scalar-valued problems, we now move on to a vector-valued problem. The equations of linear elasticity. Here, we'll treat the isotropic case.
For small deformations, the governing equation is:
$$ -\nabla \cdot \sigma = f \text{ in } \Omega, $$ with $$ \DeclareMathOperator{\Tr}{Tr} \text{the stress tensor}\quad \sigma := \lambda \Tr(\epsilon)\mathbb{I} + 2\mu\epsilon\\ \text{and the symmetric strain tensor}\quad \epsilon := \frac{1}{2}\left(\nabla u + (\nabla u)^T\right), $$ where $u$ is the unknown vector displacement field, and $\mu$ and $\lambda$ are the Lamè parameters.
As before, the variational formulation consists of multiplying by a test function in some suitable finite element space, $v \in V$, and integrating. Note that this time, the solution $u$, and hence the test space $V$ are vector-valued (so multiplication actually means taking the inner product).
We obtain
$$ -\int_\Omega (\nabla \cdot \sigma)\cdot v\,\mathrm{d}x = \int_\Omega f \cdot v\,\mathrm{d}x. $$
Since $\sigma$ is actually a function of derivatives of $u$, we must integrate this term by parts, resulting in
$$ \int_\Omega \sigma : \nabla v\,\mathrm{d} x - \int_\Gamma (\sigma \cdot n)\cdot v\,\mathrm{d} s. = \int_\Omega f \cdot v\,\mathrm{d}x.$$
We also need to specify boundary conditions. We can do so either by prescribing the displacement $u$ on the boundary, or the traction $\sigma \cdot n$. The former is a strong or Dirichlet condition, the latter a weak or Neumann condition.
Let us decide on a concrete setting. We will solve for the displacement of a beam under its own weight clamped at one end. That is, we will take $\Omega = [0, L] \times [0, W]$, we set $u = (0, 0)$ on $\Gamma_D$, the plane $x = 0$. On all other boundaries, we have traction-free conditions.
We start, as usual, by import Firedrake and defining a mesh
%matplotlib inline
import matplotlib.pyplot as plt
# Load Firedrake on Colab
try:
import firedrake
except ImportError:
!wget "https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh" -O "/tmp/firedrake-install.sh" && bash "/tmp/firedrake-install.sh"
import firedrake
from gadopt import *
from firedrake.pyplot import triplot
length = 1
width = 0.2
mesh = RectangleMesh(40, 20, length, width)
We need a function space for the solution variable $u$, this space is vector valued, so we use a VectorFunctionSpace
. By default, this constructs a space where the vectors have as many components as the geometric dimension of the mesh (two in this case).
V = VectorFunctionSpace(mesh, "Lagrange", 1)
We need a boundary condition object for $\Gamma_D$.
bc = DirichletBC(V, Constant([0, 0]), 1)
Now let's define the material parameters. The deformation due to gravity can be obtained by setting the load vector $f = (0, -\rho g)$ where $\rho$ is the material density and $g$ the acceleration due to gravity. We'll choose $\rho = 0.01$, and $g = 1$.
rho = Constant(0.01)
g = Constant(1)
f = as_vector([0, -rho*g])
mu = Constant(1)
lambda_ = Constant(0.25)
Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor
Now we'll define functions that construct the symbolic expressions for the stress and strain.
def epsilon(u):
return 0.5*(grad(u) + grad(u).T)
def sigma(u):
return lambda_*div(u)*Id + 2*mu*epsilon(u)
The variational problem can now be solved. Passing the solver parameter "ksp_monitor": None
tells PETSc to print the progress of the linear solver to screen. Firedrake uses a direct solver by default, so it should converge in one iteration.
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx
uh = Function(V)
solve(a == L, uh, bcs=bc, solver_parameters={"ksp_monitor": None})
Residual norms for firedrake_0_ solve. 0 KSP Residual norm 6.916415658096e-05 1 KSP Residual norm 5.934600681034e-15
Visualising the solution¶
We solved the equations in displacement formulation. That is, the $u_h$ we obtain is a perturbation to the original coordinate field of the mesh. If we want to view the output, we can either use the original mesh, or we can create a new mesh with the displaced coordinates: $\hat{X} = X + u_h$.
displaced_coordinates = interpolate(SpatialCoordinate(mesh) + uh, V)
displaced_mesh = Mesh(displaced_coordinates)
/home/firedrake/firedrake/src/firedrake/firedrake/interpolation.py:385: FutureWarning: The use of `interpolate` to perform the numerical interpolation is deprecated. This feature will be removed very shortly. Instead, import `interpolate` from the `firedrake.__future__` module to update the interpolation's behaviour to return the symbolic `ufl.Interpolate` object associated with this interpolation. You can then assemble the resulting object to get the interpolated quantity of interest. For example, ``` from firedrake.__future__ import interpolate ... assemble(interpolate(expr, V)) ``` Alternatively, you can also perform other symbolic operations on the interpolation operator, such as taking the derivative, and then assemble the resulting form. warnings.warn("""The use of `interpolate` to perform the numerical interpolation is deprecated.
This created a new mesh with a coordinate field provided by displaced_coordinates
that nevertheless shares a topology with the original, regular, mesh. We could, if we wanted to, go ahead and solve variational problems on this new mesh, however, we've only done this so we can see the elastic deformation.
Here we're using the matplotlib API for Axes objects to make the horizontal and vertical axes equally spaced.
Solving bigger problems¶
Up to now, we've only really solved quite small problems, and therefore haven't really had to worry about tuning the solver. As we increase the size of the problem we're solving, the direct solver approach will no longer be good enough. Firedrake uses PETSc to provide solvers, and uses PETSc solver parameters to control them.
Let's dive straight in. We'll write a function that solves the same elasticity problem, but takes parameters for the number of cells in the mesh, as well as a dictionary of solver options.
def solve_elasticity(nx, ny, options=None, **kwargs):
length = 1
width = 0.2
mesh = RectangleMesh(nx, ny, length, width)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
rho = Constant(0.01)
g = Constant(1)
f = as_vector([0, -rho*g])
mu = Constant(1)
lambda_ = Constant(0.25)
Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor
bc = DirichletBC(V, Constant([0, 0]), 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx
uh = Function(V)
solve(a == L, uh, bcs=bc, solver_parameters=options, **kwargs)
return uh
The problem is solved with a Krylov method, so let's limit ourselves to just 100 iterations. Moreover, it is symmetric positive definite, so let's use the conjugate gradient method.
Let's not worry about preconditioning for now and set "pc_type": "none"
.
uh = solve_elasticity(100, 100, options={"ksp_max_it": 100, "ksp_type": "cg", "pc_type": "none"})
--------------------------------------------------------------------------- ConvergenceError Traceback (most recent call last) Cell In[10], line 1 ----> 1 uh = solve_elasticity(100, 100, options={"ksp_max_it": 100, "ksp_type": "cg", "pc_type": "none"}) Cell In[9], line 20, in solve_elasticity(nx, ny, options, **kwargs) 17 L = dot(f, v)*dx 19 uh = Function(V) ---> 20 solve(a == L, uh, bcs=bc, solver_parameters=options, **kwargs) 21 return uh File petsc4py/PETSc/Log.pyx:188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func() File petsc4py/PETSc/Log.pyx:189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func() File /home/firedrake/firedrake/src/firedrake/firedrake/adjoint_utils/solving.py:57, in annotate_solve.<locals>.wrapper(*args, **kwargs) 54 tape.add_block(block) 56 with stop_annotating(): ---> 57 output = solve(*args, **kwargs) 59 if annotate: 60 if hasattr(args[1], "create_block_variable"): File /home/firedrake/firedrake/src/firedrake/firedrake/solving.py:128, in solve(*args, **kwargs) 126 # Call variational problem solver if we get an equation 127 if isinstance(args[0], ufl.classes.Equation): --> 128 _solve_varproblem(*args, **kwargs) 129 else: 130 # Solve pre-assembled system 131 return _la_solve(*args, **kwargs) File /home/firedrake/firedrake/src/firedrake/firedrake/solving.py:164, in _solve_varproblem(*args, **kwargs) 157 # Create solver and call solve 158 solver = vs.LinearVariationalSolver(problem, solver_parameters=solver_parameters, 159 nullspace=nullspace, 160 transpose_nullspace=nullspace_T, 161 near_nullspace=near_nullspace, 162 options_prefix=options_prefix, 163 appctx=appctx) --> 164 solver.solve() 166 # Solve nonlinear variational problem 167 else: 168 if eq.rhs != 0: File petsc4py/PETSc/Log.pyx:188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func() File petsc4py/PETSc/Log.pyx:189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func() File /home/firedrake/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py:89, in NonlinearVariationalSolverMixin._ad_annotate_solve.<locals>.wrapper(self, **kwargs) 86 tape.add_block(block) 88 with stop_annotating(): ---> 89 out = solve(self, **kwargs) 91 if annotate: 92 block.add_output(self._ad_problem._ad_u.create_block_variable()) File /home/firedrake/firedrake/src/firedrake/firedrake/variational_solver.py:286, in NonlinearVariationalSolver.solve(self, bounds) 284 work.copy(u) 285 self._setup = True --> 286 solving_utils.check_snes_convergence(self.snes) 288 # Grab the comm associated with the `_problem` and call PETSc's garbage cleanup routine 289 comm = self._problem.u.function_space().mesh()._comm File /home/firedrake/firedrake/src/firedrake/firedrake/solving_utils.py:138, in check_snes_convergence(snes) 136 else: 137 msg = reason --> 138 raise ConvergenceError(r"""Nonlinear solve failed to converge after %d nonlinear iterations. 139 Reason: 140 %s""" % (snes.getIterationNumber(), msg)) ConvergenceError: Nonlinear solve failed to converge after 0 nonlinear iterations. Reason: DIVERGED_LINEAR_SOLVE
Oh no! We didn't get a solution. This is because the linear system is ill conditioned. Fortunately, good preconditioning strategies exist for elasticity problems. For example, smoothed aggregation multigrid. We can access PETSc's implementation of this scheme, GAMG, by specifying the appropriate preconditioner. In order to reduce the verbosity of the linear solver progress monitoring, we use "ksp_converged_reason": None
, instead of "ksp_monitor": None
.
uh = solve_elasticity(200, 200, options={"ksp_type": "cg",
"ksp_max_it": 100,
"pc_type": "gamg",
"mg_levels_pc_type": "sor",
"mat_type": "aij",
"ksp_converged_reason": None})
Linear firedrake_2_ solve converged due to CONVERGED_RTOL iterations 86
This is still not ideal, taking quite a few iterations. It turns out that for smoothed aggregation to work well, the preconditioner needs access to the near nullspace of the operator. That is, the null modes of the operator if no Dirichlet conditions were applied. For elasticity, these are the rigid body modes of translation and rotation. We must build these and supply them to the solver. To do so, we must create a VectorSpaceBasis
:
def solve_elasticity(nx, ny, options=None, **kwargs):
length = 1
width = 0.2
mesh = RectangleMesh(nx, ny, length, width)
V = VectorFunctionSpace(mesh, "CG", 1)
rho = Constant(0.01)
g = Constant(1)
f = as_vector([0, -rho*g])
mu = Constant(1)
lambda_ = Constant(0.25)
Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor
bc = DirichletBC(V, Constant([0, 0]), 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx
# create rigid body modes
x, y = SpatialCoordinate(mesh)
b0 = Function(V)
b1 = Function(V)
b2 = Function(V)
b0.interpolate(Constant([1, 0]))
b1.interpolate(Constant([0, 1]))
b2.interpolate(as_vector([-y, x]))
nullmodes = VectorSpaceBasis([b0, b1, b2])
# Make sure they're orthonormal.
nullmodes.orthonormalize()
uh = Function(V)
solve(a == L, uh, bcs=bc, solver_parameters=options, near_nullspace=nullmodes)
return uh
With this done, the problem is solved in a reasonably small number of Krylov iterations. This does require that the GAMG smoothing is switched back to the default of Jacobi, rather than SOR as above.
uh = solve_elasticity(200, 200, options={"ksp_type": "cg",
"ksp_max_it": 100,
"pc_type": "gamg",
"mat_type": "aij",
"ksp_converged_reason": None})
Linear firedrake_3_ solve converged due to CONVERGED_RTOL iterations 41