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 = form(inner(e_W, e_W) * ufl.dx + inner(nabla_grad(e_W), nabla_grad(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_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):
return 0.5*(nabla_grad(w) + nabla_grad(w).T)
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])
b_vec.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
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
