Laplacian of trace kernel for internal dofs TNT quadrilateral

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. To this end tabulation which also generates the derivatives is used.

python/demo/demo_tnt-elements.py

@@ -168,7 +168,7 @@
 
 
 def create_tnt_quad(degree):
-    assert degree > 1
+    assert degree > 0
     # Polyset
     ndofs = (degree + 1) ** 2 + 4
     npoly = (degree + 2) ** 2
@@ -218,10 +218,18 @@ def create_tnt_quad(degree):
         x[2].append(np.zeros([0, 2]))
         M[2].append(np.zeros([0, 1, 0, 1]))
     else:
-        pts, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 1)
-        poly = basix.tabulate_polynomials(
-            basix.PolynomialType.legendre, basix.CellType.quadrilateral, degree - 2, pts
+        pts, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree)
+        u = pts[:, 0]
+        v = pts[:, 1]
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree-2, 2, pts
         )
+        # this assumes the conventional [0 to 1][0 to 1] domain of the reference element, 
+        # and takes the Laplacian of (1-u)*u*(1-v)*v*pol_set[0], 
+        # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
+        poly = (pol_set[5]+pol_set[3])*(u-1)*u*(v-1)*v+ \
+                2*(pol_set[2]*(u-1)*u*(2*v-1)+pol_set[1]*(v-1)*v*(2*u-1)+ \
+                   pol_set[0]*((u-1)*u+(v-1)*v))
         face_ndofs = poly.shape[0]
         x[2].append(pts)
         mat = np.zeros((face_ndofs, 1, len(pts), 1))
@@ -300,8 +308,8 @@ def poisson_error(V: fem.FunctionSpace):
 tnt_degrees.append(2)
 tnt_ndofs.append(V.dofmap.index_map.size_global)
 tnt_errors.append(poisson_error(V))
-print(f"TNT degree 2 error: {tnt_errors[-1]}")
-for degree in range(2, 9):
+print(f"TNT degree 1 error: {tnt_errors[-1]}")
+for degree in range(1, 9):
     V = fem.functionspace(msh, create_tnt_quad(degree))
     tnt_degrees.append(degree + 1)
     tnt_ndofs.append(V.dofmap.index_map.size_global)