Ordering of local DoFs in LagrangeRefSpace

[sbadrian]

We are currently implementing the general p-order case for the Lagrange elements. For p=2 there is a (partial) implementation that implies a certain ordering of the local functions:

function (f::LagrangeRefSpace{T,2,3})(t) where T
    u,v,w, = barycentric(t)

    j = jacobian(t)
    p = t.patch

    #curl=(p[3]-p[2])/j),
     #   (value=v, curl=(p[1]-p[3])/j),
      #  (value=w, curl=(p[2]-p[1])/j)

    σ = sign(dot(normal(t), cross(p[1]-p[3],p[2]-p[3])))
     SVector(
        (value=u*(2*u-1), curl=σ*(p[3]-p[2])*(4u-1)/j),
        (value=v*(2*v-1), curl=σ*(p[1]-p[3])*(4v-1)/j),
        (value=w*(2*w-1), curl=σ*(p[2]-p[1])*(4w-1)/j),
        (value=4*v*w, curl=4*σ*(w*(p[1]-p[3])+v*(p[2]-p[1]))/j),
        (value=4*w*u, curl=4*σ*(w*(p[3]-p[2])+u*(p[2]-p[1]))/j),
        (value=4*u*v, curl=4*σ*(u*(p[1]-p[3])+v*(p[3]-p[2]))/j),
    )
end

It seems that the first three functions are on the vertices and then come the functions on the edges.

For the general p-case, it is more convenient to use the multiindex notation (i,j,k) (see Figure 4.3 in [1] P. P. Silvester and R. L. Ferrari, Finite elements for electrical engineers. Cambridge University Press, 1996). Silvester use different conversion strategy multi-index -> single-index, basically, by starting from index (3,0,0), then reducing the first index to 2 and setting the second index to the maximum (2,1,0) etc.

This conversion strategy would require to change the currently used ordering. I am happy to keep the existing one, but it is not clear to me what the general pth order case is in the current ordering system.