# Problems with convergence of the elasticity equation

This is an attempt at verifying the mesh convergence results for the Elasticity equation (Ch.42 Finite Elements II: Galerkin Approximation, Elliptic and Mixed PDEs Finite Elements II: Galerkin Approximation, Elliptic and Mixed PDEs - Archive ouverte HAL or home page of Alexandre Ern)
using a manufactured problem.

The problem is pretty straightforward and the implementation follows the demos and tutorials very closely. However the expected rates of convergence are not observed.
Any feedback, pointers to fix this are greatly appreciated. Thanks.

Since it is a check for convergence, unfortunately the code looks lengthy than a typical MWE.

Checked MWE with the docker image `dolfinx/dolfinx:v0.6.0-r1`

(note: I have checked with many different manufactured problems, I can post if you would like to have a look:
ones with nonhomogeneous boundary, and trigonometric functions insead of 2nd degree polynomials)

Imports and error computations:

``````import numpy as np
import pandas as pd

from dolfinx.mesh import (create_unit_square, CellType, GhostMode,
locate_entities_boundary)

from dolfinx.fem import (Function, FunctionSpace,
dirichletbc, form, locate_dofs_topological)
from dolfinx.fem import assemble_scalar
from dolfinx import fem

import ufl
from ufl import inner, dot, nabla_grad, nabla_div, Identity, dx

from mpi4py import MPI
from petsc4py import PETSc

pd.set_option('display.float_format', '{:1.6e}'.format)

def error_L2(uh, u_ex, W=None, degree_raise=3):
msh = uh.function_space.mesh
# Interpolate approximate solution
u_W = Function(W)
u_W.interpolate(uh)

# Interpolate exact solution
u_ex_W = Function(W)
u_ex_W.interpolate(u_ex)

# Compute the error in the higher order function space
e_W = Function(W)
e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

# Integrate the error
error = form(inner(e_W, e_W) * dx)
error_local = assemble_scalar(error)
error_global = msh.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)

def error_H1(uh, u_ex, W=None):
msh = uh.function_space.mesh
# Interpolate approximate solution
u_W = Function(W)
u_W.interpolate(uh)

# Interpolate exact solution
u_ex_W = Function(W)
u_ex_W.interpolate(u_ex)

# Compute the error in the higher order function space
e_W = Function(W)
e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

# Integrate the error
error_local = assemble_scalar(error)
error_global = msh.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)

def error_max(uh, u_ex, W=None, degree_raise=3):
msh = uh.function_space.mesh
# Interpolate approximate solution
u_W = Function(W)
u_W.interpolate(uh)

# Interpolate exact solution
u_ex_W = Function(W)
u_ex_W.interpolate(u_ex)

# Compute the error in the higher order function space
e_W = Function(W)
e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

error_max = np.max(np.abs(e_W.x.array))
error_max = msh.comm.allreduce(error_max, op=MPI.MAX)
return error_max

def epsilon(w):

def sigma(w, dim, lamda, mu):
return lamda*nabla_div(w)*Identity(dim) + 2*mu*epsilon(w)
``````

Exact solution and data are

``````class u_expression():
def __init__(self, U=[20.0, 10.0]):

self.U = U

def eval(self, x):
return (-self.U[0]*x[0]*(1-x[0])*x[1]*(1-x[1]),
-self.U[1]*x[0]*(1-x[0])*x[1]*(1-x[1]))

class f_expression():
def __init__(self, U=[20.0, 10.0], LAMBDA=20.0, MU=10.0):

self.U = U
self.LAMBDA = LAMBDA
self.MU = MU

def eval(self, x):
return (- 2*self.U[0]*(1-x[1])*x[1]*(self.LAMBDA+2*self.MU)
+ self.U[1]*((1-2*x[0])*(1-2*x[1])*(self.LAMBDA + self.MU) - 2*self.MU*(1-x[0])*x[0]),
- 2*self.U[1]*(1-x[0])*x[0]*(self.LAMBDA+2*self.MU)
+ self.U[0]*((1-2*x[0])*(1-2*x[1])*(self.LAMBDA + self.MU) - 2*self.MU*(1-x[1])*x[1]),
)
``````
``````def solve(mesh=None, degree=2, boundary=None, tags_u=None):

