Suggestions for increasing the accuracy of assemble(), or a different approach for integration

Question

I am trying to find a way to integrate a second order expression with "sufficient" accuracy using assemble(). The simplified working examples below demonstrate the issue I am having. I would have expected an answer much closer to zero (I was hopping for something less than 1e-6). As shown below, the situation improves when refining the mesh, but I am perplexed by the odd way it improves. I am new to FENICS, and would very much appreciate any insight on how to address this (I have been unsuccessful finding it elsewhere).

from __future__ import print_function
from __future__ import division
from dolfin import *
import mshr as ms
import numpy as np

seg = 100
circle = ms.Circle(Point(0,0), 1, seg)
mesh = ms.generate_mesh(circle, seg/3.14)

dx = Measure('dx', domain=mesh)

expr = Expression('2*(x[0]*x[0] + x[1]*x[1]) - 1', degree=2)

int1 = assemble(expr*dx)

V = FunctionSpace(mesh, "CG", 2)
expr_i = interpolate(expr, V)
int2 = assemble(expr_i*dx)

V = FunctionSpace(mesh, "CG", 1)
expr_i = interpolate(expr, V)
int3 = assemble(expr_i*dx)

print('int1={} int2={} int3={}'.format(int1, int2, int3))

Result:
int1=-0.00206504578602 int2=-0.00206504578602 int3=0.00127584823819

Below I just refine the mesh by changing by increasing seg from 100 to 400

seg = 400
circle = ms.Circle(Point(0,0), 1, seg)
mesh = ms.generate_mesh(circle,seg/3.14)

dx = Measure('dx', domain=mesh)

expr = Expression('2*(x[0]*x[0] + x[1]*x[1]) - 1', degree=2)

int1 = assemble(expr*dx)

V = FunctionSpace(mesh, "CG", 2)
expr_i = interpolate(expr, V)
int2 = assemble(expr_i*dx)

V = FunctionSpace(mesh, "CG", 1)
expr_i = interpolate(expr, V)
int3 = assemble(expr_i*dx)

print('int1={} int2={} int3={}'.format(int1, int2, int3))

Result:
int1=-0.000129184850435 int2=-0.000129184850435 int3=7.73547968041e-05

With seg = 500 I get:
int1=-0.000129184850435 int2=-8.26801403224e-05 int3=5.02225121338e-05

Answer 1

The slow rates of convergence you see are arising due to committing the variational crime that your computation domain is not an exact representation of the geometry. I.e. $$\Omega_h \neq \Omega.$$

FEniCS has some support for curved elements with quadratic boundary representation.

Consider the following:

from dolfin import *
import mshr as ms
import numpy as np

np.set_printoptions(3)

n = 5
boundary_face_order = 2

p1_err = np.zeros(n)
p2_err = np.zeros(n)
p3_err = np.zeros(n)
h = np.zeros(n)

for j in range(n):
  mesh = UnitDiscMesh(mpi_comm_world(), 2**(j+1), boundary_face_order, 2)

  expr = Expression('2*(x[0]*x[0] + x[1]*x[1]) - 1', degree=5)

  V1 = FunctionSpace(mesh, "CG", 1)
  p1_err[j] = abs(assemble(interpolate(expr, V1)*dx))

  V2 = FunctionSpace(mesh, "CG", 2)
  p2_err[j] = abs(assemble(interpolate(expr, V2)*dx))

  V3 = FunctionSpace(mesh, "CG", 3)
  p3_err[j] = abs(assemble(interpolate(expr, V3)*dx))

  h[j] = mesh.hmax()

print "P1 error"
print p1_err
print "P1 rate of convergence"
print np.log(p1_err[1:]/p1_err[:-1])/np.log(h[1:]/h[:-1])

print "P2 error"
print p2_err
print "P2 rate of convergence"
print np.log(p2_err[1:]/p2_err[:-1])/np.log(h[1:]/h[:-1])

print "P3 error"
print p3_err
print "P3 rate of convergence"
print np.log(p3_err[1:]/p3_err[:-1])/np.log(h[1:]/h[:-1])

You should observe the following output when boundary_face_order = 1

P1 error
[ 0.317  0.087  0.022  0.006  0.001]
P1 rate of convergence
[ 2.132  2.078  2.042  2.022]
P2 error
[ 0.134  0.035  0.009  0.002  0.001]
P2 rate of convergence
[ 2.192  2.094  2.047  2.023]
P3 error
[ 0.134  0.035  0.009  0.002  0.001]
P3 rate of convergence
[ 2.192  2.094  2.047  2.023]

and the following when boundary_face_order = 2

P1 error
[ 0.388  0.112  0.03   0.008  0.002]
P1 rate of convergence
[ 2.041  2.021  2.011  2.006]
P2 error
[  6.720e-02   8.764e-03   1.111e-03   1.396e-04   1.749e-05]
P2 rate of convergence
[ 3.345  3.149  3.07   3.034]
P3 error
[  6.724e-02   8.768e-03   1.111e-03   1.396e-04   1.749e-05]
P3 rate of convergence
[ 3.346  3.149  3.07   3.034]

You can see that you will only have optimal convergence in your FEM scheme if the order of the boundary representation is equal to the polynomial order of the approximation space. This is a well known result in FEM theory.

Comments

Very helpful. Thank you.