From 105f916cb67467ae83b1bfd7a4e951648b85499e Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 11 Dec 2024 20:35:09 +0100
Subject: [PATCH 01/19] 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 | 102 ++++++++++++++++++++-----------
 1 file changed, 67 insertions(+), 35 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 06eadf560b6..2cb3a45ad39 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -60,7 +60,9 @@
 #
 # We begin by defining a basis of the polynomial space spanned by the
 # TNT element, which is defined in terms of the orthogonal Legendre
-# polynomials on the cell. For a degree 1 element, the polynomial set
+# polynomials on the cell. For a degree 2 element (here, we use the 
+# superdegree rather than the conventional subdegree to allign with
+# the definition of other elements), the polynomial set
 # contains $1$, $y$, $y^2$, $x$, $xy$, $xy^2$, $x^2$, and $x^2y$, which
 # are the first 8 polynomials in the degree 2 set of polynomials on a
 # quadrilateral. We create an $8 \times 9$  matrix (number of dofs by
@@ -146,7 +148,7 @@
 # We now create the element. Using the Basix UFL interface, we can wrap
 # this element so that it can be used with FFCx/DOLFINx.
 
-tnt_degree1 = basix.ufl.custom_element(
+tnt_degree2 = basix.ufl.custom_element(
     basix.CellType.quadrilateral,
     [],
     wcoeffs,
@@ -168,23 +170,26 @@
 
 
 def create_tnt_quad(degree):
-    assert degree > 1
+    assert degree > 0
     # Polyset
-    ndofs = (degree + 1) ** 2 + 4
-    npoly = (degree + 2) ** 2
+    ndofs = 4*degree + max(degree-2,0)**2
+    npoly = (degree + 1) ** 2
 
     wcoeffs = np.zeros((ndofs, npoly))
 
     dof_n = 0
-    for i in range(degree + 1):
-        for j in range(degree + 1):
-            wcoeffs[dof_n, i * (degree + 2) + j] = 1
+    for i in range(degree):
+        for j in range(degree):
+            wcoeffs[dof_n, i * (degree + 1) + j] = 1
             dof_n += 1
 
-    for i, j in [(degree + 1, 1), (degree + 1, 0), (1, degree + 1), (0, degree + 1)]:
-        wcoeffs[dof_n, i * (degree + 2) + j] = 1
+    for i, j in [(degree, 1), (degree, 0), (0, degree)]:
+        wcoeffs[dof_n, i * (degree + 1) + j] = 1
         dof_n += 1
 
+    if degree > 1:
+        wcoeffs[dof_n, 2 * degree + 1] = 1
+
     # Interpolation
     geometry = basix.geometry(basix.CellType.quadrilateral)
     topology = basix.topology(basix.CellType.quadrilateral)
@@ -197,31 +202,44 @@ def create_tnt_quad(degree):
         M[0].append(np.array([[[[1.0]]]]))
 
     # Edges
-    pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree)
-    poly = basix.tabulate_polynomials(
-        basix.PolynomialType.legendre, basix.CellType.interval, degree - 1, pts
-    )
-    edge_ndofs = poly.shape[0]
-    for e in topology[1]:
-        v0 = geometry[e[0]]
-        v1 = geometry[e[1]]
-        edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
-        x[1].append(edge_pts)
-
-        mat = np.zeros((edge_ndofs, 1, len(pts), 1))
-        for i in range(edge_ndofs):
-            mat[i, 0, :, 0] = wts[:] * poly[i, :]
-        M[1].append(mat)
+    if degree < 2:
+        for _ in topology[1]:
+            x[1].append(np.zeros([0, 2]))
+            M[1].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
+        poly = basix.tabulate_polynomials(
+            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
+        )
+        edge_ndofs = poly.shape[0]
+        for e in topology[1]:
+            v0 = geometry[e[0]]
+            v1 = geometry[e[1]]
+            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            x[1].append(edge_pts)
+    
+            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
+            for i in range(edge_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[1].append(mat)
 
     # Interior
-    if degree == 1:
+    if degree < 3:
         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 - 2)
+        u = pts[:, 0]
+        v = pts[:, 1]
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 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))
@@ -239,8 +257,8 @@ def create_tnt_quad(degree):
         basix.MapType.identity,
         basix.SobolevSpace.H1,
         False,
+        max(degree - 1, 1),
         degree,
-        degree + 1,
         dtype=default_real_type,
     )
 
@@ -296,12 +314,12 @@ def poisson_error(V: fem.FunctionSpace):
 tnt_degrees = []
 tnt_errors = []
 
-V = fem.functionspace(msh, tnt_degree1)
+V = fem.functionspace(msh, tnt_degree2)
 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):
+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)
@@ -317,6 +335,17 @@ def poisson_error(V: fem.FunctionSpace):
     q_ndofs.append(V.dofmap.index_map.size_global)
     q_errors.append(poisson_error(V))
     print(f"Q degree {degree} error: {q_errors[-1]}")
+
+s_ndofs = []
+s_degrees = []
+s_errors = []
+for degree in range(1, 9):
+    V = fem.functionspace(msh, ("S", degree))
+    s_degrees.append(degree)
+    s_ndofs.append(V.dofmap.index_map.size_global)
+    s_errors.append(poisson_error(V))
+    print(f"S degree {degree} error: {s_errors[-1]}")
+
 # -
 
 # We now plot the data that we have obtained. First we plot the error
@@ -326,27 +355,30 @@ def poisson_error(V: fem.FunctionSpace):
 if MPI.COMM_WORLD.rank == 0:  # Only plot on one rank
     plt.plot(q_degrees, q_errors, "bo-")
     plt.plot(tnt_degrees, tnt_errors, "gs-")
+    plt.plot(s_degrees, s_errors, "rD-")
     plt.yscale("log")
     plt.xlabel("Polynomial degree")
     plt.ylabel("Error")
-    plt.legend(["Q", "TNT"])
+    plt.legend(["Q", "TNT", "S"])
     plt.savefig("demo_tnt-elements_degrees_vs_error.png")
     plt.clf()
 
 # ![](demo_tnt-elements_degrees_vs_error.png)
 #
 # A key advantage of TNT elements is that for a given degree, they span
-# a smaller polynomial space than Q elements. This can be observed in
+# a smaller polynomial space than Q tensorproduct elements. This can be observed in
 # the following diagram, where we plot the error against the square root
 # of the number of DOFs (providing a measure of cell size in 2D)
+# S serendipity element perform even better here.
 
 if MPI.COMM_WORLD.rank == 0:  # Only plot on one rank
     plt.plot(np.sqrt(q_ndofs), q_errors, "bo-")
     plt.plot(np.sqrt(tnt_ndofs), tnt_errors, "gs-")
+    plt.plot(np.sqrt(s_ndofs), s_errors, "rD-")
     plt.yscale("log")
     plt.xlabel("Square root of number of DOFs")
     plt.ylabel("Error")
-    plt.legend(["Q", "TNT"])
+    plt.legend(["Q", "TNT", "S"])
     plt.savefig("demo_tnt-elements_ndofs_vs_error.png")
     plt.clf()
 

From 0748fe14d4f212279b4ac6b38c06a94930d887cb Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Sun, 15 Dec 2024 15:50:27 +0000
Subject: [PATCH 02/19] Update demo_tnt-elements.py

Hexahedrons
---
 python/demo/demo_tnt-elements.py | 125 +++++++++++++++++++++++++++++++
 1 file changed, 125 insertions(+)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 2cb3a45ad39..f81ef2fac2c 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -383,3 +383,128 @@ def poisson_error(V: fem.FunctionSpace):
     plt.clf()
 
 # ![](demo_tnt-elements_ndofs_vs_error.png)
