Assembling gradient jump

Question

I have a CG2 function on a triangular mesh. I need to check whether all jumps of the gradient over all interior edges are non-negative. Since the jump (in normal direction) of the gradient is linear at the edges, it is sufficient to check the two end-points.

Is there any (easy) possibility to assemble a matrix $A$, such that $A\,x \ge 0$
iff the P2-function represented by $x$ has the above property?

That would be easy if there would be a function space of functions that live on edges, are linear on each edge, and discontinuous. Then one could project the gradient jump into that space along the lines of

 u = TrialFunction("FancyBoundarySpace", triangle, 1)
 v = TestFunction("FancyBoundarySpace", triangle, 1)
 y = Coefficient("CG", triangle, 2)
 a = u*v*dS
 L = jump(grad(y),n)*v*dS

Answer 1

I don't understand what is jump of the gradient, the definition depends on the orientation of the jump. Nevertheless, after you make this clean you can do something like, i.e. test the property against CG1 function:

# Your mesh and function
mesh = UnitSquareMesh(42, 42)
u = Function(FunctionSpace(mesh, 'Lagrange', 2))

# Test the property against CG1
v = TestFunction(FunctionSpace(mesh, 'Lagrange', 1))
n = FacetNormal(mesh)
J = jump(grad(u), n)*v('+')*dS
j = assemble(J)
print j.min()

Comments

Thank you for your answer. This would be a reasonable possibility. However, this are too few inequalities. In my question, we need $2 \, n_e$ inequalities, where $n_e$ is the number of edges. Your code just produces $n_p$, with $n_p$ being the number of vertices.

---

Following fact is implied

Since the jump (in normal direction) of the gradient is linear at the edges, it is sufficient to check the two end-points.

I don't understand why you should have $2n_e$. The jump is the same on both sides of the edge.

---

Yes, the gradient's normal jump is linear on all edges. However, it is discontinuous at the vertices. Hence, we have $2$ constraints per (interior) edge (and not only $1$ per vertex).

---

Yes, now I see the problem. So it seems that it can't be tested against any cell element. You would need facet element ("FancyBoundarySpace" of your snippet) which is not yet supported as far as I know. You need to invent some workaroud.

Maybe getting DG1 projection of grad(y) and then looping around mesh entities and evaluating there. Note that you'd need to use Function.eval method with cell argument to be able distinguishing between + and - side. (For performance, the loop over mesh entities can be rewritten to C++ afterwards.)