State of the art of higher order derivatives of test functions

NGSolve
2

Hi Vincent,
you can interpolate (element-wise) nodes of an expression tree into another fe-space, which allows to build higher order (element-wise!) derivatives:

from ngsolve import *

mesh = Mesh(unit_square.GenerateMesh(maxh=2))
fes = H1(mesh, order=3)
u,v = fes.TnT()

fes2 = L2(mesh, order=2)**2

def Hesse(u):
    return u.Operator("hesse")

bf2a = BilinearForm(InnerProduct(Hesse(u), Hesse(v)) * dx).Assemble()
print (bf2a.mat.DeleteZeroElements(1e-10))

bf2b = BilinearForm(InnerProduct(grad(Interpolate(grad(u), fes2)), grad(Interpolate(grad(v), fes2))) * dx).Assemble()
print(bf2b.mat.DeleteZeroElements(1e-10))

fes3 = (L2(mesh, order=1)**2)**2

def Diff3a(u):
    return grad(Interpolate(grad(Interpolate(grad(u), fes2)), fes3))
def Diff3b(u):
    return grad(Interpolate(Hesse(u), fes3))


bf3a = BilinearForm(InnerProduct(Diff3a(u), Diff3a(v))*dx).Assemble()
print(bf3a.mat.DeleteZeroElements(1e-10))

bf3b = BilinearForm(InnerProduct(Diff3b(u), Diff3b(v))*dx).Assemble()
print(bf3b.mat.DeleteZeroElements(1e-10))

bf2a and bf2b provide the same matrices, as well as bf3a and bf3b. This feature is not heavily tested, better to double-check your results.

Joachim