Setting discontinuous function correctly to a DG Function variable

Question

I have a DG function

mesh = UnitSquareMesh(50,50)
V = FunctionSpace(mesh,'DG',1)
u = Function(V)

I want to initialize this to a step function which is 0 for x < 1/2 and 1 for x > 1/2

My mesh edges lie on x=1/2. See this code

from dolfin import *

mesh = UnitSquareMesh(50,50)
V = FunctionSpace(mesh,'DG',1)

foo = Expression("0.0 + 1.0*(x[0] > 0.5)")

u = interpolate(foo, V)
print "Interpolation = ", assemble(inner(grad(u),grad(u))*dx)

u = project(foo, V)
print "Projection = ", assemble(inner(grad(u),grad(u))*dx)

The gradient should be zero but it does not come out like that. Atleast with projection I expected to get zero gradients.

Is there a way to set the function correctly ?

Answer 1

I came up with a simple trick

from dolfin import *

mesh = UnitSquareMesh(50,50)
V0 = FunctionSpace(mesh,'DG',0)
V = FunctionSpace(mesh,'DG',1)

foo = Expression("0.0 + 1.0*(x[0] > 0.5)")

u0 = interpolate(foo, V0)
u = interpolate(u0, V)
print "Interpolation = ", assemble(inner(grad(u),grad(u))*dx)

The function u has zero gradient on each cell.

Answer 2

The problem is that FEniCS silently interpolates Expressions to CG1 if you do not specify anything. This will give a warning in the next version which I believe is appropriate. Something like this should work (not tested):

foo = Expression("0.0 + 1.0*(x[0] > 0.5)", element=V.ufl_element())

Possibly also (not tested)

foo = Expression("0.0 + 1.0*(x[0] > 0.5)", degree=0)

Comments

It does not give zero gradients even if I use this approach

from dolfin import *

mesh = UnitSquareMesh(50,50)
V = FunctionSpace(mesh,'DG',1)

foo = Expression("0.0 + 1.0*(x[0] > 0.5)",element=V.ufl_element())

u = interpolate(foo, V)
print "Interpolation = ", assemble(inner(grad(u),grad(u))*dx)

u = project(foo, V)
print "Projection = ", assemble(inner(grad(u),grad(u))*dx)

Both interpolate and project give a value of 50 for the gradient norm**2.

---

Interpolate will evaluate the Expression at the vertices for DG1. Which gives the same vertex values for cells at both sides of the discontinuity. If you give in a DG0 element to the expression then interpolate will only evaluate your expression in cell centers, and you will avoid the problem since the results will be exactly 1 or exactly 0 in each dof. You do not have access to the cell in the expression as far as I know, so an expression will inherently be single valued for any coordinate.

If you go via CG1 (which happens if you specify nothing to Expression), then there may be a value at the cell centers for the cells adjacent to the discontinuity which is why a DG0 interpolation of a CG1 expression may give wrong results even if the expression should be a step function.

To summarize: for a step function you need to specify to the expression that it represents a DG0 function.