diff --git a/scripts/generate_gn.py b/scripts/generate_gn.py
index c026b7d..f78d13f 100644
--- a/scripts/generate_gn.py
+++ b/scripts/generate_gn.py
@@ -4,12 +4,10 @@
 import sys
 import typing
 
-import numpy as np
 import sympy
 
 import symfem
-from symfem.elements.guzman_neilan import make_piecewise_lagrange
-from symfem.functions import MatrixFunction, VectorFunction
+from symfem.functions import VectorFunction, AnyFunction
 from symfem.piecewise_functions import PiecewiseFunction
 from symfem.symbols import t, x
 
@@ -17,6 +15,7 @@
 
 
 def poly(reference, k):
+    """Generate the P^perp polynomial set."""
     if k < 2:
         if reference.name == "triangle":
             return [
@@ -155,6 +154,7 @@ def poly(reference, k):
 
 
 def find_solution(mat, aim):
+    """Solve mat @ x = aim."""
     s_mat = sympy.Matrix(mat)
     solution = (s_mat.T @ s_mat).inv() @ s_mat.T @ sympy.Matrix(aim)
     assert s_mat @ solution == sympy.Matrix(aim)
@@ -259,10 +259,12 @@ def find_solution(mat, aim):
                     i += 1
 
         coeffs = find_solution(mat, aim)
-        bubble = f
+        bubble: AnyFunction = f
         for i, j in zip(coeffs, sub_basis):
             bubble -= i * j
 
+        assert isinstance(bubble, PiecewiseFunction)
+
         output += "    {\n"
         for cell, f in bubble.pieces.items():
             output += "        "
@@ -273,9 +275,7 @@ def find_solution(mat, aim):
                         "("
                         + ", ".join(
                             [
-                                f"{c}"
-                                if isinstance(c, sympy.Integer)
-                                else f"sympy.S('{c}')"
+                                f"{c}" if isinstance(c, sympy.Integer) else f"sympy.S('{c}')"
                                 for c in p
                             ]
                         )
@@ -288,7 +288,7 @@ def find_solution(mat, aim):
             output += ": (\n"
             output += ",\n".join(
                 [
-                    "            sympy.S('" + f"{c.as_sympy().expand()}".replace(" ", "") + "')"
+                    "            sympy.S('" + f"{c.as_sympy().expand()}".replace(" ", "") + "')"  # type:ignore
                     for c in f
                 ]
             )
diff --git a/symfem/basis_functions.py b/symfem/basis_functions.py
index 42859fe..447aa53 100644
--- a/symfem/basis_functions.py
+++ b/symfem/basis_functions.py
@@ -8,7 +8,6 @@
 import sympy
 
 import symfem
-
 from symfem.functions import (
     AnyFunction,
     FunctionInput,
diff --git a/symfem/elements/_guzman_neilan_tetrahedron.py b/symfem/elements/_guzman_neilan_tetrahedron.py
index ab04f28..937499e 100644
--- a/symfem/elements/_guzman_neilan_tetrahedron.py
+++ b/symfem/elements/_guzman_neilan_tetrahedron.py
@@ -4,75 +4,127 @@
 
 bubbles = [
     {
-        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('40*x**2*z+120*x*y*z+88*x*z**2-64*x*z-220*y*z**2-860*z**3/3+128*z**2+6*z'),
-            sympy.S('120*x*y*z-220*x*z**2+40*y**2*z+88*y*z**2-64*y*z-860*z**3/3+128*z**2+6*z'),
-            sympy.S('-100*x*z**2-100*y*z**2-176*z**3/3+64*z**2+6*z')),
-        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('40*x**2*y+88*x*y**2+120*x*y*z-64*x*y-860*y**3/3-220*y**2*z+128*y**2+6*y'),
-            sympy.S('-100*x*y**2-176*y**3/3-100*y**2*z+64*y**2+6*y'),
-            sympy.S('-220*x*y**2+120*x*y*z-860*y**3/3+88*y**2*z+128*y**2+40*y*z**2-64*y*z+6*y')),
-        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('-176*x**3/3-100*x**2*y-100*x**2*z+64*x**2+6*x'),
-            sympy.S('-860*x**3/3+88*x**2*y-220*x**2*z+128*x**2+40*x*y**2+120*x*y*z-64*x*y+6*x'),
-            sympy.S('-860*x**3/3-220*x**2*y+88*x**2*z+128*x**2+120*x*y*z+40*x*z**2-64*x*z+6*x')),
-        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('224*x**3/3+36*x**2*y+36*x**2*z-104*x**2-112*x*y**2-144*x*y*z+104*x*y-112*x*z**2+104*x*z+2*x-220*y**3/3-180*y**2*z+168*y**2-180*y*z**2+296*y*z-122*y-220*z**3/3+168*z**2-122*z+82/3'),
-            sympy.S('-220*x**3/3-112*x**2*y-180*x**2*z+168*x**2+36*x*y**2-144*x*y*z+104*x*y-180*x*z**2+296*x*z-122*x+224*y**3/3+36*y**2*z-104*y**2-112*y*z**2+104*y*z+2*y-220*z**3/3+168*z**2-122*z+82/3'),
-            sympy.S('-220*x**3/3-180*x**2*y-112*x**2*z+168*x**2-180*x*y**2-144*x*y*z+296*x*y+36*x*z**2+104*x*z-122*x-220*y**3/3-112*y**2*z+168*y**2+36*y*z**2+104*y*z-122*y+224*z**3/3-104*z**2+2*z+82/3')),
+        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("40*x**2*z+120*x*y*z+88*x*z**2-64*x*z-220*y*z**2-860*z**3/3+128*z**2+6*z"),
+            sympy.S("120*x*y*z-220*x*z**2+40*y**2*z+88*y*z**2-64*y*z-860*z**3/3+128*z**2+6*z"),
+            sympy.S("-100*x*z**2-100*y*z**2-176*z**3/3+64*z**2+6*z"),
+        ),
+        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("40*x**2*y+88*x*y**2+120*x*y*z-64*x*y-860*y**3/3-220*y**2*z+128*y**2+6*y"),
+            sympy.S("-100*x*y**2-176*y**3/3-100*y**2*z+64*y**2+6*y"),
+            sympy.S("-220*x*y**2+120*x*y*z-860*y**3/3+88*y**2*z+128*y**2+40*y*z**2-64*y*z+6*y"),
+        ),
+        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("-176*x**3/3-100*x**2*y-100*x**2*z+64*x**2+6*x"),
+            sympy.S("-860*x**3/3+88*x**2*y-220*x**2*z+128*x**2+40*x*y**2+120*x*y*z-64*x*y+6*x"),
+            sympy.S("-860*x**3/3-220*x**2*y+88*x**2*z+128*x**2+120*x*y*z+40*x*z**2-64*x*z+6*x"),
+        ),
+        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S(
+                "224*x**3/3+36*x**2*y+36*x**2*z-104*x**2-112*x*y**2-144*x*y*z+104*x*y-112*x*z**2+104*x*z+2*x-220*y**3/3-180*y**2*z+168*y**2-180*y*z**2+296*y*z-122*y-220*z**3/3+168*z**2-122*z+82/3"
+            ),
+            sympy.S(
+                "-220*x**3/3-112*x**2*y-180*x**2*z+168*x**2+36*x*y**2-144*x*y*z+104*x*y-180*x*z**2+296*x*z-122*x+224*y**3/3+36*y**2*z-104*y**2-112*y*z**2+104*y*z+2*y-220*z**3/3+168*z**2-122*z+82/3"
+            ),
+            sympy.S(
+                "-220*x**3/3-180*x**2*y-112*x**2*z+168*x**2-180*x*y**2-144*x*y*z+296*x*y+36*x*z**2+104*x*z-122*x-220*y**3/3-112*y**2*z+168*y**2+36*y*z**2+104*y*z-122*y+224*z**3/3-104*z**2+2*z+82/3"
+            ),
+        ),
     },
     {
-        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('-180*x*y*z+432*x*z**2+24*x*z-120*y**2*z+210*y*z**2+120*y*z-110*z**3-108*z**2-6*z'),
-            sympy.S('-60*y**2*z-438*y*z**2+84*y*z+650*z**3-138*z**2-6*z'),
-            sympy.S('150*y*z**2+2*z**3-54*z**2-6*z')),
-        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('432*x*y**2-180*x*y*z+24*x*y-110*y**3+210*y**2*z-108*y**2-120*y*z**2+120*y*z-6*y'),
-            sympy.S('2*y**3+150*y**2*z-54*y**2-6*y'),
-            sympy.S('650*y**3-438*y**2*z-138*y**2-60*y*z**2+84*y*z-6*y')),
-        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('172*x**3+150*x**2*y+150*x**2*z-84*x**2-6*x-120*y**2*z-120*y*z**2+120*y*z'),
-            sympy.S('320*x**3-258*x**2*y+330*x**2*z-138*x**2-60*x*y**2-180*x*y*z+84*x*y-6*x'),
-            sympy.S('320*x**3+330*x**2*y-258*x**2*z-138*x**2-180*x*y*z-60*x*z**2+84*x*z-6*x')),
-        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('812*x**3+2214*x**2*y+2214*x**2*z-2076*x**2+1932*x*y**2+3624*x*y*z-3624*x*y+1932*x*z**2-3624*x*z+1698*x+530*y**3+1350*y**2*z-1488*y**2+1350*y*z**2-2736*y*z+1392*y+530*z**3-1488*z**2+1392*z-434'),
-            sympy.S('-320*x**3-1218*x**2*y-630*x**2*z+822*x**2-1416*x*y**2-1596*x*y*z+2076*x*y-300*x*z**2+984*x*z-678*x-518*y**3-906*y**2*z+1194*y**2-378*y*z**2+1236*y*z-852*y+10*z**3+162*z**2-348*z+176'),
-            sympy.S('-320*x**3-630*x**2*y-1218*x**2*z+822*x**2-300*x*y**2-1596*x*y*z+984*x*y-1416*x*z**2+2076*x*z-678*x+10*y**3-378*y**2*z+162*y**2-906*y*z**2+1236*y*z-348*y-518*z**3+1194*z**2-852*z+176')),
+        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S(
+                "-180*x*y*z+432*x*z**2+24*x*z-120*y**2*z+210*y*z**2+120*y*z-110*z**3-108*z**2-6*z"
+            ),
+            sympy.S("-60*y**2*z-438*y*z**2+84*y*z+650*z**3-138*z**2-6*z"),
+            sympy.S("150*y*z**2+2*z**3-54*z**2-6*z"),
+        ),
+        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S(
+                "432*x*y**2-180*x*y*z+24*x*y-110*y**3+210*y**2*z-108*y**2-120*y*z**2+120*y*z-6*y"
+            ),
+            sympy.S("2*y**3+150*y**2*z-54*y**2-6*y"),
+            sympy.S("650*y**3-438*y**2*z-138*y**2-60*y*z**2+84*y*z-6*y"),
+        ),
+        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("172*x**3+150*x**2*y+150*x**2*z-84*x**2-6*x-120*y**2*z-120*y*z**2+120*y*z"),
+            sympy.S("320*x**3-258*x**2*y+330*x**2*z-138*x**2-60*x*y**2-180*x*y*z+84*x*y-6*x"),
+            sympy.S("320*x**3+330*x**2*y-258*x**2*z-138*x**2-180*x*y*z-60*x*z**2+84*x*z-6*x"),
+        ),
+        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S(
+                "812*x**3+2214*x**2*y+2214*x**2*z-2076*x**2+1932*x*y**2+3624*x*y*z-3624*x*y+1932*x*z**2-3624*x*z+1698*x+530*y**3+1350*y**2*z-1488*y**2+1350*y*z**2-2736*y*z+1392*y+530*z**3-1488*z**2+1392*z-434"
+            ),
+            sympy.S(
+                "-320*x**3-1218*x**2*y-630*x**2*z+822*x**2-1416*x*y**2-1596*x*y*z+2076*x*y-300*x*z**2+984*x*z-678*x-518*y**3-906*y**2*z+1194*y**2-378*y*z**2+1236*y*z-852*y+10*z**3+162*z**2-348*z+176"
+            ),
+            sympy.S(
+                "-320*x**3-630*x**2*y-1218*x**2*z+822*x**2-300*x*y**2-1596*x*y*z+984*x*y-1416*x*z**2+2076*x*z-678*x+10*y**3-378*y**2*z+162*y**2-906*y*z**2+1236*y*z-348*y-518*z**3+1194*z**2-852*z+176"
+            ),
+        ),
     },
     {
-        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('60*x**2*z+438*x*z**2-84*x*z-650*z**3+138*z**2+6*z'),
-            sympy.S('120*x**2*z+180*x*y*z-210*x*z**2-120*x*z-432*y*z**2-24*y*z+110*z**3+108*z**2+6*z'),
-            sympy.S('-150*x*z**2-2*z**3+54*z**2+6*z')),
-        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('60*x**2*y+258*x*y**2+180*x*y*z-84*x*y-320*y**3-330*y**2*z+138*y**2+6*y'),
-            sympy.S('120*x**2*z-150*x*y**2+120*x*z**2-120*x*z-172*y**3-150*y**2*z+84*y**2+6*y'),
-            sympy.S('-330*x*y**2+180*x*y*z-320*y**3+258*y**2*z+138*y**2+60*y*z**2-84*y*z+6*y')),
-        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('-2*x**3-150*x**2*z+54*x**2+6*x'),
-            sympy.S('110*x**3-432*x**2*y-210*x**2*z+108*x**2+180*x*y*z-24*x*y+120*x*z**2-120*x*z+6*x'),
-            sympy.S('-650*x**3+438*x**2*z+138*x**2+60*x*z**2-84*x*z+6*x')),
-        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('518*x**3+1416*x**2*y+906*x**2*z-1194*x**2+1218*x*y**2+1596*x*y*z-2076*x*y+378*x*z**2-1236*x*z+852*x+320*y**3+630*y**2*z-822*y**2+300*y*z**2-984*y*z+678*y-10*z**3-162*z**2+348*z-176'),
-            sympy.S('-530*x**3-1932*x**2*y-1350*x**2*z+1488*x**2-2214*x*y**2-3624*x*y*z+3624*x*y-1350*x*z**2+2736*x*z-1392*x-812*y**3-2214*y**2*z+2076*y**2-1932*y*z**2+3624*y*z-1698*y-530*z**3+1488*z**2-1392*z+434'),
-            sympy.S('-10*x**3+300*x**2*y+378*x**2*z-162*x**2+630*x*y**2+1596*x*y*z-984*x*y+906*x*z**2-1236*x*z+348*x+320*y**3+1218*y**2*z-822*y**2+1416*y*z**2-2076*y*z+678*y+518*z**3-1194*z**2+852*z-176')),
+        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("60*x**2*z+438*x*z**2-84*x*z-650*z**3+138*z**2+6*z"),
+            sympy.S(
+                "120*x**2*z+180*x*y*z-210*x*z**2-120*x*z-432*y*z**2-24*y*z+110*z**3+108*z**2+6*z"
+            ),
+            sympy.S("-150*x*z**2-2*z**3+54*z**2+6*z"),
+        ),
+        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("60*x**2*y+258*x*y**2+180*x*y*z-84*x*y-320*y**3-330*y**2*z+138*y**2+6*y"),
+            sympy.S("120*x**2*z-150*x*y**2+120*x*z**2-120*x*z-172*y**3-150*y**2*z+84*y**2+6*y"),
+            sympy.S("-330*x*y**2+180*x*y*z-320*y**3+258*y**2*z+138*y**2+60*y*z**2-84*y*z+6*y"),
+        ),
+        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("-2*x**3-150*x**2*z+54*x**2+6*x"),
+            sympy.S(
+                "110*x**3-432*x**2*y-210*x**2*z+108*x**2+180*x*y*z-24*x*y+120*x*z**2-120*x*z+6*x"
+            ),
+            sympy.S("-650*x**3+438*x**2*z+138*x**2+60*x*z**2-84*x*z+6*x"),
+        ),
+        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S(
+                "518*x**3+1416*x**2*y+906*x**2*z-1194*x**2+1218*x*y**2+1596*x*y*z-2076*x*y+378*x*z**2-1236*x*z+852*x+320*y**3+630*y**2*z-822*y**2+300*y*z**2-984*y*z+678*y-10*z**3-162*z**2+348*z-176"
+            ),
+            sympy.S(
+                "-530*x**3-1932*x**2*y-1350*x**2*z+1488*x**2-2214*x*y**2-3624*x*y*z+3624*x*y-1350*x*z**2+2736*x*z-1392*x-812*y**3-2214*y**2*z+2076*y**2-1932*y*z**2+3624*y*z-1698*y-530*z**3+1488*z**2-1392*z+434"
+            ),
+            sympy.S(
+                "-10*x**3+300*x**2*y+378*x**2*z-162*x**2+630*x*y**2+1596*x*y*z-984*x*y+906*x*z**2-1236*x*z+348*x+320*y**3+1218*y**2*z-822*y**2+1416*y*z**2-2076*y*z+678*y+518*z**3-1194*z**2+852*z-176"
+            ),
+        ),
     },
     {
-        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('-60*x**2*z-180*x*y*z-258*x*z**2+84*x*z+330*y*z**2+320*z**3-138*z**2-6*z'),
-            sympy.S('-180*x*y*z+330*x*z**2-60*y**2*z-258*y*z**2+84*y*z+320*z**3-138*z**2-6*z'),
-            sympy.S('-120*x**2*y-120*x*y**2+120*x*y+150*x*z**2+150*y*z**2+172*z**3-84*z**2-6*z')),
-        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('-60*x**2*y-438*x*y**2+84*x*y+650*y**3-138*y**2-6*y'),
-            sympy.S('150*x*y**2+2*y**3-54*y**2-6*y'),
-            sympy.S('-120*x**2*y+210*x*y**2-180*x*y*z+120*x*y-110*y**3+432*y**2*z-108*y**2+24*y*z-6*y')),
-        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('2*x**3+150*x**2*y-54*x**2-6*x'),
-            sympy.S('650*x**3-438*x**2*y-138*x**2-60*x*y**2+84*x*y-6*x'),
-            sympy.S('-110*x**3+210*x**2*y+432*x**2*z-108*x**2-120*x*y**2-180*x*y*z+120*x*y+24*x*z-6*x')),
-        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S('1/4'), sympy.S('1/4'), sympy.S('1/4'))): (
-            sympy.S('-518*x**3-906*x**2*y-1416*x**2*z+1194*x**2-378*x*y**2-1596*x*y*z+1236*x*y-1218*x*z**2+2076*x*z-852*x+10*y**3-300*y**2*z+162*y**2-630*y*z**2+984*y*z-348*y-320*z**3+822*z**2-678*z+176'),
-            sympy.S('10*x**3-378*x**2*y-300*x**2*z+162*x**2-906*x*y**2-1596*x*y*z+1236*x*y-630*x*z**2+984*x*z-348*x-518*y**3-1416*y**2*z+1194*y**2-1218*y*z**2+2076*y*z-852*y-320*z**3+822*z**2-678*z+176'),
-            sympy.S('530*x**3+1350*x**2*y+1932*x**2*z-1488*x**2+1350*x*y**2+3624*x*y*z-2736*x*y+2214*x*z**2-3624*x*z+1392*x+530*y**3+1932*y**2*z-1488*y**2+2214*y*z**2-3624*y*z+1392*y+812*z**3-2076*z**2+1698*z-434')),
+        ((0, 0, 0), (1, 0, 0), (0, 1, 0), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("-60*x**2*z-180*x*y*z-258*x*z**2+84*x*z+330*y*z**2+320*z**3-138*z**2-6*z"),
+            sympy.S("-180*x*y*z+330*x*z**2-60*y**2*z-258*y*z**2+84*y*z+320*z**3-138*z**2-6*z"),
+            sympy.S("-120*x**2*y-120*x*y**2+120*x*y+150*x*z**2+150*y*z**2+172*z**3-84*z**2-6*z"),
+        ),
+        ((0, 0, 0), (1, 0, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("-60*x**2*y-438*x*y**2+84*x*y+650*y**3-138*y**2-6*y"),
+            sympy.S("150*x*y**2+2*y**3-54*y**2-6*y"),
+            sympy.S(
+                "-120*x**2*y+210*x*y**2-180*x*y*z+120*x*y-110*y**3+432*y**2*z-108*y**2+24*y*z-6*y"
+            ),
+        ),
+        ((0, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S("2*x**3+150*x**2*y-54*x**2-6*x"),
+            sympy.S("650*x**3-438*x**2*y-138*x**2-60*x*y**2+84*x*y-6*x"),
+            sympy.S(
+                "-110*x**3+210*x**2*y+432*x**2*z-108*x**2-120*x*y**2-180*x*y*z+120*x*y+24*x*z-6*x"
+            ),
+        ),
+        ((1, 0, 0), (0, 1, 0), (0, 0, 1), (sympy.S("1/4"), sympy.S("1/4"), sympy.S("1/4"))): (
+            sympy.S(
+                "-518*x**3-906*x**2*y-1416*x**2*z+1194*x**2-378*x*y**2-1596*x*y*z+1236*x*y-1218*x*z**2+2076*x*z-852*x+10*y**3-300*y**2*z+162*y**2-630*y*z**2+984*y*z-348*y-320*z**3+822*z**2-678*z+176"
+            ),
+            sympy.S(
+                "10*x**3-378*x**2*y-300*x**2*z+162*x**2-906*x*y**2-1596*x*y*z+1236*x*y-630*x*z**2+984*x*z-348*x-518*y**3-1416*y**2*z+1194*y**2-1218*y*z**2+2076*y*z-852*y-320*z**3+822*z**2-678*z+176"
+            ),
+            sympy.S(
+                "530*x**3+1350*x**2*y+1932*x**2*z-1488*x**2+1350*x*y**2+3624*x*y*z-2736*x*y+2214*x*z**2-3624*x*z+1392*x+530*y**3+1932*y**2*z-1488*y**2+2214*y*z**2-3624*y*z+1392*y+812*z**3-2076*z**2+1698*z-434"
+            ),
+        ),
     },
 ]
