diff --git a/scripts/generate_gn.py b/scripts/generate_gn.py
index 39163399..f78d13f8 100644
--- a/scripts/generate_gn.py
+++ b/scripts/generate_gn.py
@@ -4,52 +4,172 @@
 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
 
 TESTING = "test" in sys.argv
 
-line_length = 100
-if TESTING:
-    folder = os.path.join(os.path.dirname(os.path.realpath(__file__)), "../_temp")
-else:
-    folder = os.path.join(os.path.dirname(os.path.realpath(__file__)), "../symfem/elements")
 
+def poly(reference, k):
+    """Generate the P^perp polynomial set."""
+    if k < 2:
+        if reference.name == "triangle":
+            return [
+                VectorFunction(
+                    [x[0] ** i0 * x[1] ** i1 if d2 == d else 0 for d2 in range(reference.tdim)]
+                )
+                for d in range(reference.tdim)
+                for i0 in range(k + 1)
+                for i1 in range(k + 1 - i0)
+            ]
+        else:
+            assert reference.name == "tetrahedron"
+            return [
+                VectorFunction(
+                    [
+                        x[0] ** i0 * x[1] ** i1 * x[2] ** i2 if d2 == d else 0
+                        for d2 in range(reference.tdim)
+                    ]
+                )
+                for d in range(reference.tdim)
+                for i0 in range(k + 1)
+                for i1 in range(k + 1 - i0)
+                for i2 in range(k + 1 - i0 - i1)
+            ]
+    if k == 2:
+        assert reference.name == "tetrahedron"
 
-def find_solution(mat, aim):
-    """Solve matrix-vector problem."""
-    A_data = [[float(j) for j in i] for i in mat]
-    b_data = [float(i) for i in aim]
-    A = np.asarray(A_data, dtype=np.float64)
-    b = np.asarray(b_data)
+        poly = [
+            VectorFunction(
+                [
+                    x[0] ** i0 * x[1] ** i1 * x[2] ** i2 if d2 == d else 0
+                    for d2 in range(reference.tdim)
+                ]
+            )
+            for d in range(reference.tdim)
+            for i0 in range(k + 1)
+            for i1 in range(k + 1 - i0)
+            for i2 in range(k + 1 - i0 - i1)
+        ]
 
-    res = np.linalg.solve(A, b)
+        poly[1] -= poly[0] / 4
+        poly[2] -= poly[0] / 10
+        poly[3] -= poly[0] / 4
+        poly[4] -= poly[0] / 20
+        poly[5] -= poly[0] / 10
+        poly[6] -= poly[0] / 4
+        poly[7] -= poly[0] / 20
+        poly[8] -= poly[0] / 20
+        poly[9] -= poly[0] / 10
+        poly = poly[1:]
 
-    fractions = []
-    for i in res:
-        frac = sympy.Rational(int(round(i * 162000)), 162000)
-        assert i < 10
-        assert np.isclose(float(frac), i)
-        fractions.append(frac)
+        poly[10] -= poly[9] / 4
+        poly[11] -= poly[9] / 10
+        poly[12] -= poly[9] / 4
+        poly[13] -= poly[9] / 20
+        poly[14] -= poly[9] / 10
+        poly[15] -= poly[9] / 4
+        poly[16] -= poly[9] / 20
+        poly[17] -= poly[9] / 20
+        poly[18] -= poly[9] / 10
+        poly = poly[:9] + poly[10:]
 
-    assert MatrixFunction(mat) @ VectorFunction(fractions) == VectorFunction(aim)
+        poly[19] -= poly[18] / 4
+        poly[20] -= poly[18] / 10
+        poly[21] -= poly[18] / 4
+        poly[22] -= poly[18] / 20
+        poly[23] -= poly[18] / 10
+        poly[24] -= poly[18] / 4
+        poly[25] -= poly[18] / 20
+        poly[26] -= poly[18] / 20
+        poly[27] -= poly[18] / 10
+        poly = poly[:18] + poly[19:]
 
