Laplacian of nullspace of trace in Q for internal dof functionals TNT quadrilateral

[jmv2009]

@mscroggs
According to https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py, we need to multiply by the primary bubble function, and take the Laplacian for the interior dofs. To this end tabulation which also generates the derivatives is used.
This will aid in generating 3D TNT elements in Basix.

[jhale]

Could you tighten up the PR title and description and move the extra details to a post further down?

[jmv2009]

Details:
The numerical results do not change at this point. This is because the exact definition of internal dofs does not really matter: As long as internal dof functionals are sufficiently independent from the point and edge dof functionals, the full space is spanned.
Please check contribution for form (including python conventions/requirements). Feel free to make your own version of this.
Note: Also need to check the quadrature degree to get to exactness, as it appears to be higher. [Correct: adjusted by +1]
First documented use of new tabulate_polynomial_set python function?

Reference implementation mathematica shape function tabulation (orders of vertices and edges are different):
k=3
xt={0,1,1,0,0}
yt={0,0,1,1,0}
space=Join[MonomialList[(1+x)^k(1+y)^k],{x,y}^(k+1),{x,y}^k x y]
tst=space.Inverse[Join[Table[space/.{x->xt[[i]],y->yt[[i]]},{i,4}],Flatten[Table[Integrate[space LegendreP[n,2s-1]Sqrt[2 n +1]/.{x->xt[i]+xt[[i+1]]s,y->yt[i]+yt[[i+1]]s},{s,0,1}],{i,4},{n,0,k-1}],1],Flatten[Table[Integrate[space Laplacian[s(s-1)t(t-1)LegendreP[n1,2s-1]Sqrt[2 n1+1]LegendreP[n2,2t-1]Sqrt[2 n2 +1],{s,t}]/.{x->s,y->t},{s,0,1},{t,0,1}],{n1,0,k-2},{n2,0,k-2}],1]]]
Flatten[Table[tst,{y,0,1,1/(k+1)},{x,0,1,1/(k+1)}],0]//N//MatrixForm

The only difference with symfem is that the dof functionals do not correspond to monomials but to normalized Legendre polynomials, as per the dof functional tabulation.