diff --git a/symfem/elements/_guzman_neilan_triangle.py b/symfem/elements/_guzman_neilan_triangle.py
index 6246c4f..4ce6d5a 100644
--- a/symfem/elements/_guzman_neilan_triangle.py
+++ b/symfem/elements/_guzman_neilan_triangle.py
@@ -4,36 +4,42 @@
 
 bubbles = [
     {
-        ((0, 0), (1, 0), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('-3*x*y+9*y**2/2-2*y'),
-            sympy.S('3*y**2/2-2*y')),
-        ((0, 0), (0, 1), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('3*x**2/2-2*x'),
-            sympy.S('9*x**2/2-3*x*y-2*x')),
-        ((1, 0), (0, 1), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('3*x**2/2-3*x*y-x-3*y**2/2+2*y-1/2'),
-            sympy.S('-3*x**2/2-3*x*y+2*x+3*y**2/2-y-1/2')),
+        ((0, 0), (1, 0), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("-3*x*y+9*y**2/2-2*y"),
+            sympy.S("3*y**2/2-2*y"),
+        ),
+        ((0, 0), (0, 1), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("3*x**2/2-2*x"),
+            sympy.S("9*x**2/2-3*x*y-2*x"),
+        ),
+        ((1, 0), (0, 1), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("3*x**2/2-3*x*y-x-3*y**2/2+2*y-1/2"),
+            sympy.S("-3*x**2/2-3*x*y+2*x+3*y**2/2-y-1/2"),
+        ),
     },
     {
-        ((0, 0), (1, 0), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('3*y**2-4*y'),
-            sympy.S('2*y')),
-        ((0, 0), (0, 1), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('-3*x**2+2*x+6*y**2-6*y'),
-            sympy.S('-6*x**2+6*x*y+2*x')),
-        ((1, 0), (0, 1), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('9*x**2+24*x*y-14*x+15*y**2-20*y+5'),
-            sympy.S('-6*x**2-18*x*y+10*x-12*y**2+16*y-4')),
+        ((0, 0), (1, 0), (sympy.S("1/3"), sympy.S("1/3"))): (sympy.S("3*y**2-4*y"), sympy.S("2*y")),
+        ((0, 0), (0, 1), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("-3*x**2+2*x+6*y**2-6*y"),
+            sympy.S("-6*x**2+6*x*y+2*x"),
+        ),
+        ((1, 0), (0, 1), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("9*x**2+24*x*y-14*x+15*y**2-20*y+5"),
+            sympy.S("-6*x**2-18*x*y+10*x-12*y**2+16*y-4"),
+        ),
     },
     {
-        ((0, 0), (1, 0), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('-6*x*y+6*y**2-2*y'),
-            sympy.S('-6*x**2+6*x+3*y**2-2*y')),
-        ((0, 0), (0, 1), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('-2*x'),
-            sympy.S('-3*x**2+4*x')),
-        ((1, 0), (0, 1), (sympy.S('1/3'), sympy.S('1/3'))): (
-            sympy.S('12*x**2+18*x*y-16*x+6*y**2-10*y+4'),
-            sympy.S('-15*x**2-24*x*y+20*x-9*y**2+14*y-5')),
+        ((0, 0), (1, 0), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("-6*x*y+6*y**2-2*y"),
+            sympy.S("-6*x**2+6*x+3*y**2-2*y"),
+        ),
+        ((0, 0), (0, 1), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("-2*x"),
+            sympy.S("-3*x**2+4*x"),
+        ),
+        ((1, 0), (0, 1), (sympy.S("1/3"), sympy.S("1/3"))): (
+            sympy.S("12*x**2+18*x*y-16*x+6*y**2-10*y+4"),
+            sympy.S("-15*x**2-24*x*y+20*x-9*y**2+14*y-5"),
+        ),
     },
 ]