+
+# extension to hexahedrons is trivial:
+
+def create_tnt_hex(degree):
+    assert degree > 0
+    # Polyset
+    ndofs = 12 * degree - 4 + (k + 4) * max(degree-2, 0) ** 2
+    npoly = (degree + 1) ** 3
+
+    wcoeffs = np.zeros((ndofs, npoly))
+
+    dof_n = 0
+    for i in range(degree):
+        for j in range(degree):
+            for k in range(degree):
+                wcoeffs[dof_n, (i * (degree + 2) + j) * (degree + 1) + k ] = 1
+                dof_n += 1
+
+    for i, j, k in [(degree, 0, 0), (0, degree, 0), (0, 0, degree), (1, 1, degree), (1, 0, degree), (0, 1, degree), (0, degree, 1)]:
+        wcoeffs[dof_n, (i * (degree + 2) + j) * (degree + 1) + k] = 1
+        dof_n += 1
+
+    if degree > 1:
+        for i, j in [(1, degree, 1), (degree, 1, 1), (degree, 0, 1), (degree, 0, 1), (0, 1, degree)]:
+            wcoeffs[dof_n, (i * (degree + 2) + j) * (degree + 1) + k] = 1
+            dof_n += 1
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.hexahedral)
+    topology = basix.topology(basix.CellType.hexahedral)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
+
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.array(geometry[v]))
+        M[0].append(np.array([[[[1.0]]]]))
+
+    # Edges
+    if degree < 2:
+        for _ in topology[1]:
+            x[1].append(np.zeros([0, 2]))
+            M[1].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
+        poly = basix.tabulate_polynomials(
+            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
+        )
+        edge_ndofs = poly.shape[0]
+        for e in topology[1]:
+            v0 = geometry[e[0]]
+            v1 = geometry[e[1]]
+            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            x[1].append(edge_pts)
+    
+            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
+            for i in range(edge_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[1].append(mat)
+
+    # Interior
+    if degree < 3:
+        for _ in topology[2]:
+            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 - 2)
+        u = pts[:, 0]
+        v = pts[:, 1]
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 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))
+        for i in range(face_ndofs):
+            mat[i, 0, :, 0] = wts[:] * poly[i, :]
+        M[2].append(mat)
+
+    if degree < 3:
+        x[3].append(np.zeros([0, 2]))
+        M[3].append(np.zeros([0, 1, 0, 1]))
+
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.hexahedral, 2 * degree - 2)
+        u = pts[:, 0]
+        v = pts[:, 1]
+        w = pts[:  2]
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.hexahedral, basix.PolynomialType.legendre, degree - 3, 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*(1-w)*w*pol_set[0], 
+        # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
+        poly = (pol_set[9]+pol_set[7]+pol_set[4])*(u-1)*u*(v-1)*v*(w-1)*w+ \
+                2*(pol_set[3]*(u-1)*u*(v-1)*v*(2*w-1)+pol_set[2]*(u-1)*u*(w-1)*w*(2*v-1)+pol_set[1]*(v-1)*v*(w-1)*w*(2*u-1)+ \
+                   pol_set[0]*((u-1)*u*(v-1)*v+(u-1)*u*(w-1)*w+(v-1)*v))
+        face_ndofs = poly.shape[0]
+        x[3].append(pts)
+        mat = np.zeros((face_ndofs, 1, len(pts), 1))
+        for i in range(face_ndofs):
+            mat[i, 0, :, 0] = wts[:] * poly[i, :]
+        M[3].append(mat)
+        
+    return basix.ufl.custom_element(
+        basix.CellType.hexahedral,
+        [],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.identity,
+        basix.SobolevSpace.H1,
+        False,
+        max(degree - 1,1),
+        degree,
+        dtype=default_real_type,
+    )

From c3c3634525a9e7385f79925ee6915e58375c11d1 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Sun, 15 Dec 2024 19:56:17 +0000
Subject: [PATCH 03/19] Update demo_tnt-elements.py

Fix hex
---
 python/demo/demo_tnt-elements.py | 82 ++++++++++++++++----------------
 1 file changed, 41 insertions(+), 41 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index f81ef2fac2c..d523b980127 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -386,33 +386,33 @@ def poisson_error(V: fem.FunctionSpace):
 
 # extension to hexahedrons is trivial:
 
+
 def create_tnt_hex(degree):
     assert degree > 0
     # Polyset
-    ndofs = 12 * degree - 4 + (k + 4) * max(degree-2, 0) ** 2
+    ndofs = 12 * degree - 4 + (degree + 4) * max(degree-2, 0) ** 2
     npoly = (degree + 1) ** 3
-
     wcoeffs = np.zeros((ndofs, npoly))
 
     dof_n = 0
     for i in range(degree):
         for j in range(degree):
             for k in range(degree):
-                wcoeffs[dof_n, (i * (degree + 2) + j) * (degree + 1) + k ] = 1
+                wcoeffs[dof_n, (i * (degree + 1) + j) * (degree + 1) + k ] = 1
                 dof_n += 1
 
-    for i, j, k in [(degree, 0, 0), (0, degree, 0), (0, 0, degree), (1, 1, degree), (1, 0, degree), (0, 1, degree), (0, degree, 1)]:
-        wcoeffs[dof_n, (i * (degree + 2) + j) * (degree + 1) + k] = 1
+    for i, j, k in [(degree, 0, 0), (0, degree, 0), (0, 0, degree), (1, 1, degree), (1, 0, degree), (0, 1, degree), (degree, 1, 0)]:
+        wcoeffs[dof_n, (i * (degree + 1) + j) * (degree + 1) + k] = 1
         dof_n += 1
 
     if degree > 1:
-        for i, j in [(1, degree, 1), (degree, 1, 1), (degree, 0, 1), (degree, 0, 1), (0, 1, degree)]:
-            wcoeffs[dof_n, (i * (degree + 2) + j) * (degree + 1) + k] = 1
+        for i, j, k  in [(1, degree, 1), (degree, 1, 1), (degree, 0, 1), (0, degree, 1), (1, degree, 0)]:
+            wcoeffs[dof_n, (i * (degree + 1) + j) * (degree + 1) + k] = 1
             dof_n += 1
 
     # Interpolation
-    geometry = basix.geometry(basix.CellType.hexahedral)
-    topology = basix.topology(basix.CellType.hexahedral)
+    geometry = basix.geometry(basix.CellType.hexahedron)
+    topology = basix.topology(basix.CellType.hexahedron)
     x = [[], [], [], []]
     M = [[], [], [], []]
 
@@ -424,7 +424,7 @@ def create_tnt_hex(degree):
     # Edges
     if degree < 2:
         for _ in topology[1]:
-            x[1].append(np.zeros([0, 2]))
+            x[1].append(np.zeros([0, 3]))
             M[1].append(np.zeros([0, 1, 0, 1]))
     else:
         pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
@@ -443,59 +443,59 @@ def create_tnt_hex(degree):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
             M[1].append(mat)
 
-    # Interior
+    # Faces
     if degree < 3:
         for _ in topology[2]:
-            x[2].append(np.zeros([0, 2]))
+            x[2].append(np.zeros([0, 3]))
             M[2].append(np.zeros([0, 1, 0, 1]))
-
     else:
-        pts, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2)
-        u = pts[:, 0]
-        v = pts[:, 1]
-        pol_set = basix.polynomials.tabulate_polynomial_set(
-            basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, pts
-        )
+            ptsr, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2)
+            pol_set = basix.polynomials.tabulate_polynomial_set(
+                basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, ptsr
+                )
         # 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))
-        for i in range(face_ndofs):
-            mat[i, 0, :, 0] = wts[:] * poly[i, :]
-        M[2].append(mat)
+            # 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
+            u = ptsr[:, 0]
+            v = ptsr[:, 1]
+            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]
+            for f in topology[2]:
+                x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
+                mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
+                for i in range(face_ndofs):
+                    mat[i, 0, :, 0] = wts[:] * poly[i, :]
+                M[2].append(mat)
 
+    # Interior
     if degree < 3:
-        x[3].append(np.zeros([0, 2]))
+        x[3].append(np.zeros([0, 3]))
         M[3].append(np.zeros([0, 1, 0, 1]))
-
     else:
-        pts, wts = basix.make_quadrature(basix.CellType.hexahedral, 2 * degree - 2)
-        u = pts[:, 0]
-        v = pts[:, 1]
-        w = pts[:  2]
+        pts, wts = basix.make_quadrature(basix.CellType.hexahedron, 2 * degree - 2)
         pol_set = basix.polynomials.tabulate_polynomial_set(
-            basix.CellType.hexahedral, basix.PolynomialType.legendre, degree - 3, 2, pts
+            basix.CellType.hexahedron, basix.PolynomialType.legendre, degree - 3, 2, pts
         )
+        u = pts[:, 0]
+        v = pts[:, 1]
+        w = pts[:, 2]
         # 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*(1-w)*w*pol_set[0], 
         # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
         poly = (pol_set[9]+pol_set[7]+pol_set[4])*(u-1)*u*(v-1)*v*(w-1)*w+ \
                 2*(pol_set[3]*(u-1)*u*(v-1)*v*(2*w-1)+pol_set[2]*(u-1)*u*(w-1)*w*(2*v-1)+pol_set[1]*(v-1)*v*(w-1)*w*(2*u-1)+ \
                    pol_set[0]*((u-1)*u*(v-1)*v+(u-1)*u*(w-1)*w+(v-1)*v))
-        face_ndofs = poly.shape[0]
+        vol_ndofs = poly.shape[0]
         x[3].append(pts)
-        mat = np.zeros((face_ndofs, 1, len(pts), 1))
-        for i in range(face_ndofs):
+        mat = np.zeros((vol_ndofs, 1, len(pts), 1))
+        for i in range(vol_ndofs):
             mat[i, 0, :, 0] = wts[:] * poly[i, :]
         M[3].append(mat)
         
     return basix.ufl.custom_element(
-        basix.CellType.hexahedral,
+        basix.CellType.hexahedron,
         [],
         wcoeffs,
         x,

From 39ab52547f40476819e59ec6a3ba0dcb1421fe09 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Sun, 15 Dec 2024 20:01:55 +0000
Subject: [PATCH 04/19] Update demo_tnt-elements.py

fix hex
---
 python/demo/demo_tnt-elements.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index d523b980127..b0af6d303aa 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -486,7 +486,7 @@ def create_tnt_hex(degree):
         # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
         poly = (pol_set[9]+pol_set[7]+pol_set[4])*(u-1)*u*(v-1)*v*(w-1)*w+ \
                 2*(pol_set[3]*(u-1)*u*(v-1)*v*(2*w-1)+pol_set[2]*(u-1)*u*(w-1)*w*(2*v-1)+pol_set[1]*(v-1)*v*(w-1)*w*(2*u-1)+ \
-                   pol_set[0]*((u-1)*u*(v-1)*v+(u-1)*u*(w-1)*w+(v-1)*v))
+                   pol_set[0]*((u-1)*u*(v-1)*v+(u-1)*u*(w-1)*w+(v-1)*v*(w-1)*w))
         vol_ndofs = poly.shape[0]
         x[3].append(pts)
         mat = np.zeros((vol_ndofs, 1, len(pts), 1))