fdim = mesh.topology.dim - 1

element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), degree)
V = FunctionSpace(mesh, element)

# Solution Function
u_n = Function(V, name="u_n")

# boundary data
u_D = Function(V)
u_D.interpolate(u_e.eval)
facets = locate_entities_boundary(mesh, fdim, boundary)
dofs = locate_dofs_topological(V, fdim, facets)
bcs = [fem.dirichletbc(value=u_D, dofs=dofs)]

f = Function(V)
f.interpolate(f_e.eval)

# Assembly of stiffness matrix
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

block = inner(sigma(u, 2, LAMBDA, MU), epsilon(v)) * dx
rhs = inner(f, v) * dx

A_form = form(block)
L_form = form(rhs)

A_mat = fem.petsc.assemble_matrix(A_form, bcs=bcs)
A_mat.assemble()

b_vec = fem.petsc.create_vector(L_form)

# Direct solver
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A_mat)
solver.setType("preonly")
solver.getPC().setType("lu")

with b_vec.localForm() as loc_b:
loc_b.set(0)

fem.petsc.assemble_vector(b_vec, L_form)

fem.petsc.apply_lifting(b_vec, [A_form], bcs=[bcs])

fem.petsc.set_bc(b_vec, bcs)

solver.solve(b_vec, u_n.vector)
u_n.x.scatter_forward()

return u_n, solver
``````

## Convergence check

``````boundary = lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 1) | np.isclose(x[1], 0) | np.isclose(x[1], 1)

ns_list = [2**i for i in range(3,8)]

degree = 2
LAMBDA = 20.0
MU = 10.0
U = [20.0, 10.0]
u_e = u_expression(U=U)
f_e = f_expression(U=U, LAMBDA=LAMBDA, MU=MU)

df = pd.DataFrame([], columns=["ns", "error_L2", "error_H1", "error_max"])

for cnt, ns in enumerate(ns_list):
mesh = create_unit_square(MPI.COMM_WORLD,
ns, ns,
CellType.triangle, GhostMode.none)

element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), degree+3)
V_ex = FunctionSpace(mesh, element)

# Function
u_ref = Function(V_ex)
u_ref.interpolate(u_e.eval)

u_n, solver = solve(mesh, degree=degree, boundary=boundary)

df.loc[cnt, "ns"] = ns
df.loc[cnt, "error_L2"] = error_L2(u_n, u_ref, V_ex)
df.loc[cnt, "error_H1"] = error_H1(u_n, u_ref, V_ex)
df.loc[cnt, "error_max"] = error_max(u_n, u_ref, V_ex)

print(df)
``````

The result:

``````ns	error_L2	error_H1	error_max
0	8	9.993293e-02	4.757606e-01	1.288192e-01
1	16	1.000108e-01	4.743987e-01	1.289061e-01
2	32	1.000236e-01	4.744044e-01	1.289178e-01
3	64	1.000255e-01	4.744200e-01	1.289195e-01
4	128	1.000258e-01	4.744235e-01	1.289197e-01
``````

I ran into the classic problem, I guess I was too hasty to put my question here.
I had written down the Navier-Cauchy equations wrong all along since a long long time, which worked for a manufactured problem with certain symmetry.

I overlooked checking that up. I guess that should have been the first step instead of hallucinating about how the dofs are organized in VectorFunctionSpaces .

``````class f_expression():
def __init__(self, U=[20.0, 10.0], LAMBDA=20.0, MU=10.0):

self.U = U
self.LAMBDA = LAMBDA
self.MU = MU

def eval(self, x):
return (- 2*self.U[0]*((1-x[1])*x[1]*(self.LAMBDA+2*self.MU) + self.MU*x[0]*(1-x[0]))
+ self.U[1]*(1-2*x[0])*(1-2*x[1])*(self.LAMBDA + self.MU),
- 2*self.U[1]*((1-x[0])*x[0]*(self.LAMBDA+2*self.MU) + self.MU*x[1]*(1-x[1]))
+ self.U[0]*(1-2*x[0])*(1-2*x[1])*(self.LAMBDA + self.MU),
)

``````