diff --git a/symfem/elements/abf.py b/symfem/elements/abf.py
index 864e7a4..12c014c 100644
--- a/symfem/elements/abf.py
+++ b/symfem/elements/abf.py
@@ -6,6 +6,8 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange
+from symfem.elements.q import Nedelec
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     IntegralMoment,
@@ -17,8 +19,6 @@
 from symfem.moments import make_integral_moment_dofs
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import Lagrange
-from symfem.elements.q import Nedelec
 
 __all__ = ["ArnoldBoffiFalk"]
 
diff --git a/symfem/elements/ac.py b/symfem/elements/ac.py
index d34c6d1..b439fe5 100644
--- a/symfem/elements/ac.py
+++ b/symfem/elements/ac.py
@@ -6,6 +6,7 @@
 
 import typing
 
+from symfem.elements.dpc import DPC
 from symfem.finite_element import CiarletElement
 from symfem.functionals import IntegralAgainst, ListOfFunctionals, NormalIntegralMoment
 from symfem.functions import FunctionInput
@@ -13,7 +14,6 @@
 from symfem.polynomials import Hdiv_serendipity, polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.dpc import DPC
 
 __all__ = ["AC"]
 
diff --git a/symfem/elements/aw.py b/symfem/elements/aw.py
index 6477459..d50b084 100644
--- a/symfem/elements/aw.py
+++ b/symfem/elements/aw.py
@@ -9,6 +9,7 @@
 
 import sympy
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     InnerProductIntegralMoment,
@@ -20,7 +21,6 @@
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["ArnoldWinther", "NonConformingArnoldWinther"]
 