From c48a1eec6e96780f2a6425a925641b33186ddc1e Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Mon, 16 Dec 2024 21:32:41 +0000
Subject: [PATCH 05/19] Provisional TNT Prism 0Form

---
 python/demo/demo_tnt-elements.py | 179 ++++++++++++++++++++++++++++---
 1 file changed, 162 insertions(+), 17 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index b0af6d303aa..709caed2dee 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -449,25 +449,25 @@ def create_tnt_hex(degree):
             x[2].append(np.zeros([0, 3]))
             M[2].append(np.zeros([0, 1, 0, 1]))
     else:
-            ptsr, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2)
-            pol_set = basix.polynomials.tabulate_polynomial_set(
-                basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, ptsr
+        ptsr, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2)
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, ptsr
                 )
         # 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
-            u = ptsr[:, 0]
-            v = ptsr[:, 1]
-            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]
-            for f in topology[2]:
-                x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
-                mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
-                for i in range(face_ndofs):
-                    mat[i, 0, :, 0] = wts[:] * poly[i, :]
-                M[2].append(mat)
+        # 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
+        u = ptsr[:, 0]
+        v = ptsr[:, 1]
+        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]
+        for f in topology[2]:
+            x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
+            mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
+            for i in range(face_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[2].append(mat)
 
     # Interior
     if degree < 3:
@@ -508,3 +508,148 @@ def create_tnt_hex(degree):
         degree,
         dtype=default_real_type,
     )
+
+def create_tnt_prism(degree):
+    assert degree > 0
+    # Polyset
+    ndofs = 9 * degree - 3 + (round((degree + 5) * (degree - 2) / 2) + 1) * max(degree - 2, 0)
+    npoly = round((degree + 1) * (degree + 2) / 2) * (degree + 1)
+
+    wcoeffs = np.zeros((ndofs, npoly))
+
+    dof_n = 0
+    for i in range(round((degree + 1) * degree / 2)):
+        for j in range(degree):
+            wcoeffs[dof_n, i * (degree + 1) + j] = 1
+            dof_n += 1
+
+    for i in range(3):
+        if degree > 1 or i == 0:
+            wcoeffs[dof_n, i * (degree + 1) + degree] = 1
+            dof_n += 1
+    
+    for i in range(degree+1):
+        for j in range(2):
+            wcoeffs[dof_n, (i + round((degree + 1) * degree / 2)) * (degree + 1) + j] = 1
+            dof_n += 1
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.prism)
+    topology = basix.topology(basix.CellType.prism)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
+
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.array(geometry[v]))
+        M[0].append(np.array([[[[1.0]]]]))
+
+    # Edges
+    if degree < 2:
+        for _ in topology[1]:
+            x[1].append(np.zeros([0, 3]))
+            M[1].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
+        poly = basix.tabulate_polynomials(
+            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
+        )
+        edge_ndofs = poly.shape[0]
+        for e in topology[1]:
+            v0 = geometry[e[0]]
+            v1 = geometry[e[1]]
+            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            x[1].append(edge_pts)
+    
+            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
+            for i in range(edge_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[1].append(mat)
+
+    # Faces
+    if degree < 3:
+        for _ in topology[2]:
+            x[2].append(np.zeros([0, 3]))
+            M[2].append(np.zeros([0, 1, 0, 1]))
+    else:
+        ptsr_t, wts_t = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2) # or triangle for face 0 or 5
+        pol_set_t = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 2, ptsr_t
+        )
+        ptsr_q, wts_q = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2) # or triangle for face 0 or 5
+        pol_set_q = basix.polynomials.tabulate_polynomial_set(
+             basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, ptsr_q
+            )
+        # 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
+        u = ptsr_t[:, 0]
+        v = ptsr_t[:, 1]
+        poly_t = (pol_set_t[5]+pol_set_t[3])*(u-1)*u*(v-1)*v+ \
+                2*(pol_set_t[2]*(u-1)*u*(2*v-1)+pol_set_t[1]*(v-1)*v*(2*u-1)+ \
+                   pol_set_t[0]*((u-1)*u+(v-1)*v))
+        u = ptsr_q[:, 0]
+        v = ptsr_q[:, 1]
+        poly_q = (pol_set_q[5]+pol_set_q[3])*(u-1)*u*(v-1)*v+ \
+                2*(pol_set_q[2]*(u-1)*u*(2*v-1)+pol_set_q[1]*(v-1)*v*(2*u-1)+ \
+                   pol_set_q[0]*((u-1)*u+(v-1)*v))
+        for f in topology[2]:
+            if geometry[f].shape[0] == 3:
+                poly = poly_t
+                wts = wts_t
+                ptsr = ptsr_t
+            else:
+                poly = poly_q
+                wts = wts_q
+                ptsr = ptsr_q
+            face_ndofs = poly.shape[0]
+            x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
+            mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
+            for i in range(face_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[2].append(mat)
+
+    # Interior
+    if degree < 3:
+        x[3].append(np.zeros([0, 3]))
+        M[3].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.prism, 2 * degree - 2)
+        #The dimension of the left over space tells us we should reduce the xy space by 1, so we are making the appopriate selection
+        sel=[]
+        for i in range(round((degree - 2) * (degree - 3)/ 2)):
+            for j in range(degree - 2):
+                sel.append(i * (degree  - 2) + j)
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.prism, basix.PolynomialType.legendre, degree - 3, 2, pts)[:,sel,:]
+        u = pts[:, 0]
+        v = pts[:, 1]
+        w = pts[:, 2]    
+        # this assumes the conventional [0 to 1][0 to 1] domain of the reference element, 
+        # and takes the Laplacian of (u+v-1)*u*v*(w-1)*w*pol_set[0], 
+        # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
+        poly = (pol_set[9]+pol_set[7]+pol_set[4])*(u+v-1)*u*v*(w-1)*w+ \
+                2*(pol_set[3]*(u+v-1)*u*v*(2*w-1)+pol_set[2]*u*(w-1)*w*(2*v+u-1)+pol_set[1]*v*(w-1)*w*(2*u+v-1)+ \
+                   pol_set[0]*((u+v-1)*u*v+u*(w-1)*w+v*(w-1)*w))
+        poly=pol_set[0]
+        vol_ndofs = poly.shape[0]
+        x[3].append(pts)
+        mat = np.zeros((vol_ndofs, 1, len(pts), 1))
+        for i in range(vol_ndofs):
+            mat[i, 0, :, 0] = wts[:] * poly[i, :]
+        M[3].append(mat)
+    
+    return basix.ufl.custom_element(
+        basix.CellType.prism,
+        [],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.identity,
+        basix.SobolevSpace.H1,
+        False,
+        max(degree - 1,1),
+        degree,
+        dtype=default_real_type,
+    )

From 170660ca1773198f81644442edaeddbda669d263 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Mon, 16 Dec 2024 22:45:29 +0100
Subject: [PATCH 06/19] Update demo_tnt-elements.py

---
 python/demo/demo_tnt-elements.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 709caed2dee..d04242c62f9 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -610,12 +610,13 @@ def create_tnt_prism(degree):
             M[2].append(mat)
 
     # Interior
-    if degree < 3:
+    if degree < 4:
         x[3].append(np.zeros([0, 3]))
         M[3].append(np.zeros([0, 1, 0, 1]))
     else:
         pts, wts = basix.make_quadrature(basix.CellType.prism, 2 * degree - 2)
         #The dimension of the left over space tells us we should reduce the xy space by 1, so we are making the appopriate selection
+        #This corresponds to tetrahedrons not needing an internal dof at degree 3
         sel=[]
         for i in range(round((degree - 2) * (degree - 3)/ 2)):
             for j in range(degree - 2):

From 57022e7095d7947e93ac6685970ae5d3486bc13a Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Tue, 17 Dec 2024 16:34:24 +0000
Subject: [PATCH 07/19] Update demo_tnt-elements.py

Fix prism
---
 python/demo/demo_tnt-elements.py | 32 +++++++++++++++++++-------------
 1 file changed, 19 insertions(+), 13 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index d04242c62f9..8b1aa5c3a9b 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -523,16 +523,20 @@ def create_tnt_prism(degree):
             wcoeffs[dof_n, i * (degree + 1) + j] = 1
             dof_n += 1
 
