Interpolation of Discontinous Function Using DG

1

I am trying to interpolate the characteristic function

\begin{align}f(x) = \begin{cases} 0 & 0 \leq x < 10 \\ 1 & 10 \leq x \leq 25 \\ 0 & 25 < x \leq 400 \end{cases} \end{align}

I use the IntervalMesh() command to generate the required spatial mesh. I try to make sure that the nodes will line up with the discontinuities. So I pick the number of subintervals so that the size of each one is 5 units. With the interval being [0, 400], this means I need 80 subintervals. I then use V as a DG0 function space on this interval mesh. I interpolate f in this space, but it doesn’t want to give me the correct value at 10. It outputs 0.0, not 1.0. I don’t want this to be the case, as I want to increase the degree of V, but this will lead to a piece being a line with nonzero slope due to some ‘overhang’. It is critical to my work that this ‘overhang’ not happen. Any help with this is greatly appreciated.

My minimum working code example is as follows.

from fenics import *
from numpy import *

p1 = Expression(‘1.0’, degree=4)
p2 = Expression(‘0.0’, degree=4)
func = Expression(‘x[0] < 10.0 ? p2 : x[0] > 25.0 ? p2 : p1’, p1=p1, p2=p2, degree=4)

mesh = IntervalMesh(80, 0, 400)

V = FunctionSpace(mesh, ‘DG’, 0)

f = interpolate(func, V)

print('f(9.99999) = ', f(9.99999))
print('f(10.0) = ', f(10.0))
print('f(10.00001) = ', f(10.00001))
print('f(15.0) = ', f(15.0))
print('f(20.0) = ', f(20.0))
print('f(25.0) = ', f(25.0))
print('f(25.00001) = ', f(25.00001))

2

Hello,
I am sorry but evaluating f (being a DG0 function) at x=10 is not well defined since the function is precisely discontinuous. Ignoring fixed machine precision x=10 belongs to neither left or right element.
What happens here is that due to finite arithmetic 10.0 gets tested to belong to the left element, hence f = 0.
But use 10.000000000000001 and then you get f=1.

For DG0 function, you will have the correct behaviour if you look at the vector of degrees of freedom (f.vector().get_local()) since the expression will be evaluated at the cell centers which do not lie on the mesh boundary.
For DG1 or higher it may be trickier. I suggest interpolating your expression first on a DG0 space, then interpolate this function on a DG1 space. You should then obtain the correct values.

1 Like