diff --git a/symfem/elements/bddf.py b/symfem/elements/bddf.py
index 7d0cdd5..bc9d76a 100644
--- a/symfem/elements/bddf.py
+++ b/symfem/elements/bddf.py
@@ -6,6 +6,7 @@
 
 import typing
 
+from symfem.elements.dpc import DPC, VectorDPC
 from symfem.finite_element import CiarletElement
 from symfem.functionals import IntegralMoment, ListOfFunctionals, NormalIntegralMoment
 from symfem.functions import FunctionInput, VectorFunction
@@ -13,7 +14,6 @@
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.dpc import DPC, VectorDPC
 
 __all__ = ["bddf_polyset", "BDDF"]
 
diff --git a/symfem/elements/bdfm.py b/symfem/elements/bdfm.py
index a69589e..ce611f0 100644
--- a/symfem/elements/bdfm.py
+++ b/symfem/elements/bdfm.py
@@ -7,6 +7,9 @@
 
 import typing
 
+from symfem.elements.dpc import DPC, VectorDPC
+from symfem.elements.lagrange import Lagrange
+from symfem.elements.nedelec import NedelecFirstKind
 from symfem.finite_element import CiarletElement
 from symfem.functionals import IntegralMoment, ListOfFunctionals, NormalIntegralMoment
 from symfem.functions import FunctionInput