-    for i in range(3):
-        if degree > 1 or i == 0:
-            wcoeffs[dof_n, i * (degree + 1) + degree] = 1
-            dof_n += 1
-    
     for i in range(degree+1):
         for j in range(2):
             wcoeffs[dof_n, (i + round((degree + 1) * degree / 2)) * (degree + 1) + j] = 1
             dof_n += 1
 
+    wcoeffs[dof_n, degree] = 1
+    dof_n += 1
+
+    if degree > 1:
+        for i in range(1,3):
+            wcoeffs[dof_n, i * (degree + 1) + degree] = 1
+            dof_n += 1
+
+
     # Interpolation
     geometry = basix.geometry(basix.CellType.prism)
     topology = basix.topology(basix.CellType.prism)
@@ -563,7 +567,7 @@ def create_tnt_prism(degree):
     
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
-                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+                mat[i, 0, :, 0] = wts[:] * poly[i, :] * np.linalg.norm(v1 - v0)
             M[1].append(mat)
 
     # Faces
@@ -581,13 +585,13 @@ def create_tnt_prism(degree):
              basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, ptsr_q
             )
         # 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], 
+        # and takes the Laplacian of (1-u)*u*(1-v)*v*pol_set[0] or (u+v-1)*u*v*pol_set[0], , 
         # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
         u = ptsr_t[:, 0]
         v = ptsr_t[:, 1]
-        poly_t = (pol_set_t[5]+pol_set_t[3])*(u-1)*u*(v-1)*v+ \
-                2*(pol_set_t[2]*(u-1)*u*(2*v-1)+pol_set_t[1]*(v-1)*v*(2*u-1)+ \
-                   pol_set_t[0]*((u-1)*u+(v-1)*v))
+        poly_t = (pol_set_t[5]+pol_set_t[3])*(u+v-1)*u*v+ \
+                2*(pol_set_t[2]*u*(2*v+u-1)+pol_set_t[1]*v*(2*u+v-1)+ \
+                   pol_set_t[0]*(u+v))
         u = ptsr_q[:, 0]
         v = ptsr_q[:, 1]
         poly_q = (pol_set_q[5]+pol_set_q[3])*(u-1)*u*(v-1)*v+ \
@@ -606,17 +610,16 @@ def create_tnt_prism(degree):
             x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
             mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
             for i in range(face_ndofs):
-                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+                mat[i, 0, :, 0] = wts[:] * poly[i, :] * np.linalg.norm(np.cross(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))
             M[2].append(mat)
 
     # Interior
-    if degree < 4:
+    if degree < 3:
         x[3].append(np.zeros([0, 3]))
         M[3].append(np.zeros([0, 1, 0, 1]))
     else:
         pts, wts = basix.make_quadrature(basix.CellType.prism, 2 * degree - 2)
         #The dimension of the left over space tells us we should reduce the xy space by 1, so we are making the appopriate selection
-        #This corresponds to tetrahedrons not needing an internal dof at degree 3
         sel=[]
         for i in range(round((degree - 2) * (degree - 3)/ 2)):
             for j in range(degree - 2):
@@ -639,6 +642,8 @@ def create_tnt_prism(degree):
         for i in range(vol_ndofs):
             mat[i, 0, :, 0] = wts[:] * poly[i, :]
         M[3].append(mat)
+
+#    print("x: ", x)
     
     return basix.ufl.custom_element(
         basix.CellType.prism,
@@ -654,3 +659,4 @@ def create_tnt_prism(degree):
         degree,
         dtype=default_real_type,
     )
+

From f5639bd51c60eec00f692f1dab7cc01a31ad7a14 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Tue, 17 Dec 2024 19:36:21 +0000
Subject: [PATCH 08/19] Fix things, add triangle and interval

---
 python/demo/demo_tnt-elements.py | 153 +++++++++++++++++++++++++++++--
 1 file changed, 146 insertions(+), 7 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 8b1aa5c3a9b..422d43635fe 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -315,13 +315,13 @@ def poisson_error(V: fem.FunctionSpace):
 tnt_errors = []
 
 V = fem.functionspace(msh, tnt_degree2)
-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]}")
+#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(1, 9):
     V = fem.functionspace(msh, create_tnt_quad(degree))
-    tnt_degrees.append(degree + 1)
+    tnt_degrees.append(degree)
     tnt_ndofs.append(V.dofmap.index_map.size_global)
     tnt_errors.append(poisson_error(V))
     print(f"TNT degree {degree} error: {tnt_errors[-1]}")
@@ -567,7 +567,7 @@ def create_tnt_prism(degree):
     
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
-                mat[i, 0, :, 0] = wts[:] * poly[i, :] * np.linalg.norm(v1 - v0)
+                mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(v1 - v0)
             M[1].append(mat)
 
     # Faces
@@ -610,7 +610,7 @@ def create_tnt_prism(degree):
             x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
             mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
             for i in range(face_ndofs):
-                mat[i, 0, :, 0] = wts[:] * poly[i, :] * np.linalg.norm(np.cross(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))
+                mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(np.cross(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))
             M[2].append(mat)
 
     # Interior
@@ -660,3 +660,142 @@ def create_tnt_prism(degree):
         dtype=default_real_type,
     )
 
+# Here is the matching triangle and interval elements. Their boundary dofs are moments, not point values.
+
+def create_tnt_triangle(degree):
+    assert degree > 0
+    # Polyset
+    ndofs = round((degree+1)*(degree+2)/2)
+    npoly = round((degree+1)*(degree+2)/2)
+    
+    wcoeffs = np.eye(ndofs)
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.triangle)
+    topology = basix.topology(basix.CellType.triangle)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
+
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.array(geometry[v]))
+        M[0].append(np.array([[[[1.0]]]]))
+
+    # Edges
+    if degree < 2:
+        for _ in topology[1]:
+            x[1].append(np.zeros([0, 2]))
+            M[1].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
+        poly = basix.tabulate_polynomials(
+            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
+        )
+        edge_ndofs = poly.shape[0]
+        for e in topology[1]:
+            v0 = geometry[e[0]]
+            v1 = geometry[e[1]]
+            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            x[1].append(edge_pts)
+    
+            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
+            for i in range(edge_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(v1 - v0)
+            M[1].append(mat)
+
+    # Faces
+    if degree < 3:
+        for _ in topology[2]:
+            x[2].append(np.zeros([0, 2]))
+            M[2].append(np.zeros([0, 1, 0, 1]))
+    else:
+        ptsr, wts = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2) # or triangle for face 0 or 5
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 2, ptsr
+        )
+        # this assumes the conventional [0 to 1][0 to 1] domain of the reference element, 
+        # and takes the Laplacian of (u+v-1)*u*v*pol_set[0], , 
+        # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
+        u = ptsr[:, 0]
+        v = ptsr[:, 1]
+        poly = (pol_set[5]+pol_set[3])*(u+v-1)*u*v+ \
+                2*(pol_set[2]*u*(2*v+u-1)+pol_set[1]*v*(2*u+v-1)+ \
+                   pol_set[0]*(u+v))
+
+        face_ndofs = poly.shape[0]
+        x[2].append(ptsr)
+        mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
+        for i in range(face_ndofs):
+            mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(np.cross(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))
+        M[2].append(mat)
+    
+    return basix.ufl.custom_element(
+        basix.CellType.triangle,
+        [],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.identity,
+        basix.SobolevSpace.H1,
+        False,
+        degree,
+        degree,
+        dtype=default_real_type,
+    )
+
+def create_tnt_interval(degree):
+    assert degree > 0
+    # Polyset
+    ndofs = degree + 1
+    npoly = degree + 1
+    
+    wcoeffs = np.eye(ndofs)
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.interval)
+    topology = basix.topology(basix.CellType.interval)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
+
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.array(geometry[v]))
+        M[0].append(np.array([[[[1.0]]]]))
+
+    # Edges
+    if degree < 2:
+        for _ in topology[1]:
+            x[1].append(np.zeros([0, 1]))
+            M[1].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
+        poly = basix.tabulate_polynomials(
+            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
+        )
+        edge_ndofs = poly.shape[0]
+        for e in topology[1]:
+            v0 = geometry[e[0]]
+            v1 = geometry[e[1]]
+            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            x[1].append(edge_pts)
+    
+            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
+            for i in range(edge_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[1].append(mat)
+    
+    return basix.ufl.custom_element(
+        basix.CellType.interval,
+        [],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.identity,
+        basix.SobolevSpace.H1,
+        False,
+        degree,
+        degree,
+        dtype=default_real_type,
+    )

From eef0d8fe6168c5bda7585b979035b7f496fffcef Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 18 Dec 2024 16:08:24 +0000
Subject: [PATCH 09/19] Update demo_tnt-elements.py

Tetrahedron, misc.
---
 python/demo/demo_tnt-elements.py | 162 +++++++++++++++++++++++++++----
 1 file changed, 144 insertions(+), 18 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 422d43635fe..1af034c0880 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -61,7 +61,7 @@
 # We begin by defining a basis of the polynomial space spanned by the
 # TNT element, which is defined in terms of the orthogonal Legendre
 # polynomials on the cell. For a degree 2 element (here, we use the 
-# superdegree rather than the conventional subdegree to allign with
+# superdegree rather than the conventional subdegree to align with
 # the definition of other elements), the polynomial set
 # contains $1$, $y$, $y^2$, $x$, $xy$, $xy^2$, $x^2$, and $x^2y$, which
 # are the first 8 polynomials in the degree 2 set of polynomials on a
@@ -263,10 +263,10 @@ def create_tnt_quad(degree):
     )
 
 