-    return fractions
+        poly[1] -= poly[0] * 2 / 3
+        poly[2] += poly[0] / 3
+        poly[3] -= poly[0] / 9
+        poly[4] += poly[0] * 2 / 9
+        poly[5] += poly[0] / 3
+        poly[6] -= poly[0] / 9
+        poly[7] += poly[0] / 9
+        poly[8] += poly[0] * 2 / 9
+        poly[18] -= poly[0] / 3
+        poly[19] -= poly[0] * 2 / 9
+        poly[20] -= poly[0] / 3
+        poly[21] -= poly[0] / 9
+        poly[22] -= poly[0] * 2 / 9
+        poly[23] += poly[0]
+        poly[24] += poly[0] / 9
+        poly[25] += poly[0] / 9
+        poly[26] += poly[0] * 2 / 3
+        poly = poly[1:]
 
+        poly[9] -= poly[8] * 2 / 3
+        poly[10] += poly[8] / 3
+        poly[11] -= poly[8] / 9
+        poly[12] += poly[8] * 2 / 9
+        poly[13] += poly[8] / 3
+        poly[14] -= poly[8] / 9
+        poly[15] += poly[8] / 9
+        poly[16] += poly[8] * 2 / 9
+        poly[17] -= poly[8] / 3
+        poly[18] -= poly[8] * 2 / 9
+        poly[19] += poly[8]
+        poly[20] += poly[8] / 9
+        poly[21] += poly[8] * 2 / 3
+        poly[22] -= poly[8] / 3
+        poly[23] -= poly[8] / 9
+        poly[24] += poly[8] / 9
+        poly[25] -= poly[8] * 2 / 9
+        poly = poly[:8] + poly[9:]
 
-for ref in ["triangle", "tetrahedron"]:
-    # TODO: work out why tetrahedron fails on github but passes locally
-    if TESTING and ref == "tetrahedron":
-        break
+        poly[2] -= poly[1] / 6
+        poly[3] -= poly[1] * 2 / 3
+        poly[4] += poly[1] / 2
+        poly[5] += poly[1] / 12
+        poly[6] -= poly[1] / 12
+        poly[7] += poly[1] / 3
+        poly[9] -= poly[1] / 2
+        poly[10] -= poly[1] / 12
+        poly[11] -= poly[1] / 3
+        poly[12] += poly[1]
+        poly[13] += poly[1] / 6
+        poly[14] += poly[1] / 12
+        poly[15] += poly[1] * 2 / 3
+        poly[18] -= poly[1] / 2
+        poly[19] -= poly[1] / 12
+        poly[20] -= poly[1] / 3
+        poly[21] += poly[1] / 2
+        poly[22] += poly[1] / 12
+        poly[24] += poly[1] / 3
+        poly = poly[:1] + poly[2:]
 
+        return poly
+
+    raise NotImplementedError()
+
+
+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)
+    return list(solution)
+
+
+if TESTING:
+    folder = os.path.join(os.path.dirname(os.path.realpath(__file__)), "../_temp")
+else:
+    folder = os.path.join(os.path.dirname(os.path.realpath(__file__)), "../symfem/elements")
+
+for ref in ["triangle", "tetrahedron"]:
     reference = symfem.create_reference(ref)
     br = symfem.create_element(ref, "Bernardi-Raugel", 1)
-    mid = tuple(sympy.Rational(sum(i), len(i)) for i in zip(*reference.vertices))
+    mid = reference.midpoint()
 
     sub_cells: typing.List[symfem.geometry.SetOfPoints] = []
 
@@ -62,6 +182,16 @@ def find_solution(mat, aim):
         ]
         xx, yy, zz = x
         terms = [1, xx, yy]
