# Orders of convergence too low for simple linear-elastic problem

Dear community,

I am trying to conduct a convergence analysis for a simple two-dimensional linear-elastic problem. I consider the convergence of the L2 error with respect to a high-resolution reference solution. However, the order of convergence I get is much lower than expected. The order of convergence also seems too low to me if you consider that I use a high-resolution approximation instead of an analytical solution to calculate the L2 error.

Here is a brief description of the problem setup:
The domain is rectangular with an edge length of 100 mm. The left edge is clamped (displacement in the x and y directions are 0.0, respectively) and a traction force of 100 N is applied to the right edge in the x direction. The plate is assumed to be in plane stress state.

Code:
I apologize that the code example is so lengthy. A further reduction of the code would have simplified the problem too much and impaired readability. Thanks in advance for reading through the code!

``````import dolfinx
import numpy as np
import ufl
from dolfinx import default_scalar_type
from dolfinx.fem import (
Constant,
Function,
assemble_scalar,
create_nonmatching_meshes_interpolation_data,
dirichletbc,
form,
functionspace,
locate_dofs_topological,
)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, locate_entities_boundary, meshtags
from mpi4py import MPI
from scipy import stats

# Setup
edge_length = 100.0
E = 210000.0
nu = 0.3
# FEM
fem_element_family = "Lagrange"
fem_element_degree = 1
fem_element_shape = (2,)
fem_num_elements_reference = 256
fem_num_elements_analysis = [8, 16, 32]

def solve_problem(num_elements):
# Mesh
mesh = create_rectangle(
comm=MPI.COMM_WORLD,
points=[[0.0, 0.0], [edge_length, edge_length]],
n=[num_elements, num_elements],
cell_type=fem_cell_type,
)

# Function space
element = (fem_element_family, fem_element_degree, fem_element_shape)
func_space = functionspace(mesh, element)
trial_function = ufl.TrialFunction(func_space)
test_function = ufl.TestFunction(func_space)

# Marking facets
tag_left = 0
tag_right = 1
facets_dim = mesh.topology.dim - 1
locate_left_facet = lambda x: np.isclose(x[0], 0.0)
locate_right_facet = lambda x: np.isclose(x[0], edge_length)
left_facets = locate_entities_boundary(mesh, facets_dim, locate_left_facet)
right_facets = locate_entities_boundary(mesh, facets_dim, locate_right_facet)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack(
[np.full_like(left_facets, tag_left), np.full_like(right_facets, tag_right)]
)
sorted_facet_indices = np.argsort(marked_facets)
facet_tags = meshtags(
mesh,
facets_dim,
marked_facets[sorted_facet_indices],
marked_values[sorted_facet_indices],
)

# Variational problem
mu_ = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
# plane stress
lambda_ = 2 * mu_ * lambda_ / (lambda_ + 2 * mu_)

def epsilon(u):

def sigma(u):
return 2.0 * mu_ * epsilon(u) + lambda_ * ufl.tr(epsilon(u)) * ufl.Identity(2)

ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags)

volume_force = Constant(mesh, (default_scalar_type((0.0, 0.0))))
traction_right = Constant(mesh, default_scalar_type((100.0, 0.0)))

lhs = ufl.inner(sigma(trial_function), epsilon(test_function)) * ufl.dx
rhs = ufl.dot(volume_force, test_function) * ufl.dx + ufl.dot(
traction_right, test_function
) * ds(tag_right)

u_left = Constant(mesh, default_scalar_type((0.0, 0.0)))
dofs_facet_left = locate_dofs_topological(func_space, facets_dim, left_facets)
u_bc_left = dirichletbc(u_left, dofs_facet_left, func_space)

problem = LinearProblem(
lhs,
rhs,
bcs=[u_bc_left],
petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
)

return problem.solve()

def L2_error(u_approx, u_exact, degree_raise=3):
# Create higher order function space
degree = u_approx.function_space.ufl_element().degree()
family = u_approx.function_space.ufl_element().family()
mesh = u_approx.function_space.mesh
element = (family, degree + degree_raise, fem_element_shape)
W = functionspace(mesh, element)

# Interpolate approximate solution
u_W = Function(W)
u_W.interpolate(u_approx)

# Interpolate exact solution
u_ex_W = Function(W)
interpolation_data = create_nonmatching_meshes_interpolation_data(
W.mesh._cpp_object, W.element, u_exact.function_space.mesh._cpp_object
)
u_ex_W.interpolate(u_exact, nmm_interpolation_data=interpolation_data)

# 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(ufl.inner(e_W, e_W) * ufl.dx)
error_local = assemble_scalar(error)
error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)

def calculate_convegrence_order(errors, element_sizes):
log_element_sizes = np.log(np.array(element_sizes))
log_errors = np.log(np.array(errors))
slope, _, _, _, _ = stats.linregress(log_element_sizes, log_errors)
return slope

u_exact = solve_problem(fem_num_elements_reference)

element_sizes = []
l2_errors = []

for num_elements in fem_num_elements_analysis:
u_approx = solve_problem(num_elements)
element_sizes.append(edge_length / num_elements)
l2_errors.append(L2_error(u_approx, u_exact))

convergence_order = calculate_convegrence_order(l2_errors, element_sizes)
print(f"Convergence order: {convergence_order}")
``````

Development environment
I run the code in a docker container built from dolfinx/dolfinx:stable. In the container, DOLFINx version 0.7.0 is installed.

Results:
With the above code I got the following convergence orders for the L2 error:

• Triangle elements, element degree = 1: 1.5016
• Triangle elements, element degree = 2: 1.5520
• Quadrilateral elements, element degree = 1: 1.4264
• Quadrilateral elements, element degree = 2: 1.5654

As an example, the following plot shows for triangular elements of order 1 that the calculated error values are well aligned on a line in a double logarithmic plot.

In addition, I observed that the approximated displacement field in y direction is asymmetrical when using triangle elements. The absolute displacement is smaller in the top right-hand corner than in the bottom left-hand corner. This effect disappears as the number of elements increases. Can this observation be explained by the fact that I have used triangular elements?

Since I am not calculating the L2 error with respect to the analytical solution, but a high-resolution approximation, I am aware that I cannot expect perfect convergence orders. Nevertheless, the convergence orders seem to me to be significantly too low.

The results in unexpectedly-high-l2-error-when-comparing-coarse-and-fine-resolution-approximation show for the poisson problem that even if you use a high-resolution approximation to calculate the L2 error you get convergence order very close to the expected ones (2 for linear ansatz functions and 3 for quadratic ansatz functions, respectively).

Does anyone have any idea why my convergence orders deviate so significantly from the expected values? Thank you in advance for your help! Greetings!