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?

[apalha]

We tried this idea in this way:

Jump_gradu = jump(dot(grad(u_h), n))
Vjump = FunctionSpace(mesh, "HDiv Trace", p)
jump_vals = project(Jump_gradu, Vjump)
linf_error = max(abs(jump_vals.dat.data_ro))

But we get the following error at the project line:

"ValueError: Cannot restrict an expression twice."

But if we change this line

Jump_gradu = jump(dot(grad(u_h), n))

to

Jump_gradu = dot(grad(u_h), n)

it does not generate any error and we see that this quantity converges to zero. Therefore it could be the jump (or half of the jump), but we are not sure.

Are we doing it correctly?

[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?

[cpraveen]

Isnt

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