Projection degree error-norm

Question

Hello.

In the following code, where I try to take the projection of a piecewise constant function and explore some different error-norm measures, one of the norm-calculations does not work when I choose degree 1 in the the line with projection. With degree 2 it works.
I was wondering if anybody knows why this is? The first code section below shows the norm-calculation that does not work (the norm is manually calculated).

from dolfin import *

N = 100
mesh =UnitIntervalMesh(N)
V = FunctionSpace(mesh, "Lagrange", 1)
u_ex = Expression('(x[0] <= 1./3 || x[0] >= 2./3) ? f : g', f=0., g=1., degree=1) 
u_ex_int = interpolate(Expression('(x[0] <= 1./3 || x[0] >= 2./3) ? f : g', f=0., g=1., degree=1), V)
h = project(Expression('(x[0] <= 1./3 || x[0] >= 2./3) ? f : g', f=0., g=1., degree=1), V)

e = h - u_ex_int
error_L2_manual = sqrt(assemble(inner(e, e)*dx)) 
print "L2 error_norm manual, using u_ex_int  as exact  %e " % (error_L2_manual)

The above code gives the error messages:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 699, in runfile
    execfile(filename, namespace)
  File "/usr/lib/python2.7/dist-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 81, in execfile
    builtins.execfile(filename, *where)
  File "/home/....py", line 22, in <module>
    error_L2_manual = sqrt(assemble(inner(e, e)*dx)) 
  File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 577, in sqrt
    return _mathfunction(f, Sqrt)
  File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 569, in _mathfunction
    r = cls(f)
  File "/usr/lib/python2.7/dist-packages/ufl/mathfunctions.py", line 86, in __new__
    return FloatValue(math.sqrt(float(argument)))
ValueError: math domain error

If I instead use the built-in expression for error-norm, in the code below, it works also for degree 1 in the projection command.

e2 = errornorm(u_ex_int, h, 'L2')
print "L2 error_norm error_norm with u_interpolpolated as exact  %e " % (e2)

Answer 1

The assembled value error_L2_manual is negative (maybe due to a round-off error). You can calculate the exact value of the difference e in this way:

e = Function(V)
e.vector()[:] = h.vector() - u_ex_int.vector()

Comments

The commands you gave above showed that "e" was basically equal to zero, so I guess you're right about round-off errors causing trouble. Thanks for helping!

Choosing degree 1 in the projection command gives the exact solution it seems, while choosing degree 2 gives some sort of Gibb's phenomenon. Not sure what degree 2 really implies, when the function space is of degree 1. Guess I would have to read up on this.