Dear NGSolve community,
I am currently dealing with a function that solves an advectiondominated equation,
which I let NGSolve approximate in L2 space with DG.
Analytically, this function is smooth enough to have a gradient which I want to calculate.
In this context I played a little with the L2 space and came up with behaviour I do not understand,
so I kindly ask for elaboration on this:

Gradients of functions from L2 with order=1 seem to evaluate to zero, but increasing the order
gives nonzero gradients. I would understand that the gradient isn’t implemented since it isn’t expected to exists analytically  although in a discrete sense one could always calculate it elementwise, right? , but am puzzled with the behaviour for higher orders. Is this related with different basis functions? 
I can evaluate L2 functions on vertices and also the gradient of higher order L2 functions although
it isn’t continous across elements (I found discontinuous gradient  Kunena). How can I understand the value at the vertex (which is shared by multiple elements)? 
I found that when doing mesh refinement (code below), the gradient of higher order L2 functions
vanishes. Is this intended and if yes, is there a proper way to work around this?
from ngsolve import *
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))
fes = L2(mesh, order=4)
phi = GridFunction(fes, nested = True)
phi.Set(atan(sqrt((x0.5)**2 + (y0.5)**2)0.2))
print("phi: ", max(phi.vec), ", ", min(phi.vec))
for v in mesh.vertices:
print("grad ", cf_gphi(mesh(v.point[0], v.point[1])))
for i in range(3):
mesh.Refine()
phi.Update()
fes.Update()
cf_gphi = grad(phi)
for v in mesh.vertices:
print("grad ", cf_gphi(mesh(v.point[0], v.point[1])))
print("phi: ", max(phi.vec), ", ", min(phi.vec))
I would be very grateful if you could help me out here: DG is kind of new to me.
All the best,
Philipp