@@ -14,9 +17,6 @@
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.dpc import DPC, VectorDPC
-from symfem.elements.lagrange import Lagrange
-from symfem.elements.nedelec import NedelecFirstKind
 
 __all__ = ["bdfm_polyset", "BDFM"]
 
diff --git a/symfem/elements/bdm.py b/symfem/elements/bdm.py
index 05d0f0e..3b71ba2 100644
--- a/symfem/elements/bdm.py
+++ b/symfem/elements/bdm.py
@@ -6,14 +6,14 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange
+from symfem.elements.nedelec import NedelecFirstKind
 from symfem.finite_element import CiarletElement
 from symfem.functionals import IntegralMoment, ListOfFunctionals, NormalIntegralMoment
 from symfem.functions import FunctionInput
 from symfem.moments import make_integral_moment_dofs
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
-from symfem.elements.lagrange import Lagrange
-from symfem.elements.nedelec import NedelecFirstKind
 
 __all__ = ["BDM"]
 
diff --git a/symfem/elements/bernardi_raugel.py b/symfem/elements/bernardi_raugel.py
index 6f6ece6..2bcdd71 100644
--- a/symfem/elements/bernardi_raugel.py
+++ b/symfem/elements/bernardi_raugel.py
@@ -8,6 +8,7 @@
 
 import sympy
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     DivergenceIntegralMoment,
@@ -20,7 +21,6 @@
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["BernardiRaugel"]
 
diff --git a/symfem/elements/bubble.py b/symfem/elements/bubble.py
index 28cf61c..5fd2f54 100644
--- a/symfem/elements/bubble.py
+++ b/symfem/elements/bubble.py
@@ -9,11 +9,11 @@
 
 import sympy
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import DotPointEvaluation, ListOfFunctionals, PointEvaluation
 from symfem.functions import FunctionInput
 from symfem.references import Reference
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["Bubble", "BubbleEnrichedLagrange", "BubbleEnrichedVectorLagrange"]
 
diff --git a/symfem/elements/direct_serendipity.py b/symfem/elements/direct_serendipity.py
index 9694551..97d7905 100644
--- a/symfem/elements/direct_serendipity.py
+++ b/symfem/elements/direct_serendipity.py
@@ -6,10 +6,10 @@
 
 import typing
 
+from symfem.elements.dpc import DPC
 from symfem.finite_element import DirectElement
 from symfem.references import Reference
 from symfem.symbols import x
-from symfem.elements.dpc import DPC
 
 __all__ = ["DirectSerendipity"]
 
diff --git a/symfem/elements/dpc.py b/symfem/elements/dpc.py
index 4c704c9..af1e9fa 100644
--- a/symfem/elements/dpc.py
+++ b/symfem/elements/dpc.py
@@ -5,6 +5,7 @@
 
 import sympy
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     DotPointEvaluation,
@@ -15,7 +16,6 @@
 from symfem.functions import FunctionInput
 from symfem.polynomials import polynomial_set_1d, polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["DPC", "VectorDPC"]
 
diff --git a/symfem/elements/enriched_galerkin.py b/symfem/elements/enriched_galerkin.py
index 51b30fa..0c34ae9 100644
--- a/symfem/elements/enriched_galerkin.py
+++ b/symfem/elements/enriched_galerkin.py
@@ -6,10 +6,10 @@
 
 import typing
 
-from symfem.finite_element import EnrichedElement
-from symfem.references import Reference
 from symfem.elements.lagrange import Lagrange
 from symfem.elements.q import Q
+from symfem.finite_element import EnrichedElement
+from symfem.references import Reference
 
 __all__ = ["EnrichedGalerkin"]
 
