errornorm and discontinuous solutions

Question

I just figured out a problem in my code that results from using errornorm with discontinuous solutions:

from dolfin import *
mesh = UnitSquareMesh(10, 10)
Q = FunctionSpace(mesh, 'DG', 0)
P = Expression('x[1] > .5')
p = interpolate(P, Q)
print errornorm(P, p)

Above, I interpolate a step-function to a piecewise constant function space and I make sure that the discontinuity is resolved by the mesh. However, executing the code yields

0.0872871560944

which is clearly not the right answer. From the description of errornorm I could find out that it projects both functions to a higher degree discontinuous space, subtracts the coefficients and then square-integrates the resulting difference on the mesh associated with the discrete function space. Now it seems to me that some of the quadrature points are located at the discontinuity and therefore the projection to higher degree discontinuous spaces always smears out the jump on one side of the step. Is it possible to work around this issue?

Answer 1

Choose appropriate element for expression (default is CG1). For instance

elem = FiniteEleement('DG', triangle, 0)
P = Expression('x[1] > .5', element=elem)

Comments

Thanks for your answer. Even when I specify the element, the problem still persists.

However, your comment directed me to a possible solution... The piecewise defined function needs to be defined in the following way, by checking if the cell midpoint is 'above' or 'below' the interface...

from dolfin import *

mesh = UnitSquareMesh(10, 10)
Q = FunctionSpace(mesh, 'DG', 0)

class Step(Expression):
   def __init__(self, mesh):
      self.mesh = mesh
   def eval_cell(self, value, x, ufc_cell):
      cell = Cell(self.mesh, ufc_cell.index)
      value[0] = (cell.midpoint().y() > 0.5)

P = Step(mesh)
p = interpolate(P, Q)

print errornorm(P, p)

I don't know if there is a shorter way, but I am happy with this.

---

Ok. Nevertheless, note that DG0 P is a proper representation (on this matching mesh). The discrepancy is given by the (relatively complicated) logic in site-packages/dolfin/fem/norms.py.

from dolfin import *
mesh = UnitSquareMesh(10, 10)
Q = FunctionSpace(mesh, 'DG', 0)
elem = FiniteElement('DG', triangle, 0)
P = Expression('x[1] > .5', element=elem)
p = interpolate(P, Q)
print errornorm(P, p, degree_rise=0) # returns 0.0