-# ## Comparing TNT elements and Q elements
+# ## Comparing TNT elements, Q and S elements
 #
 # We now use the code above to compare TNT elements and
-# [Q](https://defelement.com/elements/lagrange.html) elements on
+# [Q and S](https://defelement.com/elements/lagrange.html) elements on
 # quadrilaterals. The following function takes a DOLFINx function space
 # as input, and solves a Poisson problem and returns the $L_2$ error of
 # the solution.
@@ -303,9 +303,29 @@ def poisson_error(V: fem.FunctionSpace):
 
 
 # We create a mesh, then solve the Poisson problem using our TNT
-# elements of degree 1 to 8. We then do the same with Q elements of
-# degree 1 to 9. For the TNT elements, we store a number 1 larger than
-# the degree as this is the highest degree polynomial in the space.
+# elements of degree 1 to 9. We then do the same with Q elements of
+# degree 1 to 9.
+
+# 0 Form S elements have the following polynomial power pattern in 2D:
+# 00000
+# 0000
+# 000
+# 00
+# 0
+
+# 0 Form TNT element have the following polynomial power pattern in 2D:
+# 00000
+# 00000
+# 0000
+# 0000
+# 00
+
+# 0 elements have the following polynomial power pattern in 2D:
+# 00000
+# 00000
+# 00000
+# 00000
+# 00000
 
 # +
 msh = mesh.create_unit_square(MPI.COMM_WORLD, 15, 15, mesh.CellType.quadrilateral)
@@ -443,18 +463,20 @@ def create_tnt_hex(degree):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
             M[1].append(mat)
 
-    # Faces
+
     if degree < 3:
         for _ in topology[2]:
             x[2].append(np.zeros([0, 3]))
             M[2].append(np.zeros([0, 1, 0, 1]))
+        x[3].append(np.zeros([0, 3]))
+        M[3].append(np.zeros([0, 1, 0, 1]))
     else:
+        # Faces
         ptsr, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2)
         pol_set = basix.polynomials.tabulate_polynomial_set(
             basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, ptsr
                 )
-        # 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], 
+        # This 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
         u = ptsr[:, 0]
         v = ptsr[:, 1]
@@ -469,11 +491,8 @@ def create_tnt_hex(degree):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
             M[2].append(mat)
 
-    # Interior
-    if degree < 3:
-        x[3].append(np.zeros([0, 3]))
-        M[3].append(np.zeros([0, 1, 0, 1]))
-    else:
+        # Interior
+
         pts, wts = basix.make_quadrature(basix.CellType.hexahedron, 2 * degree - 2)
         pol_set = basix.polynomials.tabulate_polynomial_set(
             basix.CellType.hexahedron, basix.PolynomialType.legendre, degree - 3, 2, pts
@@ -481,8 +500,7 @@ def create_tnt_hex(degree):
         u = pts[:, 0]
         v = pts[:, 1]
         w = pts[:, 2]
-        # 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*(1-w)*w*pol_set[0], 
+        # this takes as the Laplacian of (1-u)*u*(1-v)*v*(1-w)*w*pol_set[0], 
         # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
         poly = (pol_set[9]+pol_set[7]+pol_set[4])*(u-1)*u*(v-1)*v*(w-1)*w+ \
                 2*(pol_set[3]*(u-1)*u*(v-1)*v*(2*w-1)+pol_set[2]*(u-1)*u*(w-1)*w*(2*v-1)+pol_set[1]*(v-1)*v*(w-1)*w*(2*u-1)+ \
@@ -614,7 +632,7 @@ def create_tnt_prism(degree):
             M[2].append(mat)
 
     # Interior
-    if degree < 3:
+    if degree < 4:
         x[3].append(np.zeros([0, 3]))
         M[3].append(np.zeros([0, 1, 0, 1]))
     else:
@@ -660,8 +678,116 @@ def create_tnt_prism(degree):
         dtype=default_real_type,
     )
 
-# Here is the matching triangle and interval elements. Their boundary dofs are moments, not point values.
+# Here is the matching tetrahedron, triangle and interval elements. Their boundary dofs are moments, not point values, except at the vertices.
+
+def create_tnt_tetrahedron(degree):
+    assert degree > 0
+    # Polyset
+    ndofs = round((degree+3)*(degree+1)*(degree+2)/6)
+    npoly = round((degree+3)*(degree+1)*(degree+2)/6)
+    
+    wcoeffs = np.eye(ndofs)
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.tetrahedron)
+    topology = basix.topology(basix.CellType.tetrahedron)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
+
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.array(geometry[v]))
+        M[0].append(np.array([[[[1.0]]]]))
+
+    # Edges
+    if degree < 2:
+        for _ in topology[1]:
+            x[1].append(np.zeros([0, 3]))
+            M[1].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
+        poly = basix.tabulate_polynomials(
+            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
+        )
+        edge_ndofs = poly.shape[0]
+        for e in topology[1]:
+            v0 = geometry[e[0]]
+            v1 = geometry[e[1]]
+            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            x[1].append(edge_pts)
+    
+            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
+            for i in range(edge_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[1].append(mat)
+
+    # Faces
+    if degree < 3:
+        for _ in topology[2]:
+            x[2].append(np.zeros([0, 3]))
+            M[2].append(np.zeros([0, 1, 0, 1]))
+    else:
+        ptsr, wts = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2)
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 2, ptsr
+        )
+        # this assumes the conventional [0 to 1][0 to 1] domain of the reference element, 
+        # and takes the Laplacian of (u+v-1)*u*v*pol_set[0], , 
+        # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
+        u = ptsr[:, 0]
+        v = ptsr[:, 1]
+        poly = (pol_set[5]+pol_set[3])*(u+v-1)*u*v+ \
+                2*(pol_set[2]*u*(2*v+u-1)+pol_set[1]*v*(2*u+v-1)+ \
+                   pol_set[0]*(u+v))
+
+        for f in topology[2]:
+            face_ndofs = poly.shape[0]
+            x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
+            mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
+            for i in range(face_ndofs):
+                mat[i, 0, :, 0] = wts[:] * poly[i, :]
+            M[2].append(mat)
+
+    # Interior
+    if degree < 4:
+        x[3].append(np.zeros([0, 3]))
+        M[3].append(np.zeros([0, 1, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.tetrahedron, 2 * degree - 2)
+        pol_set = basix.polynomials.tabulate_polynomial_set(
+            basix.CellType.tetrahedron, basix.PolynomialType.legendre, degree - 4, 2, pts)
+        u = pts[:, 0]
+        v = pts[:, 1]
+        w = pts[:, 2]    
+        # this assumes the conventional [0 to 1][0 to 1] domain of the reference element, 
+        # and takes the Laplacian of u*v*w*(u+v+w-1)*pol_set[0], 
+        # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
+        poly = (pol_set[9]+pol_set[7]+pol_set[4])*u*v*w*(u+v+w-1)+ \
+                2*(pol_set[3]*(2*w+u+v-1)*u*v+pol_set[2]*u*w*(2*v+u+w-1)+pol_set[1]*v*w*(2*u+v+w-1)+ \
+                   pol_set[0]*(u*v+u*w+v*w))
+        poly=pol_set[0]
+        vol_ndofs = poly.shape[0]
+        x[3].append(pts)
+        mat = np.zeros((vol_ndofs, 1, len(pts), 1))
+        for i in range(vol_ndofs):
+            mat[i, 0, :, 0] = wts[:] * poly[i, :]
+        M[3].append(mat)
 
+    return basix.ufl.custom_element(
+        basix.CellType.tetrahedron,
+        [],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.identity,
+        basix.SobolevSpace.H1,
+        False,
+        degree,
+        degree,
+        dtype=default_real_type,
+    )
+    
 def create_tnt_triangle(degree):
     assert degree > 0
     # Polyset

From 0c8b6d0bd17cbfdf2d731c1069d4a73cd4b5eea3 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 18 Dec 2024 16:15:53 +0000
Subject: [PATCH 10/19] Update demo_tnt-elements.py

comments
---
 python/demo/demo_tnt-elements.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 1af034c0880..b683de1b9db 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -306,21 +306,22 @@ def poisson_error(V: fem.FunctionSpace):
 # elements of degree 1 to 9. We then do the same with Q elements of
 # degree 1 to 9.
 
-# 0 Form S elements have the following polynomial power pattern in 2D:
+# 0-Form S elements have the following polynomial power pattern for a parallelogram or quadrilateral:
 # 00000
 # 0000
 # 000
 # 00
 # 0
+# (For serendipity on general quadrilateral, precision is lost, unless one uses "direct" serendipity.)
 
-# 0 Form TNT element have the following polynomial power pattern in 2D:
+# 0-Form TNT element have the following polynomial power pattern for a parallelogram or quadrilateral:
 # 00000
 # 00000
 # 0000
 # 0000
 # 00
 
-# 0 elements have the following polynomial power pattern in 2D:
+# 0-Form Q elements have the following polynomial power pattern for a parallelogram or quadrilateral:
 # 00000
 # 00000
 # 00000

From efc0d1c160b96c091fd3f1157524bbb51fee6689 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 18 Dec 2024 16:35:38 +0000
Subject: [PATCH 11/19] use numpy.linalg.dot more

---
 python/demo/demo_tnt-elements.py | 18 ++++++------------
 1 file changed, 6 insertions(+), 12 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index b683de1b9db..30e086204f2 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -213,9 +213,7 @@ def create_tnt_quad(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            v0 = geometry[e[0]]
-            v1 = geometry[e[1]]
-            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
             x[1].append(edge_pts)
     
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
@@ -419,7 +417,7 @@ def create_tnt_hex(degree):
     for i in range(degree):
         for j in range(degree):
             for k in range(degree):
-                wcoeffs[dof_n, (i * (degree + 1) + j) * (degree + 1) + k ] = 1
+                wcoeffs[dof_n, (i * (degree + 1) + j) * (degree + 1) + k] = 1
                 dof_n += 1
 
     for i, j, k in [(degree, 0, 0), (0, degree, 0), (0, 0, degree), (1, 1, degree), (1, 0, degree), (0, 1, degree), (degree, 1, 0)]:
@@ -456,7 +454,7 @@ def create_tnt_hex(degree):
         for e in topology[1]:
             v0 = geometry[e[0]]
             v1 = geometry[e[1]]
-            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
             x[1].append(edge_pts)
     
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
@@ -579,9 +577,7 @@ def create_tnt_prism(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            v0 = geometry[e[0]]
-            v1 = geometry[e[1]]
-            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
             x[1].append(edge_pts)
     
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
@@ -714,7 +710,7 @@ def create_tnt_tetrahedron(degree):
         for e in topology[1]:
             v0 = geometry[e[0]]
             v1 = geometry[e[1]]
-            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
             x[1].append(edge_pts)
     
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
@@ -820,9 +816,7 @@ def create_tnt_triangle(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            v0 = geometry[e[0]]
-            v1 = geometry[e[1]]
-            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
+            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
             x[1].append(edge_pts)
     
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))

From 851a70ac12849da24d2bafa1fd4ca16521e2b67c Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 18 Dec 2024 17:17:53 +0000
Subject: [PATCH 12/19] Simplify edge dofs

---
 python/demo/demo_tnt-elements.py | 28 ++++++----------------------
 1 file changed, 6 insertions(+), 22 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 30e086204f2..c968e71a116 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -213,9 +213,7 @@ def create_tnt_quad(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
-            x[1].append(edge_pts)
-    
+            x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
@@ -454,9 +452,7 @@ def create_tnt_hex(degree):
         for e in topology[1]:
             v0 = geometry[e[0]]
             v1 = geometry[e[1]]
-            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
-            x[1].append(edge_pts)
-    
+            x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
@@ -577,9 +573,7 @@ def create_tnt_prism(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
-            x[1].append(edge_pts)
-    
+            x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(v1 - v0)
@@ -708,11 +702,7 @@ def create_tnt_tetrahedron(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            v0 = geometry[e[0]]
-            v1 = geometry[e[1]]
-            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
-            x[1].append(edge_pts)
-    
+            x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
@@ -816,9 +806,7 @@ def create_tnt_triangle(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            edge_pts = np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
-            x[1].append(edge_pts)
-    
+            x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(v1 - v0)
@@ -896,11 +884,7 @@ def create_tnt_interval(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            v0 = geometry[e[0]]
-            v1 = geometry[e[1]]
-            edge_pts = np.array([v0 + p * (v1 - v0) for p in pts])
-            x[1].append(edge_pts)
-    
+            x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))    
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]

From a18aab1f7e5e12f1a5811623dbde51c06955c2f0 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 18 Dec 2024 18:03:19 +0000
Subject: [PATCH 13/19] Relation to legendre serendipity

---
 python/demo/demo_tnt-elements.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index c968e71a116..978eae97640 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -670,6 +670,7 @@ def create_tnt_prism(degree):
     )
 
 # Here is the matching tetrahedron, triangle and interval elements. Their boundary dofs are moments, not point values, except at the vertices.
+# The interval elements match the Legendre variant of the serendipity elements.
 
 def create_tnt_tetrahedron(degree):
     assert degree > 0

From eee575b1a9931004ddbb84f933fb8e7fed5a8695 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 18 Dec 2024 18:58:56 +0000
Subject: [PATCH 14/19] Clean up further

---
 python/demo/demo_tnt-elements.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 978eae97640..0b0fed10b2f 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -480,7 +480,7 @@ def create_tnt_hex(degree):
                     pol_set[0]*((u-1)*u+(v-1)*v))
         face_ndofs = poly.shape[0]
         for f in topology[2]:
-            x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
+            x[2].append(np.dot(ptsr,(geometry[f[1]]-geometry[f[0]],geometry[f[2]]-geometry[f[0]]))+geometry[f[0]])
             mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
             for i in range(face_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
@@ -585,11 +585,11 @@ def create_tnt_prism(degree):
             x[2].append(np.zeros([0, 3]))
             M[2].append(np.zeros([0, 1, 0, 1]))
     else:
-        ptsr_t, wts_t = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2) # or triangle for face 0 or 5
+        ptsr_t, wts_t = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2)
         pol_set_t = basix.polynomials.tabulate_polynomial_set(
             basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 2, ptsr_t
         )
-        ptsr_q, wts_q = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2) # or triangle for face 0 or 5
+        ptsr_q, wts_q = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2)
         pol_set_q = basix.polynomials.tabulate_polynomial_set(
              basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 2, ptsr_q
             )
