Evaluate the jump discontinuity along the interface of two elements

[apalha]

Dear all,

We need to evaluate the jump discontinuity along the interface of two elements, i.e., evaluate the following expression at points

u_h = FunctionSpace(mesh, "DG L2", p)
n = FacetNormal(mesh)
jump(dot(grad(u_h), n))

We are stuck at evaluating the last expression. We can create it, but it is only usable in an assembly, not evaluation.

Is there any way to evaluate an expression on the + and - side of an element?

Thank you for your help.

Kind regards,

-artur palha

[dham]

I'm not sure that the term "evaluate" is well-defined here. jump(dot(grad(u_h), n)) is a function from the skeleton of the mesh to R. What would make sense is to project it into a finite element space which is defined on the skeleton. The space for that is HDiv Trace which is effectively DG on the facets and so is defined exactly where that expression is defined. Is that what you are after?

[dham]

project when projecting into a trace space takes the average on the interior facets. I think in this case you need to formulate the projection yourself. You also need to make a decision about what you want to have happen on the exterior facets of the mesh where jump isn't defined. One answer would be to simply impose a zero DirichletBC everywhere on the boundary.

However, I'm not sure if this is going to get you what you want, because we don't expose an API for evaluating trace spaces at arbitrary points. Maybe we shoud take a step back here and ask what the actually mathematical problem that you are trying to solve is?

[apalha]

In this case, if we use

dot(grad(u_h), n)

We will get the 1/2 of the jump, right? Since the normal will take a different sign on each side of the interface, making the average effectively a difference. Indeed, we are using homogeneous boundary conditions to avoid the issue you mention.

To clarify our objective. We are working on a new approach that, while using standard C0 spaces, provides an improved error in the jump of the derivative across interfaces (something akin to approximate C1 spaces). We want to compare our results to other existing approaches. One such approach is a mixed approach that we implemented in Firedrake.

We are trying to compute an H\infty norm of the jump on the whole skeleton. As a start we are computing a crude estimation with

linf_error = max(abs(jump_vals.dat.data_ro))

in the example above.

[dham]

Ah, this misinterprets facet integrals in Firedrake. The FacetNormal is relative to the intrinsic orientation of the facet. It is not a cell outward normal. So no, you will get the average of the gradient evaluated from each side. To evaluate the jump you really need to take the jump operator.

I think you can estimate your h_inf norm by just projecting the jump into a fairly high degree trace space and then taking the approach you are taking. Obviously that's not a guaranteed h_inf (for grad(u_h) of degree greater than 1) but it does at least converge as degree increases.

[apalha]

Thanks for clarifying that.

The issue is that

Jump_gradu = jump(dot(grad(u_h), n))
Vjump = FunctionSpace(mesh, "HDiv Trace", p)
jump_vals = project(Jump_gradu, Vjump)

Does not work because it stats that

"ValueError: Cannot restrict an expression twice."

How can we project the jump operator?

[dham]

I was suggesting above that you just do it yourself:

from firedrake import *
m = UnitSquareMesh(1,1)
fs = FunctionSpace(m, "CG", 1)
fs_t = FunctionSpace(m, "HDiv Trace", 0)
t = Function(fs_t)
v = TestFunction(fs_t)
f = Function(fs)
x = SpatialCoordinate(m)
f.interpolate(x[0]**2 * x[1]**2)
bc= DirichletBC(fs_t, 0.0, "on_boundary")
n = FacetNormal(m)

solve(inner(avg(t)-jump(grad(f), n), avg(v))*dS == 0, t, bcs=bc)

print(t.dat.data)

[cpraveen]

Isnt

avg(dot(grad(u),n)) = 0.5 * jump in normal derivative of u

[dham]

Perhaps I'm confused. I thought n did not change sign in UFL but maybe I have that definition wrong.

[J15525]

I work together with @apalha and we made a simple example to check the definitions.

We started with a simple function (on quads, since we are more familiar with those), using the following script.

from firedrake import *

import numpy as np
np.set_printoptions(precision=2, suppress=True)   # To reduce visual noise

# Create a simple unit square mesh with 2x1 elements, using quadrilaterals.
mesh = UnitSquareMesh(2, 1, quadrilateral=True)

# Create a function in the discontinuous space.
Wh = FunctionSpace(mesh, "DQ", 1)
vals = np.zeros(Wh.dim())
vals[0] = 1.0
vals[1] = 1.0
vals[6] = 1.0
vals[7] = 1.0
u_h = Function(Wh, val=vals, name="u_h")

This creates an element-wise linear (in x) function which matches at the interface. Its gradient (in the x-direction) is constant on each element, where the value of the gradient on the left element is about 3.5 and on the left -3.5. This is also visible in the two figures (first the function, then the gradient).

If we now compute the 'jump' of the gradient like this

# Compute (half) the jump of the gradient across facets.
n = FacetNormal(mesh)
Jump_gradu = dot(grad(u_h), n)
Vjump = FunctionSpace(mesh, "HDiv Trace", 0)
jump_vals = project(Jump_gradu, Vjump)
print(jump_vals.dat.data_ro[0])
print(jump_vals.dat.data_ro)

the output is

3.4641016151377553
[ 3.46  0.   -3.46  0.   -3.46  0.    0.  ]

Note that the actual jump is about 7.0, while the average is 0. The other two non-zero entries are related to the boundary facets, which we ignore here (as @dham mentioned, these are probably best set to 0, though that would be too much for this example. In our original case, this is already the case.). So the numbering of facets seems to first number internal facets, then boundary facets.

---

If we instead define the function on the same mesh as

Wh = FunctionSpace(mesh, "DQ", 1)
vals = np.zeros(Wh.dim())
vals[0] = 1.0
vals[1] = 1.0
vals[6] = -1.0
vals[7] = -1.0
u_h = Function(Wh, val=vals, name="u_h")

the function and its gradient look like

If we then compute the 'jump' of the gradient like this again

# Compute (half) the jump of the gradient across facets.
n = FacetNormal(mesh)
Jump_gradu = dot(grad(u_h), n)
Vjump = FunctionSpace(mesh, "HDiv Trace", 0)
jump_vals = project(Jump_gradu, Vjump)
print(jump_vals.dat.data_ro[0])
print(jump_vals.dat.data_ro)

we instead get

-3.3053706985194236e-17
[-0.    0.   -3.46  0.    3.46 -0.    0.  ]

Note that the actual jump is 0, while the average is -3.5. The other two non-zero entries are again related to the boundary facets, one of which has now changed sign (as expected).

---

In conclusion:

• It appears that the normals do indeed change sign. Which is why example 1 leads to a computed jump and example two to zero.
• This also means that what we computed is actually half the jump, in line with @cpraveen's answer.

Thank you both for your help in figuring this out.