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

183

views

0

```
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)
a = dot(grad(u), grad(v))*dx #new for for a
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

Hi Dave,

I am looking at the error at individual vertices.

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

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

Please login to add an answer/comment or follow this question.

L^{2}error, by includingat the end of the script, it's much larger than 10e-15, even with many more elements in the UnitSquareMesh.