### Problem using Gmesh - solution to manufactured problem gives errors that are too large.

183
views
0
6 months ago by
from fenics import *  #replace dolfin with fenics
import numpy as np    #use np

# Create mesh from internal method
mesh = UnitSquareMesh(3, 3)

#Create mesh from XML file
#mesh = Mesh("mesh.xml")

#define function space
V = FunctionSpace(mesh, 'Lagrange', 1)

# Define boundary conditions
u_D = Expression('1. + x[0]*x[0] + 2*x[1]*x[1]', degree=2)    #u_D, degree = 2

def u0_boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, u_D, u0_boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)​
Hi I am new to fenics and have been doing the examples.  In the example above, I solved the Poisson equation for a manufactured solution using a mesh created with UnitSquareMesh.  I get the expected error of 10e-15.  Next I created a mesh with Gmesh for a unit square and I generated a msh file.  I used convert-mesh to generate the xml file.  Then I read the xml file and generated the mesh using Mesh(‘name.xml’).  Now, checking the error between the computed solution and the exact solution gives errors in the range of 10e-5.  Can anyone explain what I am doing wrong?  I am using fenics installed with Docker, python3 and Gmesh 3.0 on a machine running Windows 7.
Community: FEniCS Project
What are you measuring as the error?  If I check, e.g.,  $L^2$L2 error, by including

import math
L2error = math.sqrt(assemble(((u-u_D)**2)*dx))
print(L2error)​

at the end of the script, it's much larger than 10e-15, even with many more elements in the UnitSquareMesh.
written 6 months ago by David Kamensky
Hi Dave,

I am looking at the error at individual vertices.

#Dump solution to the screen
u_nodal_values = u.vector()
u_array = u_nodal_values.array()

# Verification
ue = interpolate(u_D, V)
ue_array = ue.vector().array()

print('\n\nCompute error at each vertice\n')
for i in range(len(ue_array)):
u_abs_error = np.abs(u_array[i] - ue_array[i])
print('u({:8.6f},{:8.6f}) = {:6.4f}  {:.8e}'.format(coor[i][0], coor[i][1],
u_array[i], u_abs_error))
​
written 6 months ago by Larry Christman
Okay, it seems that there is a superconvergence phenomenon, where, because of your specific problem and data, on a structured mesh, you happen to get the nodal interpolant.  (This is relatively easy to show for the Poisson problem in 1D with linear elements, and the result can probably be extended to the current case using linearity, since the exact solution is the sum of 1D solutions in each direction.)  However that is not typical of the finite element method.  My guess is that the Gmsh mesh is something unstructured, like this:

http://openfoamwiki.net/index.php/File:Gmsh_1layerPrismsMesh.png

and you no longer get nodal interpolation.  For instance, if you change the mesh in FEniCS to be unstructured, like

from mshr import *
r = Rectangle(Point(0,0),Point(1,1))
mesh = generate_mesh(r,10)​

or you manufacture a more complicated solution with the structured mesh, like

u_D = Expression('1. + x[0]*x[0] + 2*x[1]*x[1] + x[0]*x[0]*x[0]*x[1]', degree=4)

# ...

f = Expression('-6.0*x[0]*x[1] - 6.0',degree=2)

then you will still see convergence with mesh refinement, but the finite element solution will no longer match the exact solution at the nodes to machine precision.

written 6 months ago by David Kamensky
Thank you Dave.  You cleared up something that was really confusing me.
written 6 months ago by Larry Christman