Update demo_tnt-elements.py

Shift degree by 1 for TNT

python/demo/demo_tnt-elements.py

@@ -60,7 +60,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 1 element, the polynomial set
+# polynomials on the cell. For a degree 2 element, 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 +146,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,25 +168,25 @@
 
 
 def create_tnt_quad(degree):
-    assert degree > -1
+    assert degree > 0
     # Polyset
-    ndofs = 4*(degree+1) + max(degree-1,0)**2
-    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), (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 > 0:
-        wcoeffs[dof_n, 2 * degree + 3] = 1
+    if degree > 1:
+        wcoeffs[dof_n, 2 * degree + 1] = 1
 
     # Interpolation
     geometry = basix.geometry(basix.CellType.quadrilateral)
@@ -200,14 +200,14 @@ def create_tnt_quad(degree):
         M[0].append(np.array([[[[1.0]]]]))
 
     # Edges
-    if degree < 1:
+    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)
+        pts, wts = basix.make_quadrature(basix.CellType.interval, 2 * degree - 2)
         poly = basix.tabulate_polynomials(
-            basix.PolynomialType.legendre, basix.CellType.interval, degree - 1, pts
+            basix.PolynomialType.legendre, basix.CellType.interval, degree - 2, pts
         )
         edge_ndofs = poly.shape[0]
         for e in topology[1]:
@@ -222,15 +222,15 @@ def create_tnt_quad(degree):
             M[1].append(mat)
 
     # Interior
-    if degree < 2:
+    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)
+        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-2, 2, pts
+            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], 
@@ -255,8 +255,8 @@ def create_tnt_quad(degree):
         basix.MapType.identity,
         basix.SobolevSpace.H1,
         False,
+        degree - 1,
         degree,
-        degree + 1,
         dtype=default_real_type,
     )
 
@@ -312,12 +312,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 1 error: {tnt_errors[-1]}")
-for degree in range(0, 8):
+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_ndofs.append(V.dofmap.index_map.size_global)