@@ -616,10 +616,10 @@ def create_tnt_prism(degree):
                 wts = wts_q
                 ptsr = ptsr_q
             face_ndofs = poly.shape[0]
-            x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
+            x[2].append(np.dot(ptsr,(geometry[f[1]]-geometry[f[0]],geometry[f[2]]-geometry[f[0]]))+geometry[f[0]])
             mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
             for i in range(face_ndofs):
-                mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(np.cross(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))
+                mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(np.cross(geometry[f[1]]-geometry[f[0]],geometry[f[2]]-geometry[f[0]]))
             M[2].append(mat)
 
     # Interior
@@ -628,7 +628,8 @@ def create_tnt_prism(degree):
         M[3].append(np.zeros([0, 1, 0, 1]))
     else:
         pts, wts = basix.make_quadrature(basix.CellType.prism, 2 * degree - 2)
-        #The dimension of the left over space tells us we should reduce the xy space by 1, so we are making the appopriate selection
+        #The dimension of the left over space tells us we should reduce the xy space by 1 as the xy bubble is 3rd order,
+        #so we are making the appopriate selection.
         sel=[]
         for i in range(round((degree - 2) * (degree - 3)/ 2)):
             for j in range(degree - 2):
@@ -671,6 +672,7 @@ def create_tnt_prism(degree):
 
 # Here is the matching tetrahedron, triangle and interval elements. Their boundary dofs are moments, not point values, except at the vertices.
 # The interval elements match the Legendre variant of the serendipity elements.
+# The performance of the triangle element matches the P type element, as the same space is spanned.
 
 def create_tnt_tetrahedron(degree):
     assert degree > 0
@@ -730,7 +732,7 @@ def create_tnt_tetrahedron(degree):
 
         for f in topology[2]:
             face_ndofs = poly.shape[0]
-            x[2].append(np.dot(ptsr,(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))+geometry[f][0])
+            x[2].append(np.dot(ptsr,(geometry[f[1]]-geometry[f[0]],geometry[f[2]]-geometry[f[0]]))+geometry[f[0]])
             mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
             for i in range(face_ndofs):
                 mat[i, 0, :, 0] = wts[:] * poly[i, :]
@@ -819,7 +821,7 @@ def create_tnt_triangle(degree):
             x[2].append(np.zeros([0, 2]))
             M[2].append(np.zeros([0, 1, 0, 1]))
     else:
-        ptsr, wts = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2) # or triangle for face 0 or 5
+        ptsr, wts = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2)
         pol_set = basix.polynomials.tabulate_polynomial_set(
             basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 2, ptsr
         )

From e49bc23a6e94650ccf2eb2955ebb08fec789a84c Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 18 Dec 2024 19:06:31 +0000
Subject: [PATCH 15/19] Remove superfluous

---
 python/demo/demo_tnt-elements.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 0b0fed10b2f..fc0901972a7 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -450,8 +450,6 @@ def create_tnt_hex(degree):
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
-            v0 = geometry[e[0]]
-            v1 = geometry[e[1]]
             x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))
             mat = np.zeros((edge_ndofs, 1, len(pts), 1))
             for i in range(edge_ndofs):

From ab999c4230b278ee9ea554e1ff0282a5ca77d4de Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Thu, 19 Dec 2024 17:29:43 +0000
Subject: [PATCH 16/19] 0Form and 1Form interval

