Use firedrake to solve the equations of 2D linear elasticity

firedrake_elasticity.py

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

V = VectorFunctionSpace(mesh, "Lagrange", 1)
bc = DirichletBC(V, Constant([0, 0]), 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)

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

plt.figure(figsize=(8, 8))
ax = plt.gca()
ax.set_aspect("equal")
triplot(mesh, axes=ax, interior_kw=dict(alpha=0.05))
quiver(uh, axes=ax)