v('+'), v('-'), and jump(v) give wrong values in 1D

Question

In the following program I define the piecewise constant function v which is equal to 1 on
(0, 1/2) and to 0 on (1/2, 1), and then I print out v('+'), v('-'), and jump(v) at the point x = 1/2. In common usage these should be 0., 1., and -1., respectively. But FEniCS reports 1., 0.0, 1.0. It seems that FEniCS reports the value from the left for v('+') and the value from the right for
v('-'), opposite of the normal thing. Then the jump, which is v('+') - v('-') has the wrong sign.

Have I misunderstood, or is this a bug (Dolfin version 1.4.0)?

from dolfin import *
from numpy import array
mesh = UnitIntervalMesh(2)
V = FunctionSpace(mesh, 'DG', 0)
v = Function(V)
v.vector()[:] = array([1., 0.])
F = v('+')*dS
print ( assemble( v('+')*dS ) , assemble( v('-')*dS ) , assemble( jump(v)*dS ) )

Answer 1

Hi,

I don't think that there is any bug

If you look at that question, it is said:

"n('+') is the outward unit normal vector from the cell on the + side.
It's not possible to say just from a picture of a mesh which direction
this will be. What is sure is that n('+') = -n('-')"

If you assume that n('+') has been chosen by dolfin to be (1.0,) - i.e. pointing to the right, then v('+') is the value on the left of the interface and v('-') on the right, which in your case is respectively 1 and 0.

Then the jump of v is indeed 1.0

Comments

Thanks for your comments and the pointer to the relevant previous discussion. I guess
the outcome is that for a function of 1 variable, one cannot count on v('+') representing
the limit of evaluating v from the left or from the right. So what if someone needs
to access the value from the right (to properly upwind, for instance). How does one code that?

---

If you look at demo_dg-advection-diffusion.py, it shows how to upwind:

In your case, if you want to have a direction of flow from the right to the left, then you could have the following code:

v   = TestFunction(V_dg)
phi = TrialFunction(V_dg)
...
u = as_vector((-1.0,))
n = FacetNormal(mesh)
un = (dot(u, n) + abs(dot(u, n)))/2.0
...
a  = dot(grad(v),- u*phi)*dx # volumetric term
a += dot(jump(v), un('+')*phi('+') - un('-')*phi('-') )*dS  # interior facet term
a += dot(v, un*phi)*ds # exterior facet term

It does not really matter how dolfin decides that n will the (1.0,) or (-1.0,):
- if n = (1.0,), then un('+') = 0 and un('-') = 1 so the value selected on the interior facet is phi('-') i.e. the value on the right
- if n = (-1.0,), then un('+') = 1 and un('-') = 0 so the value selected on the interior facet is phi('+') i.e. the value on the right

Same thing for the exterior facets.

Does it make sense?

---

Yes, got it. Many thanks!