---
 python/demo/demo_tnt-elements.py | 55 ++++++++++++++++++++++++++++----
 1 file changed, 48 insertions(+), 7 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index fc0901972a7..c75ecbe299d 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -669,7 +669,6 @@ def create_tnt_prism(degree):
     )
 
 # Here is the matching tetrahedron, triangle and interval elements. Their boundary dofs are moments, not point values, except at the vertices.
-# The interval elements match the Legendre variant of the serendipity elements.
 # The performance of the triangle element matches the P type element, as the same space is spanned.
 
 def create_tnt_tetrahedron(degree):
@@ -854,7 +853,10 @@ def create_tnt_triangle(degree):
         dtype=default_real_type,
     )
 
-def create_tnt_interval(degree):
+# The interval elements match the Legendre variant of the serendipity elements for the 0-form
+# The interval elements match the Legendre variant of the discontinuous P elements for the 1-Form of 1 lower degree
+
+def create_tnt_interval_0Form(degree):
     assert degree > 0
     # Polyset
     ndofs = degree + 1
@@ -883,13 +885,9 @@ def create_tnt_interval(degree):
         poly = basix.tabulate_polynomials(
             basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
         )
-        edge_ndofs = poly.shape[0]
         for e in topology[1]:
             x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))    
-            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
-            for i in range(edge_ndofs):
-                mat[i, 0, :, 0] = wts[:] * poly[i, :]
-            M[1].append(mat)
+            M[1].append((wts*[[poly]]).transpose([2,1,3,0]))
     
     return basix.ufl.custom_element(
         basix.CellType.interval,
@@ -905,3 +903,46 @@ def create_tnt_interval(degree):
         degree,
         dtype=default_real_type,
     )
+
+def create_tnt_interval_1Form(degree):
+    assert degree > 0
+    # Polyset
+    ndofs = degree
+    npoly = degree
+    
+    wcoeffs = np.eye(ndofs)
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.interval)
+    topology = basix.topology(basix.CellType.interval)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
+
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.zeros([0, 1]))
+        M[0].append(np.zeros([0, 1, 0, 1]))
+
+    # Edges
+    pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 1)
+    poly = basix.tabulate_polynomials(
+        basix.PolynomialType.legendre, basix.CellType.interval, degree - 1, pts
+    )
+    for e in topology[1]:
+        x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))   
+        M[1].append((wts*[[poly]]).transpose([2,1,3,0]))
+
+    return basix.ufl.custom_element(
+        basix.CellType.interval,
+        [],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.identity,
+        basix.SobolevSpace.H1,
+        False,
+        degree-1,
+        degree-1,
+        dtype=default_real_type,
+    )

From dd688ba80d3a55a96e5e960c52fb156cd300c5fd Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Tue, 24 Dec 2024 14:42:39 +0100
Subject: [PATCH 17/19] TNT Triangle style N1E 1-Form

---
 python/demo/demo_tnt-elements.py | 120 +++++++++++++++++++++++++------
 1 file changed, 97 insertions(+), 23 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index c75ecbe299d..a827f452339 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -669,7 +669,8 @@ def create_tnt_prism(degree):
     )
 
 # Here is the matching tetrahedron, triangle and interval elements. Their boundary dofs are moments, not point values, except at the vertices.
-# The performance of the triangle element matches the P type element, as the same space is spanned.
+# The performance of the triangle element 0-form matches the P type element, as the same space is spanned. The interior spaces are differently parametrized.
+# The performance of the triangle element 1-Form matches the Legendre variant of Nedelec1, as the same space is spanned. The interior spaces are differently parametrized.
 
 def create_tnt_tetrahedron(degree):
     assert degree > 0
@@ -775,7 +776,7 @@ def create_tnt_tetrahedron(degree):
         dtype=default_real_type,
     )
     
-def create_tnt_triangle(degree):
+def create_tnt_triangle_0Form(degree):
     assert degree > 0
     # Polyset
     ndofs = round((degree+1)*(degree+2)/2)
@@ -801,16 +802,11 @@ def create_tnt_triangle(degree):
             M[1].append(np.zeros([0, 1, 0, 1]))
     else:
         pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
-        poly = basix.tabulate_polynomials(
-            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
-        )
-        edge_ndofs = poly.shape[0]
+        poly = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts)
         for e in topology[1]:
-            x[1].append(np.array(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]])))
-            mat = np.zeros((edge_ndofs, 1, len(pts), 1))
-            for i in range(edge_ndofs):
-                mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(v1 - v0)
-            M[1].append(mat)
+            x[1].append(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
+            M[1].append((wts[None,None,:]*[[poly]]).transpose([2,1,3,0]))
+            print(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
 
     # Faces
     if degree < 3:
@@ -818,25 +814,20 @@ def create_tnt_triangle(degree):
             x[2].append(np.zeros([0, 2]))
             M[2].append(np.zeros([0, 1, 0, 1]))
     else:
-        ptsr, wts = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2)
-        pol_set = basix.polynomials.tabulate_polynomial_set(
-            basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 2, ptsr
-        )
+        pts, wts = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2)
+        pol_set = basix.polynomials.tabulate_polynomial_set(basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 2, pts)
         # this assumes the conventional [0 to 1][0 to 1] domain of the reference element, 
         # and takes the Laplacian of (u+v-1)*u*v*pol_set[0], , 
         # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py
-        u = ptsr[:, 0]
-        v = ptsr[:, 1]
+        u = pts[:, 0]
+        v = pts[:, 1]
         poly = (pol_set[5]+pol_set[3])*(u+v-1)*u*v+ \
                 2*(pol_set[2]*u*(2*v+u-1)+pol_set[1]*v*(2*u+v-1)+ \
                    pol_set[0]*(u+v))
 
-        face_ndofs = poly.shape[0]
-        x[2].append(ptsr)
-        mat = np.zeros((face_ndofs, 1, len(ptsr), 1))
-        for i in range(face_ndofs):
-            mat[i, 0, :, 0] = wts[:] * poly[i, :] #* np.linalg.norm(np.cross(geometry[f][1]-geometry[f][0],geometry[f][2]-geometry[f][0]))
-        M[2].append(mat)
+        x[2].append(pts)
+        M[2].append((wts[None,None,:]*[[poly]]).transpose([2,1,3,0]))
+
     
     return basix.ufl.custom_element(
         basix.CellType.triangle,
@@ -852,7 +843,90 @@ def create_tnt_triangle(degree):
         degree,
         dtype=default_real_type,
     )
