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!