*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 rate 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

In [1]:

```
%matplotlib notebook
import matplotlib.pyplot as plt
from firedrake import *
length = 1
width = 0.2
mesh = RectangleMesh(40, 20, length, width)
```

`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).

In [2]:

```
V = VectorFunctionSpace(mesh, "Lagrange", 1)
```

We need a boundary condition object for $\Gamma_D$.

In [3]:

```
bc = DirichletBC(V, Constant([0, 0]), 1)
```

In [4]:

```
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.

In [5]:

```
def epsilon(u):
return 0.5*(grad(u) + grad(u).T)
def sigma(u):
return lambda_*div(u)*Id + 2*mu*epsilon(u)
```

`"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.

In [6]:

```
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})
```

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$.

In [7]:

```
displaced_coordinates = interpolate(SpatialCoordinate(mesh) + uh, V)
displaced_mesh = Mesh(displaced_coordinates)
```

`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.

In [8]:

```
# NBVAL_IGNORE_OUTPUT
fig, axes = plt.subplots()
triplot(displaced_mesh, axes=axes)
axes.set_aspect("equal");
```

Modify the problem so that the material density $\rho$ is not constant, but rather varies in space. For example, you could try setting $\rho(x, y) = 0.01 + xy$. That is, the material gets denser the futher away from the clamped end you get.

- Hint 1: The values for the components of
`as_vector`

can be arbitrary UFL expressions. - Hint 2: You can get symbolic expressions for the coordinates with
`x, y = SpatialCoordinate(mesh)`

.

In [ ]:

```
```

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.

In [9]:

```
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
```

`"pc_type": "none"`

.

In [10]:

```
# NBVAL_RAISES_EXCEPTION
uh = solve_elasticity(100, 100, options={"ksp_max_it": 100, "ksp_type": "cg", "pc_type": "none"})
```

`"ksp_converged_reason": None`

, instead of `"ksp_monitor": None`

.

In [11]:

```
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 85

`VectorSpaceBasis`

:

In [12]:

```
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
def epsilon(u):
return 0.5*(grad(u) + grad(u).T)
def sigma(u):
return lambda_*div(u)*Id + 2*mu*epsilon(u)
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.

In [13]:

```
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 28

Study what happens to the number of iterations for this last setup as you change the mesh resolution. Try, perhaps, 10, 50, 100, and 200.

In [ ]:

```
```