+    
+def cross2d(x):
+    return [x[1],-x[0]] 
+
+def create_tnt_triangle_1Form (degree):
+    assert degree > 0
+    # Polyset
+    ndofs = (degree+2)*degree
+    npoly = (degree+1)*(degree+2)
+    
+    wcoeffs = np.zeros((ndofs,npoly))
+    wcoeffs[:round((degree+1)*degree/2),:round((degree+1)*degree/2)] = np.eye(round((degree+1)*degree/2))
+    wcoeffs[round((degree+1)*degree/2):(degree+1)*degree,round((degree+1)*(degree+2)/2):round((degree+1)*(degree+2)/2)+round((degree+1)*degree/2)] = np.eye(round((degree+1)*degree/2))
+
+    pts, wts = basix.make_quadrature(basix.CellType.triangle, 2*degree)
+    poly  = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.triangle, degree, pts)
+    poly0 = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.triangle, degree-1, pts)
+    u = pts[:, 0]
+    v = pts[:, 1]
+    for j in range(degree):
+        for i in range(round((degree+1)*(degree+2)/2)):
+            wcoeffs[(degree+1)*degree + j, i]                                  = sum( v * poly0[round((degree-1)*(degree)/2)+j, :] * poly[i, :] * wts)
+            wcoeffs[(degree+1)*degree + j, round((degree+1)*(degree+2)/2) + i] = sum(-u * poly0[round((degree-1)*(degree)/2)+j, :] * poly[i, :] * wts)
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.triangle)
+    topology = basix.topology(basix.CellType.triangle)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
 
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.zeros([0, 2]))
+        M[0].append(np.zeros([0, 2, 0, 1]))
+
+    # Edges
+    pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2 )
+    poly = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.interval, degree - 1, pts)
+    for e in topology[1]:
+        x[1].append(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
+        mat = np.ones((degree, 2, len(wts), 1))
+        mat0 = np.multiply.outer(geometry[e[1]]-geometry[e[0]],wts*[poly]).transpose([2,0,3,1])
+        mat1=np.zeros(mat0.shape)
+        mat1[:,:,:,:]=mat0
+        M[1].append(mat1)
+        
+    # Faces
+    if degree < 2:
+        for _ in topology[2]:
+            x[2].append(np.zeros([0, 2]))
+            M[2].append(np.zeros([0, 2, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.triangle, 2 * degree - 2)
+        x[2].append(pts)
+        pol_set = basix.polynomials.tabulate_polynomial_set(basix.CellType.triangle, basix.PolynomialType.legendre, degree - 1, 1, pts)
+        mat0=(wts[None,None,None,:]*[cross2d(pol_set[1:,1:])]).transpose([2,1,3,0])
+        mat1=np.zeros(mat0.shape)
+        mat1[:,:,:,:]=mat0
+        if degree<3:
+            M[2].append(mat1)
+        else:
+            pol_set = basix.polynomials.tabulate_polynomial_set(basix.CellType.triangle, basix.PolynomialType.legendre, degree - 3, 1, pts)
+            u = pts[:, 0]
+            v = pts[:, 1]
+            poly = cross2d(cross2d([v*(pol_set[1]*(u+v-1)*u+pol_set[0]*(2*u+v-1)),u*(pol_set[2]*(u+v-1)*v+pol_set[0]*(2*v+u-1))]))
+            mat2=(wts[None,None,None,:]*[poly]).transpose([2,1,3,0])
+            M[2].append(np.concatenate((mat1,mat2)))
+
+
+    return basix.ufl.custom_element(
+        basix.CellType.triangle,
+        [2],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.covariantPiola, 
+        basix.SobolevSpace.HCurl,
+        False,
+        degree-1,
+        degree,
+        basix.PolysetType.standard,
+    )
+    
 # The interval elements match the Legendre variant of the serendipity elements for the 0-form
 # The interval elements match the Legendre variant of the discontinuous P elements for the 1-Form of 1 lower degree
 

From 87699f92a0d1cdf4f7cb54919ae14e4a79e79867 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 25 Dec 2024 11:12:05 +0100
Subject: [PATCH 18/19] Quadrilateral 1Form

---
 python/demo/demo_tnt-elements.py | 99 ++++++++++++++++++++++++++++++--
 1 file changed, 95 insertions(+), 4 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index a827f452339..211bf84360c 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -169,7 +169,7 @@
 # arbitrary degree TNT elements.
 
 
-def create_tnt_quad(degree):
+def create_tnt_quadrilateral_0Form(degree):
     assert degree > 0
     # Polyset
     ndofs = 4*degree + max(degree-2,0)**2
@@ -337,7 +337,7 @@ def poisson_error(V: fem.FunctionSpace):
 #tnt_errors.append(poisson_error(V))
 #print(f"TNT degree 2 error: {tnt_errors[-1]}")
 for degree in range(1, 9):
-    V = fem.functionspace(msh, create_tnt_quad(degree))
+    V = fem.functionspace(msh, create_tnt_quadrilateral_0Form(degree))
     tnt_degrees.append(degree)
     tnt_ndofs.append(V.dofmap.index_map.size_global)
     tnt_errors.append(poisson_error(V))
@@ -401,10 +401,101 @@ def poisson_error(V: fem.FunctionSpace):
 
 # ![](demo_tnt-elements_ndofs_vs_error.png)
 
-# extension to hexahedrons is trivial:
+# We can also generate the 1Form version:
+
+def cross2d(x):
+    return [x[1],-x[0]] 
+
+def create_tnt_quadrilateral_1Form(degree):
+    assert degree > 0
+    # Polyset
+    ndofs = 4*degree+degree**2-1+max(0,(degree-2))**2
+    npoly = 2*(degree+1)**2
+
+    wcoeffs = np.zeros((ndofs,npoly))
+    for i in range(degree):
+        for j in range(degree):
+            wcoeffs[i*degree+j,i*(degree+1)+j]=1
+            wcoeffs[degree**2+i*degree+j,(degree+1)**2+i*(degree+1)+j]=1
+    wcoeffs[2*degree**2,degree**2+degree-1]=1
+    wcoeffs[2*degree**2,2*(degree+1)**2-2]=-1
+    
+    pts, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2*degree)
+    poly = basix.polynomials.tabulate_polynomial_set(basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree, 1, pts)
+#    u = pts[:, 0]
+#    v = pts[:, 1]
+    for i in range((degree+1)**2):
+# alternative calculation of this row for checking
+#        wcoeffs[2*degree**2,i]              =sum( v*poly[0,(degree+1)**2-degree-3, :] * poly[0,i, :] * wts)
+#        wcoeffs[2*degree**2,(degree+1)**2+i]=sum(-u*poly[0,(degree+1)**2-degree-3, :] * poly[0,i, :] * wts)
+        wcoeffs[2*degree**2+1,i]               = sum( poly[1,2*degree+1, :] * poly[0,i, :] * wts)
+        wcoeffs[2*degree**2+1,(degree+1)**2+i] = sum( poly[2,2*degree+1, :] * poly[0,i, :] * wts)
+        if degree>1:
+            wcoeffs[2*degree**2+2,i]               = sum( poly[1,(degree+1)**2-degree, :] * poly[0,i, :] * wts)
+            wcoeffs[2*degree**2+2,(degree+1)**2+i] = sum( poly[2,(degree+1)**2-degree, :] * poly[0,i, :] * wts)
+
+    # Interpolation
+    geometry = basix.geometry(basix.CellType.quadrilateral)
+    topology = basix.topology(basix.CellType.quadrilateral)
+    x = [[], [], [], []]
+    M = [[], [], [], []]
 
+    # Vertices
+    for v in topology[0]:
+        x[0].append(np.zeros([0, 2]))
+        M[0].append(np.zeros([0, 2, 0, 1]))
+
+    # Edges
+    pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2 )
+    poly = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.interval, degree - 1, pts)
+    for e in topology[1]:
+        x[1].append(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
+        mat = np.ones((degree, 2, len(wts), 1))
+        mat0 = np.multiply.outer(geometry[e[1]]-geometry[e[0]],wts*[poly]).transpose([2,0,3,1])
+        mat1=np.zeros(mat0.shape)
+        mat1[:,:,:,:]=mat0
+        M[1].append(mat1)
+
+    # Faces
+    if degree < 2:
+        for _ in topology[2]:
+            x[2].append(np.zeros([0, 2]))
+            M[2].append(np.zeros([0, 2, 0, 1]))
+    else:
+        pts, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 2)
+        x[2].append(pts)
+        pol_set = basix.polynomials.tabulate_polynomial_set(basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 1, 1, pts)
+        mat0=(wts[None,None,None,:]*[cross2d(pol_set[1:,1:])]).transpose([2,1,3,0])
+        mat1=np.zeros(mat0.shape)
+        mat1[:,:,:,:]=mat0
+        if degree<3:
+            M[2].append(mat1)
+        else:
+            pol_set = basix.polynomials.tabulate_polynomial_set(basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree - 3, 1, pts)
+            u = pts[:, 0]
+            v = pts[:, 1]
+            poly = cross2d(cross2d([v*(v-1)*(pol_set[1]*(u-1)*u+pol_set[0]*(2*u-1)),u*(u-1)*(pol_set[2]*(v-1)*v+pol_set[0]*(2*v-1))]))
+            mat2=(wts[None,None,None,:]*[poly]).transpose([2,1,3,0])
+            M[2].append(np.concatenate((mat1,mat2)))
+
+    return basix.ufl.custom_element(
+        basix.CellType.quadrilateral,
+        [2],
+        wcoeffs,
+        x,
+        M,
+        0,
+        basix.MapType.covariantPiola, 
+        basix.SobolevSpace.HCurl,
+        False,
+        degree-1,
+        degree,
+        basix.PolysetType.standard,
+    )
+
+# extension to hexahedrons is trivial:
 
-def create_tnt_hex(degree):
+def create_tnt_hexahedron(degree):
     assert degree > 0
     # Polyset
     ndofs = 12 * degree - 4 + (degree + 4) * max(degree-2, 0) ** 2

From 8c565124d56e27246102647efe1d420cfbad9d53 Mon Sep 17 00:00:00 2001
From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com>
Date: Wed, 25 Dec 2024 16:10:36 +0100
Subject: [PATCH 19/19] Update demo_tnt-elements.py

---
 python/demo/demo_tnt-elements.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py
index 211bf84360c..a3962da322d 100644
--- a/python/demo/demo_tnt-elements.py
+++ b/python/demo/demo_tnt-elements.py
@@ -450,7 +450,6 @@ def create_tnt_quadrilateral_1Form(degree):
     poly = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.interval, degree - 1, pts)
     for e in topology[1]:
         x[1].append(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
-        mat = np.ones((degree, 2, len(wts), 1))
         mat0 = np.multiply.outer(geometry[e[1]]-geometry[e[0]],wts*[poly]).transpose([2,0,3,1])
         mat1=np.zeros(mat0.shape)
         mat1[:,:,:,:]=mat0
@@ -974,7 +973,6 @@ def create_tnt_triangle_1Form (degree):
     poly = basix.tabulate_polynomials(basix.PolynomialType.legendre, basix.CellType.interval, degree - 1, pts)
     for e in topology[1]:
         x[1].append(geometry[e[0]] + np.dot(pts,[geometry[e[1]]-geometry[e[0]]]))
-        mat = np.ones((degree, 2, len(wts), 1))
         mat0 = np.multiply.outer(geometry[e[1]]-geometry[e[0]],wts*[poly]).transpose([2,0,3,1])
         mat1=np.zeros(mat0.shape)
         mat1[:,:,:,:]=mat0