diff --git a/symfem/elements/guzman_neilan.py b/symfem/elements/guzman_neilan.py
index b6e3b7c..4c6a8a0 100644
--- a/symfem/elements/guzman_neilan.py
+++ b/symfem/elements/guzman_neilan.py
@@ -8,6 +8,8 @@
 
 import sympy
 
+from symfem.elements.bernardi_raugel import BernardiRaugel
+from symfem.elements.lagrange import Lagrange, VectorLagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import DotPointEvaluation, ListOfFunctionals, NormalIntegralMoment
 from symfem.functions import FunctionInput, VectorFunction
@@ -15,8 +17,6 @@
 from symfem.moments import make_integral_moment_dofs
 from symfem.piecewise_functions import PiecewiseFunction
 from symfem.references import NonDefaultReferenceError, Reference
-from symfem.elements.bernardi_raugel import BernardiRaugel
-from symfem.elements.lagrange import Lagrange, VectorLagrange
 
 __all__ = ["GuzmanNeilanFirstKind", "GuzmanNeilanSecondKind", "make_piecewise_lagrange"]
 
@@ -46,11 +46,11 @@ def __init__(self, reference: Reference, order: int):
             assert order == 2 and reference.name == "tetrahedron"
             dofs = br.dofs[:-3]
 
-        super().__init__(reference, order, poly, dofs, reference.tdim, reference.tdim)
+        super().__init__(reference, order, poly, dofs, reference.tdim, reference.tdim)  # type: ignore
 
     def _make_polyset_triangle(
         self, reference: Reference, order: int
-    ) -> typing.List[FunctionInput]:
+    ) -> typing.List[PiecewiseFunction]:
         """Make the polyset for a triangle.
 
         Args:
@@ -73,14 +73,13 @@ def _make_polyset_triangle(
         ]
 
         lagrange = VectorLagrange(reference, order)
-        basis: typing.List[FunctionInput] = [
+        return [
             PiecewiseFunction({i: p for i in sub_tris}, 2) for p in lagrange.get_polynomial_basis()
-        ] + [PiecewiseFunction(b, 2) for b in bubbles]
-        return basis
+        ] + [PiecewiseFunction(b, 2) for b in bubbles]  # type: ignore
 
     def _make_polyset_tetrahedron(
         self, reference: Reference, order: int
-    ) -> typing.List[FunctionInput]:
+    ) -> typing.List[PiecewiseFunction]:
         """Make the polyset for a tetrahedron.
 
         Args:
@@ -91,7 +90,7 @@ def _make_polyset_tetrahedron(
             The polynomial set
         """
         assert order in [1, 2]
-        from symfem.elements._guzman_neilan_tetrahedron import coeffs
+        from symfem.elements._guzman_neilan_tetrahedron import bubbles
 
         mid = reference.midpoint()
 
@@ -103,27 +102,9 @@ def _make_polyset_tetrahedron(
         ]
 
         lagrange = VectorLagrange(reference, order)
-        basis: typing.List[FunctionInput] = [
+        return [
             PiecewiseFunction({i: p for i in sub_tets}, 3) for p in lagrange.get_polynomial_basis()
-        ]
-
-        sub_basis = make_piecewise_lagrange(sub_tets, "tetrahedron", reference.tdim, True)
-        fs = BernardiRaugel(reference, 1).get_basis_functions()[-4:]
-        for c, f in zip(coeffs, fs):
-            assert isinstance(f, VectorFunction)
-            pieces: typing.Dict[SetOfPointsInput, FunctionInput] = {}
-            for tet in sub_tets:
-                function: typing.List[sympy.core.expr.Expr] = []
-                for j in range(reference.tdim):
-                    component = f[j]
-                    for k, b in zip(c, sub_basis):
-                        component -= k * b.pieces[tet][j]
-                    c_sym = component.as_sympy()
-                    assert isinstance(c_sym, sympy.core.expr.Expr)
-                    function.append(c_sym)
-                pieces[tet] = VectorFunction(tuple(function))
-            basis.append(PiecewiseFunction(pieces, 3))
-        return basis
+        ] + [PiecewiseFunction(b, 3) for b in bubbles]  # type: ignore
 
     @property
     def lagrange_subdegree(self) -> int:
@@ -227,11 +208,11 @@ def __init__(self, reference: Reference, order: int):
                 DotPointEvaluation(reference, mid, direction, entity=(tdim, 0), mapping="identity")
             )
 
-        super().__init__(reference, order, poly, dofs, tdim, tdim)
+        super().__init__(reference, order, poly, dofs, tdim, tdim)  # type: ignore
 
     def _make_polyset_triangle(
         self, reference: Reference, order: int
-    ) -> typing.List[FunctionInput]:
+    ) -> typing.List[PiecewiseFunction]:
         """Make the polyset for a triangle.
 
         Args:
@@ -243,7 +224,7 @@ def _make_polyset_triangle(
         """
         assert order == 1
 
-        from symfem.elements._guzman_neilan_triangle import coeffs
+        from symfem.elements._guzman_neilan_triangle import bubbles
 
         mid = reference.midpoint()
 
@@ -253,30 +234,14 @@ def _make_polyset_triangle(
             (reference.vertices[1], reference.vertices[2], mid),
         ]
 
-        basis: typing.List[FunctionInput] = []
-        basis += make_piecewise_lagrange(sub_tris, "triangle", order)
-
-        sub_basis = make_piecewise_lagrange(sub_tris, "triangle", reference.tdim, True)
-        fs = BernardiRaugel(reference, 1).get_basis_functions()[-3:]
-        for c, f in zip(coeffs, fs):
-            assert isinstance(f, VectorFunction)
-            pieces: typing.Dict[SetOfPointsInput, FunctionInput] = {}
-            for tri in sub_tris:
-                function: typing.List[sympy.core.expr.Expr] = []
-                for j in range(reference.tdim):
-                    component = f[j]
-                    for k, b in zip(c, sub_basis):
-                        component -= k * b.pieces[tri][j]
-                    c_sym = component.as_sympy()
-                    assert isinstance(c_sym, sympy.core.expr.Expr)
-                    function.append(c_sym)
-                pieces[tri] = VectorFunction(tuple(function))
-            basis.append(PiecewiseFunction(pieces, 2))
-        return basis
+        return make_piecewise_lagrange(sub_tris, "triangle", order) + [
+            PiecewiseFunction(b, 2)
+            for b in bubbles  # type: ignore
+        ]
 
     def _make_polyset_tetrahedron(
         self, reference: Reference, order: int
-    ) -> typing.List[FunctionInput]:
+    ) -> typing.List[PiecewiseFunction]:
         """Make the polyset for a tetrahedron.
 
         Args:
@@ -287,7 +252,7 @@ def _make_polyset_tetrahedron(
             The polynomial set
         """
         assert order in [1, 2]
-        from symfem.elements._guzman_neilan_tetrahedron import coeffs
+        from symfem.elements._guzman_neilan_tetrahedron import bubbles
 
         mid = reference.midpoint()
 
@@ -298,26 +263,10 @@ def _make_polyset_tetrahedron(
             (reference.vertices[1], reference.vertices[2], reference.vertices[3], mid),
         ]
 
-        basis: typing.List[FunctionInput] = []
-        basis += make_piecewise_lagrange(sub_tets, "tetrahedron", order)
-
-        sub_basis = make_piecewise_lagrange(sub_tets, "tetrahedron", reference.tdim, True)
-        fs = BernardiRaugel(reference, 1).get_basis_functions()[-4:]
-        for c, f in zip(coeffs, fs):
-            assert isinstance(f, VectorFunction)
-            pieces: typing.Dict[SetOfPointsInput, FunctionInput] = {}
-            for tet in sub_tets:
-                function: typing.List[sympy.core.expr.Expr] = []
-                for j in range(reference.tdim):
-                    component = f[j]
-                    for k, b in zip(c, sub_basis):
-                        component -= k * b.pieces[tet][j]
-                    c_sym = component.as_sympy()
-                    assert isinstance(c_sym, sympy.core.expr.Expr)
-                    function.append(c_sym)
-                pieces[tet] = VectorFunction(tuple(function))
-            basis.append(PiecewiseFunction(pieces, 3))
-        return basis
+        return make_piecewise_lagrange(sub_tets, "tetrahedron", order) + [
+            PiecewiseFunction(b, 3)
+            for b in bubbles  # type: ignore
+        ]
 
     @property
     def lagrange_subdegree(self) -> int:
diff --git a/symfem/elements/hhj.py b/symfem/elements/hhj.py
index c97dfb4..66fcdee 100644
--- a/symfem/elements/hhj.py
+++ b/symfem/elements/hhj.py
@@ -13,6 +13,7 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import IntegralMoment, ListOfFunctionals, NormalInnerProductIntegralMoment
 from symfem.functions import FunctionInput
@@ -20,7 +21,6 @@
 from symfem.moments import make_integral_moment_dofs
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["HellanHerrmannJohnson"]
 
diff --git a/symfem/elements/huang_zhang.py b/symfem/elements/huang_zhang.py
index 41fdede..230118b 100644
--- a/symfem/elements/huang_zhang.py
+++ b/symfem/elements/huang_zhang.py
@@ -6,6 +6,7 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     IntegralAgainst,
@@ -17,7 +18,6 @@
 from symfem.moments import make_integral_moment_dofs
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["HuangZhang"]
 
diff --git a/symfem/elements/mtw.py b/symfem/elements/mtw.py
index 147d149..63b4fd6 100644
--- a/symfem/elements/mtw.py
+++ b/symfem/elements/mtw.py
@@ -7,6 +7,8 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange
+from symfem.elements.nedelec import NedelecFirstKind
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     IntegralMoment,
@@ -19,8 +21,6 @@
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import Lagrange
-from symfem.elements.nedelec import NedelecFirstKind
 
 __all__ = ["MardalTaiWinther"]
 
diff --git a/symfem/elements/nedelec.py b/symfem/elements/nedelec.py
index 29502cf..1c0c52f 100644
--- a/symfem/elements/nedelec.py
+++ b/symfem/elements/nedelec.py
@@ -6,14 +6,14 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange, VectorLagrange
+from symfem.elements.rt import RaviartThomas
 from symfem.finite_element import CiarletElement
 from symfem.functionals import IntegralMoment, ListOfFunctionals, TangentIntegralMoment
 from symfem.functions import FunctionInput
 from symfem.moments import make_integral_moment_dofs
 from symfem.polynomials import Hcurl_polynomials, polynomial_set_vector
 from symfem.references import Reference
-from symfem.elements.lagrange import Lagrange, VectorLagrange
-from symfem.elements.rt import RaviartThomas
 
 __all__ = ["NedelecFirstKind", "NedelecSecondKind"]
 
diff --git a/symfem/elements/nedelec_prism.py b/symfem/elements/nedelec_prism.py
index 64f66c8..224fdde 100644
--- a/symfem/elements/nedelec_prism.py
+++ b/symfem/elements/nedelec_prism.py
@@ -2,6 +2,8 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange, VectorLagrange
+from symfem.elements.q import RaviartThomas as QRT
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     IntegralAgainst,
@@ -14,8 +16,6 @@
 from symfem.polynomials import Hcurl_polynomials, polynomial_set_1d, polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import Lagrange, VectorLagrange
-from symfem.elements.q import RaviartThomas as QRT
 
 __all__ = ["Nedelec"]
 
diff --git a/symfem/elements/regge.py b/symfem/elements/regge.py
index 47f04b8..e8db48d 100644
--- a/symfem/elements/regge.py
+++ b/symfem/elements/regge.py
@@ -12,6 +12,7 @@
 
 import sympy
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     InnerProductIntegralMoment,
@@ -25,7 +26,6 @@
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import Reference
 from symfem.symbols import t, x
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["Regge", "ReggeTP"]
 
diff --git a/symfem/elements/rt.py b/symfem/elements/rt.py
index a924b46..d90f06b 100644
--- a/symfem/elements/rt.py
+++ b/symfem/elements/rt.py
@@ -6,13 +6,13 @@
 
 import typing
 
+from symfem.elements.lagrange import Lagrange, VectorLagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import IntegralMoment, ListOfFunctionals, NormalIntegralMoment
 from symfem.functions import FunctionInput
 from symfem.moments import make_integral_moment_dofs
 from symfem.polynomials import Hdiv_polynomials, polynomial_set_vector
 from symfem.references import Reference
-from symfem.elements.lagrange import Lagrange, VectorLagrange
 
 __all__ = ["RaviartThomas"]
 
diff --git a/symfem/elements/serendipity.py b/symfem/elements/serendipity.py
index 453de53..cb744a7 100644
--- a/symfem/elements/serendipity.py
+++ b/symfem/elements/serendipity.py
@@ -6,6 +6,7 @@
 
 import typing
 
+from symfem.elements.dpc import DPC, VectorDPC
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     IntegralMoment,
@@ -24,7 +25,6 @@
     serendipity_set_1d,
 )
 from symfem.references import NonDefaultReferenceError, Reference
-from symfem.elements.dpc import DPC, VectorDPC
 
 __all__ = ["Serendipity", "SerendipityCurl", "SerendipityDiv"]
 
diff --git a/symfem/elements/taylor.py b/symfem/elements/taylor.py
index 4e67c4a..379fa40 100644
--- a/symfem/elements/taylor.py
+++ b/symfem/elements/taylor.py
@@ -3,13 +3,13 @@
 import typing
 from itertools import product
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import DerivativePointEvaluation, IntegralMoment, ListOfFunctionals
 from symfem.functions import FunctionInput
 from symfem.moments import make_integral_moment_dofs
 from symfem.polynomials import polynomial_set_1d
 from symfem.references import Reference
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["Taylor"]
 
diff --git a/symfem/elements/tnt.py b/symfem/elements/tnt.py
index 5fb4da2..2a56dfd 100644
--- a/symfem/elements/tnt.py
+++ b/symfem/elements/tnt.py
@@ -9,6 +9,7 @@
 
 import sympy
 
+from symfem.elements.q import Q
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     DerivativeIntegralMoment,
@@ -23,7 +24,6 @@
 from symfem.polynomials import orthogonal_basis, quolynomial_set_1d, quolynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import t, x
-from symfem.elements.q import Q
 
 __all__ = ["p", "b", "TNT", "TNTcurl", "TNTdiv"]
 
diff --git a/symfem/elements/transition.py b/symfem/elements/transition.py
index 39a3abd..0d69db6 100644
--- a/symfem/elements/transition.py
+++ b/symfem/elements/transition.py
@@ -5,6 +5,7 @@
 
 import sympy
 
+from symfem.elements.lagrange import Lagrange
 from symfem.finite_element import CiarletElement
 from symfem.functionals import ListOfFunctionals, PointEvaluation
 from symfem.functions import FunctionInput, ScalarFunction
@@ -12,7 +13,6 @@
 from symfem.quadrature import get_quadrature
 from symfem.references import Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import Lagrange
 
 __all__ = ["Transition"]
 
diff --git a/symfem/elements/trimmed_serendipity.py b/symfem/elements/trimmed_serendipity.py
index 5e050d0..16e0891 100644
--- a/symfem/elements/trimmed_serendipity.py
+++ b/symfem/elements/trimmed_serendipity.py
@@ -7,6 +7,7 @@
 
 import typing
 
+from symfem.elements.dpc import DPC, VectorDPC
 from symfem.finite_element import CiarletElement
 from symfem.functionals import (
     IntegralAgainst,
@@ -20,7 +21,6 @@
 from symfem.polynomials import polynomial_set_vector
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import t, x
-from symfem.elements.dpc import DPC, VectorDPC
 
 __all__ = ["TrimmedSerendipityHcurl", "TrimmedSerendipityHdiv"]
 
diff --git a/symfem/elements/vector_enriched_galerkin.py b/symfem/elements/vector_enriched_galerkin.py
index bdc6538..2476b0b 100644
--- a/symfem/elements/vector_enriched_galerkin.py
+++ b/symfem/elements/vector_enriched_galerkin.py
@@ -6,13 +6,13 @@
 
 import typing
 
+from symfem.elements.lagrange import VectorLagrange
+from symfem.elements.q import VectorQ
 from symfem.finite_element import CiarletElement, EnrichedElement
 from symfem.functionals import BaseFunctional, IntegralAgainst
 from symfem.functions import FunctionInput, VectorFunction
 from symfem.references import NonDefaultReferenceError, Reference
 from symfem.symbols import x
-from symfem.elements.lagrange import VectorLagrange
-from symfem.elements.q import VectorQ
 
 __all__ = ["Enrichment", "VectorEnrichedGalerkin"]
 
diff --git a/symfem/functions.py b/symfem/functions.py
index 121c941..3e9a135 100644
--- a/symfem/functions.py
+++ b/symfem/functions.py
@@ -8,7 +8,6 @@
 import sympy
 
 import symfem.references
-
 from symfem.geometry import PointType
 from symfem.symbols import AxisVariables, AxisVariablesNotSingle, t, x
 
diff --git a/symfem/piecewise_functions.py b/symfem/piecewise_functions.py
index 3d45b7c..974f1b1 100644
--- a/symfem/piecewise_functions.py
+++ b/symfem/piecewise_functions.py
@@ -7,7 +7,6 @@
 import sympy
 
 import symfem
-
 from symfem.functions import (
     AnyFunction,
     FunctionInput,
diff --git a/symfem/polynomials/lobatto.py b/symfem/polynomials/lobatto.py
index adc8dbd..10a1e34 100644
--- a/symfem/polynomials/lobatto.py
+++ b/symfem/polynomials/lobatto.py
@@ -3,9 +3,9 @@
 import typing
 
 from symfem.functions import ScalarFunction
-from symfem.symbols import x
 from symfem.polynomials.dual import l2_dual
 from symfem.polynomials.legendre import orthonormal_basis
+from symfem.symbols import x
 
 __all__: typing.List[str] = []
 
diff --git a/symfem/references.py b/symfem/references.py
index 7716b5f..d39fbce 100644
--- a/symfem/references.py
+++ b/symfem/references.py
@@ -8,7 +8,6 @@
 import sympy
 
 import symfem.functions
-
 from symfem.geometry import (
     PointType,
     PointTypeInput,
diff --git a/test/test_degrees.py b/test/test_degrees.py
index f93e075..6695274 100644
--- a/test/test_degrees.py
+++ b/test/test_degrees.py
@@ -1,12 +1,12 @@
 """Test every element."""
 
 import pytest
-
 import sympy
+
 from symfem import create_element
-from symfem.functions import VectorFunction, MatrixFunction
-from symfem.symbols import x
+from symfem.functions import MatrixFunction, VectorFunction
 from symfem.polynomials import polynomial_set, polynomial_set_1d
+from symfem.symbols import x
 
 from .utils import test_elements
 
diff --git a/test/test_guzman_neilan.py b/test/test_guzman_neilan.py
index 54c6452..88d1565 100644
--- a/test/test_guzman_neilan.py
+++ b/test/test_guzman_neilan.py
@@ -47,12 +47,12 @@ def test_guzman_neilan_triangle(order):
 def test_guzman_neilan_tetrahedron(order):
     e = symfem.create_element("tetrahedron", "Guzman-Neilan second kind", order)
 
-    mid = tuple(sympy.Rational(sum(i), len(i)) for i in zip(*e.reference.vertices))
     for p in e._basis[-4:]:
         for piece in p.pieces.values():
             float(piece.div().as_sympy().expand())
 
-        assert p.subs(x, mid) == (0, 0, 0)
+        for v in e.reference.vertices:
+            assert p.subs(x, v) == (0, 0, 0)
 
 
 @pytest.mark.parametrize("order", [1])