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_cell_type = dolfinx.mesh.CellType.triangle # dolfinx.mesh.CellType.quadrilateral
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):
return ufl.sym(ufl.grad(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!