+
+        lamb = PiecewiseFunction(
+            {sub_cells[0]: 3 * x[1], sub_cells[1]: 3 * x[0], sub_cells[2]: 3 * (1 - x[0] - x[1])}, 2
+        )
+
+        for c, f in lamb.pieces.items():
+            assert f.subs(x, c[0]) == 0
+            assert f.subs(x, c[1]) == 0
+            assert f.subs(x, c[2]) == 1
+
     if ref == "tetrahedron":
         fs = br.get_basis_functions()[-4:]
         sub_cells = [
@@ -73,16 +203,37 @@ def find_solution(mat, aim):
         xx, yy, zz = x
         terms = [1, xx, yy, zz, xx**2, yy**2, zz**2, xx * yy, xx * zz, yy * zz]
 
-    sub_basis = make_piecewise_lagrange(sub_cells, ref, br.reference.tdim, True)
+        lamb = PiecewiseFunction(
+            {
+                sub_cells[0]: 4 * x[2],
+                sub_cells[1]: 4 * x[1],
+                sub_cells[2]: 4 * x[0],
+                sub_cells[3]: 4 * (1 - x[0] - x[1] - x[2]),
+            },
+            3,
+        )
+
+        for c, f in lamb.pieces.items():
+            assert f.subs(x, c[0]) == 0
+            assert f.subs(x, c[1]) == 0
+            assert f.subs(x, c[2]) == 0
+            assert f.subs(x, c[3]) == 1
+
+    sub_basis = [
+        p * lamb**j
+        for j in range(1, reference.tdim + 1)
+        for p in poly(reference, reference.tdim - j)
+    ]
 
     filename = os.path.join(folder, f"_guzman_neilan_{ref}.py")
-    output = '"""Values for Guzman-Neilan element."""\n' "\n" "import sympy\n" "\n" "coeffs = [\n"
+    output = '"""Values for Guzman-Neilan element."""\n\n'
+    output += "import sympy\n\nbubbles = [\n"
 
     for f in fs:
         assert isinstance(f, VectorFunction)
-        output += "    [\n"
         integrand = f.div().subs(x, t)
         fun_s = (f.div() - integrand.integral(reference) / reference.volume()).as_sympy()
+
         assert isinstance(fun_s, sympy.core.expr.Expr)
         fun = fun_s.as_coefficients_dict()
 
@@ -104,53 +255,45 @@ def find_solution(mat, aim):
                     assert term == 0 or term in terms
                 for term in terms:
                     if i < len(mat):
-                        mat[i].append(d[term] if term in d else 0)
+                        mat[i].append(d[term] if term in d else 0)  # type: ignore
                     i += 1
 
-        if ref == "triangle":
-            mat = mat[:-1]
-            aim = aim[:-1]
-            fractions = find_solution(mat, aim)
-        if ref == "tetrahedron":
-            for i in range(3):
-                row: typing.List[symfem.functions.ScalarFunction] = [sympy.Integer(0)] * 45
-                row[i] = sympy.Integer(1)
-                mat.append(row)
-            subf = f.subs(x, mid)
-            assert isinstance(subf, VectorFunction)
-            aim += [i for i in subf]
-            aim += [0, 0]
-
-            fractions = None
-            for n in range(3, 45):
-                for m in range(n + 1, 45):
-                    mat2 = [i for i in mat]
-                    for i in [n, m]:
-                        row = [sympy.Integer(0)] * 45
-                        row[i] = sympy.Integer(1)
-                        mat2.append(row)
-                    try:
-                        fractions = find_solution(mat2, aim)
-                        break
-                    except AssertionError:
-                        pass
-                else:
-                    continue
-                break
-            assert fractions is not None
-
-        line = " " * 7
-        for frac in fractions:
-            if frac.denominator == 1:
-                next = f" {frac.numerator},"
-            else:
-                next = f" sympy.Rational({frac.numerator}, {frac.denominator}),"
-            if len(line + next) > line_length:
-                output += line + "\n"
-                line = " " * 7 + next
-            else:
-                line += next
-        output += line + "\n    ],\n"
+        coeffs = find_solution(mat, aim)
+        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 += "        "
+            output += (
+                "("
+                + ", ".join(
+                    [
+                        "("
+                        + ", ".join(
+                            [
+                                f"{c}" if isinstance(c, sympy.Integer) else f"sympy.S('{c}')"
+                                for c in p
+                            ]
+                        )
+                        + ")"
+                        for p in cell
+                    ]
+                )
+                + ")"
+            )
+            output += ": (\n"
+            output += ",\n".join(
+                [
+                    "            sympy.S('" + f"{c.as_sympy().expand()}".replace(" ", "") + "')"  # type:ignore
+                    for c in f
+                ]
+            )
+            output += "),\n"
+        output += "    },\n"
 
     output += "]\n"
 
diff --git a/symfem/basis_functions.py b/symfem/basis_functions.py
index 42859fe7..447aa53b 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 aaa5fa78..937499ed 100644
--- a/symfem/elements/_guzman_neilan_tetrahedron.py
+++ b/symfem/elements/_guzman_neilan_tetrahedron.py
@@ -2,193 +2,129 @@
 
 import sympy
 
-coeffs = [
-    [
-        sympy.Rational(5, 8),
-        sympy.Rational(5, 8),
-        sympy.Rational(5, 8),
-        0,
-        0,
-        sympy.Rational(125, 72),
-        sympy.Rational(17, 36),
-        sympy.Rational(17, 36),
-        sympy.Rational(71, 18),
-        sympy.Rational(-6343, 81000),
-        sympy.Rational(-16157, 81000),
-        sympy.Rational(325, 216),
-        sympy.Rational(6391, 20250),
-        sympy.Rational(742, 10125),
-        sympy.Rational(58, 27),
-        sympy.Rational(-7879, 10125),
-        sympy.Rational(325, 216),
-        sympy.Rational(-100093, 81000),
-        sympy.Rational(-43907, 40500),
-        sympy.Rational(58, 27),
-        sympy.Rational(-20242, 10125),
-        sympy.Rational(325, 216),
-        sympy.Rational(-26609, 40500),
-        sympy.Rational(-109907, 81000),
-        sympy.Rational(58, 27),
-        sympy.Rational(-34093, 40500),
-        sympy.Rational(-45391, 20250),
-        sympy.Rational(-13471, 54000),
-        sympy.Rational(15971, 54000),
-        sympy.Rational(-25, 27),
-        sympy.Rational(-7879, 20250),
-        sympy.Rational(-545, 144),
-        sympy.Rational(10133, 40500),
-        sympy.Rational(-545, 144),
-        sympy.Rational(-26609, 81000),
-        sympy.Rational(15359, 81000),
-        sympy.Rational(39149, 13500),
-        sympy.Rational(-25, 432),
-        sympy.Rational(-9631, 3375),
-        sympy.Rational(-6343, 162000),
-        sympy.Rational(-16157, 162000),
-        sympy.Rational(145, 36),
-        sympy.Rational(-25, 432),
-        sympy.Rational(63577, 27000),
-        sympy.Rational(-62327, 27000),
-    ],
-    [
-        sympy.Rational(15, 8),
-        0,
-        0,
-        0,
-        sympy.Rational(125, 72),
-        0,
-        sympy.Rational(1, 12),
-        sympy.Rational(113, 36),
-        sympy.Rational(-1, 3),
-        sympy.Rational(-10921, 8100),
-        sympy.Rational(27467, 16200),
-        sympy.Rational(-25, 18),
-        sympy.Rational(-21167, 8100),
-        sympy.Rational(24767, 8100),
-        sympy.Rational(-16, 9),
-        sympy.Rational(25217, 16200),
-        sympy.Rational(-25, 18),
-        sympy.Rational(8533, 16200),
-        sympy.Rational(6473, 2025),
-        sympy.Rational(-16, 9),
-        sympy.Rational(5833, 8100),
-        sympy.Rational(125, 72),
-        sympy.Rational(-8273, 4050),
-        sympy.Rational(13967, 16200),
-        sympy.Rational(44, 9),
-        sympy.Rational(-8948, 2025),
-        sympy.Rational(11267, 8100),
-        sympy.Rational(9671, 1800),
-        sympy.Rational(-2449, 450),
-        sympy.Rational(25, 36),
-        sympy.Rational(25217, 32400),
-        sympy.Rational(235, 48),
-        sympy.Rational(8533, 32400),
-        sympy.Rational(-235, 144),
-        sympy.Rational(-4967, 32400),
-        sympy.Rational(13967, 32400),
-        sympy.Rational(-9239, 1200),
-        sympy.Rational(-25, 144),
-        sympy.Rational(27467, 3600),
-        sympy.Rational(-10921, 16200),
-        sympy.Rational(6949, 4050),
-        sympy.Rational(-35, 12),
-        sympy.Rational(-145, 144),
-        sympy.Rational(12967, 3600),
-        sympy.Rational(-663, 400),
-    ],
-    [
-        0,
-        sympy.Rational(-15, 8),
-        0,
-        0,
-        sympy.Rational(-61387, 81000),
-        sympy.Rational(-39619, 40500),
-        sympy.Rational(1, 3),
-        sympy.Rational(-32381, 20250),
-        sympy.Rational(-32869, 20250),
-        sympy.Rational(-37, 18),
-        sympy.Rational(59137, 81000),
-        sympy.Rational(25, 18),
-        sympy.Rational(-34, 9),
-        sympy.Rational(27881, 20250),
-        sympy.Rational(16, 9),
-        sympy.Rational(2, 3),
-        sympy.Rational(-125, 72),
-        sympy.Rational(-9881, 40500),
-        sympy.Rational(5, 3),
-        sympy.Rational(-44, 9),
-        sympy.Rational(-3131, 20250),
-        sympy.Rational(25, 18),
-        sympy.Rational(-13, 72),
-        sympy.Rational(-1, 6),
-        sympy.Rational(16, 9),
-        sympy.Rational(-4, 9),
-        0,
-        sympy.Rational(-3, 4),
-        sympy.Rational(59, 72),
-        sympy.Rational(-16631, 81000),
-        sympy.Rational(1, 3),
-        sympy.Rational(-8003, 4500),
-        sympy.Rational(-11, 18),
-        sympy.Rational(35, 12),
-        sympy.Rational(-19003, 40500),
-        sympy.Rational(-46369, 81000),
-        sympy.Rational(-47, 9),
-        sympy.Rational(14032, 10125),
-        sympy.Rational(59, 18),
-        sympy.Rational(-37, 36),
-        sympy.Rational(-1, 72),
-        sympy.Rational(-13369, 9000),
-        sympy.Rational(-25, 36),
-        sympy.Rational(1844, 375),
-        sympy.Rational(-43631, 9000),
-    ],
-    [
-        0,
-        0,
-        sympy.Rational(15, 8),
-        sympy.Rational(108281, 40500),
-        0,
-        sympy.Rational(-75937, 81000),
-        sympy.Rational(101531, 20250),
-        sympy.Rational(-1, 3),
-        sympy.Rational(-36281, 20250),
-        sympy.Rational(-132187, 81000),
-        sympy.Rational(-35, 72),
-        sympy.Rational(125, 72),
-        sympy.Rational(-145687, 40500),
-        sympy.Rational(-47, 36),
-        sympy.Rational(44, 9),
-        sympy.Rational(25, 72),
-        sympy.Rational(-25, 18),
-        0,
-        sympy.Rational(13, 36),
-        sympy.Rational(-16, 9),
-        sympy.Rational(1, 12),
-        sympy.Rational(-25, 18),
-        sympy.Rational(15, 8),
-        sympy.Rational(23203, 20250),
-        sympy.Rational(-16, 9),
-        sympy.Rational(41, 12),
-        sympy.Rational(96187, 40500),
-        sympy.Rational(-355, 144),
-        sympy.Rational(635, 144),
-        sympy.Rational(-21797, 40500),
-        sympy.Rational(244687, 162000),
-        sympy.Rational(-35, 12),
-        sympy.Rational(-75937, 162000),
-        sympy.Rational(82031, 9000),
-        sympy.Rational(15, 16),
-        sympy.Rational(5, 48),
-        sympy.Rational(-64687, 18000),
-        sympy.Rational(25, 36),
-        sympy.Rational(63437, 18000),
-        sympy.Rational(25, 48),
-        sympy.Rational(-35, 144),
-        sympy.Rational(-2194, 375),
-        sympy.Rational(-52031, 81000),
-        sympy.Rational(25, 16),
-        sympy.Rational(-235, 144),
-    ],
+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(
+                "-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-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 58e8da0f..4ce6d5a2 100644
--- a/symfem/elements/_guzman_neilan_triangle.py
+++ b/symfem/elements/_guzman_neilan_triangle.py
@@ -2,35 +2,44 @@
 
 import sympy
 
-coeffs = [
-    [
-        sympy.Rational(1, 6),
-        sympy.Rational(1, 6),
-        sympy.Rational(5, 24),
-        sympy.Rational(-1, 24),
-        sympy.Rational(5, 24),
-        sympy.Rational(5, 24),
-        sympy.Rational(-1, 24),
-        sympy.Rational(5, 24),
-    ],
-    [
-        sympy.Rational(1, 3),
-        sympy.Rational(-2, 3),
-        sympy.Rational(5, 12),
-        sympy.Rational(-1, 3),
-        sympy.Rational(-1, 12),
-        sympy.Rational(-1, 3),
-        sympy.Rational(5, 12),
-        sympy.Rational(-5, 6),
-    ],
-    [
-        sympy.Rational(2, 3),
-        sympy.Rational(-1, 3),
-        sympy.Rational(5, 6),
-        sympy.Rational(-5, 12),
-        sympy.Rational(1, 3),
-        sympy.Rational(1, 12),
-        sympy.Rational(1, 3),
-        sympy.Rational(-5, 12),
-    ],
+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*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"),
+        ),
+    },
 ]
diff --git a/symfem/elements/abf.py b/symfem/elements/abf.py
index 864e7a4d..12c014c3 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 d34c6d18..b439fe52 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 64774597..d50b084f 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 7d0cdd57..bc9d76a7 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 a69589e4..ce611f01 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 05d0f0e1..3b71ba23 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 6f6ece6f..bdb3f9f5 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"]
 
@@ -132,4 +132,4 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
     max_order = {"triangle": 1, "tetrahedron": 2}
     continuity = "L2"
     value_type = "vector"
-    last_updated = "2024.10"
+    last_updated = "2024.10.1"
diff --git a/symfem/elements/bubble.py b/symfem/elements/bubble.py
index 28cf61c2..5fd2f540 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 96945514..97d7905d 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 4c704c92..af1e9fab 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 51b30fa4..0c34ae91 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 1447e38f..548af66f 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:
@@ -62,7 +62,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()
 
@@ -73,31 +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()
-        ]
-
-        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
+        ] + [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:
@@ -108,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()
 
@@ -120,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:
@@ -164,7 +128,7 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
     max_order = {"triangle": 1, "tetrahedron": 2}
     continuity = "L2"
     value_type = "vector macro"
-    last_updated = "2024.10.2"
+    last_updated = "2024.10.3"
 
 
 class GuzmanNeilanSecondKind(CiarletElement):
@@ -243,11 +207,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:
@@ -259,7 +223,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()
 
@@ -269,30 +233,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)  # type: ignore
+            for b in bubbles
+        ]
 
     def _make_polyset_tetrahedron(
         self, reference: Reference, order: int
-    ) -> typing.List[FunctionInput]:
+    ) -> typing.List[PiecewiseFunction]:
         """Make the polyset for a tetrahedron.
 
         Args:
@@ -303,7 +251,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()
 
@@ -314,26 +262,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)  # type: ignore
+            for b in bubbles
+        ]
 
     @property
     def lagrange_subdegree(self) -> int:
@@ -357,7 +289,7 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
     max_order = {"triangle": 1, "tetrahedron": 2}
     continuity = "L2"
     value_type = "vector macro"
-    last_updated = "2024.10"
+    last_updated = "2024.10.1"
 
 
 def make_piecewise_lagrange(
diff --git a/symfem/elements/hhj.py b/symfem/elements/hhj.py
index c97dfb46..66fcdee0 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 41fdeded..230118b6 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 147d1498..63b4fd65 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 29502cf0..1c0c52f5 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 64f66c8e..224fddef 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 47f04b8d..e8db48d1 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 a924b463..d90f06bb 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 453de53a..cb744a7e 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 4e67c4af..379fa404 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 5fb4da20..2a56dfd7 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 39a3abd0..0d69db63 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 5e050d07..16e08919 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 bdc65384..2476b0bb 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 121c941d..3e9a1357 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 3d45b7c9..974f1b14 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 adc8dbd1..10a1e34c 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 7716b5f2..d39fbce8 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 f93e075a..6695274a 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 373d81f0..cdda2d2c 100644
--- a/test/test_guzman_neilan.py
+++ b/test/test_guzman_neilan.py
@@ -10,25 +10,55 @@
 from symfem.symbols import x
 
 
-@pytest.mark.parametrize("order", [1])
-def test_guzman_neilan_triangle(order):
-    e = symfem.create_element("triangle", "Guzman-Neilan second kind", order)
+def test_triangle_bubbles():
+    from symfem.elements._guzman_neilan_triangle import bubbles
 
-    for p in e._basis[-3:]:
-        for piece in p.pieces.values():
-            float(piece.div().as_sympy().expand())
+    for b in bubbles:
+        div = None
+        for part in b.values():
+            value = (part[0].diff(x[0]) + part[1].diff(x[1])).expand()
+            if div is None:
+                div = value
+            assert div == value
 
 
-@pytest.mark.parametrize("order", [1, 2])
-def test_guzman_neilan_tetrahedron(order):
-    e = symfem.create_element("tetrahedron", "Guzman-Neilan second kind", order)
+def test_tetrahedron_bubbles():
+    from symfem.elements._guzman_neilan_tetrahedron import bubbles
+
+    for b in bubbles:
+        div = None
+        for part in b.values():
+            value = (part[0].diff(x[0]) + part[1].diff(x[1]) + part[2].diff(x[2])).expand()
+            if div is None:
+                div = value
+            assert div == value
 
-    mid = tuple(sympy.Rational(sum(i), len(i)) for i in zip(*e.reference.vertices))
-    for p in e._basis[-4:]:
+
+@pytest.mark.parametrize("cell,order", [("triangle", 1), ("tetrahedron", 1), ("tetrahedron", 2)])
+def test_guzman_neilan_triangle(cell, order):
+    e = symfem.create_element(cell, "Guzman-Neilan second kind", order)
+
+    if cell == "triangle":
+        nb = 3
+    elif cell == "tetrahedron":
+        nb = 4
+    else:
+        raise ValueError(f"Unsupported cell: {cell}")
+
+    mid = e.reference.midpoint()
+
+    for p in e._basis[-nb:]:
         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) == tuple(0 for _ in range(e.reference.tdim))
+
+        value = None
+        for piece in p.pieces.values():
+            if value is None:
+                value = piece.subs(x, mid)
+            assert value == piece.subs(x, mid)
 
 
 @pytest.mark.parametrize("order", [1])