From 103843d716d35480e22af4bfa614b558beb739ed Mon Sep 17 00:00:00 2001
From: Matthew Scroggs <matthew.w.scroggs@gmail.com>
Date: Fri, 18 Oct 2024 08:09:04 +0100
Subject: [PATCH] workgin on correcting G-N bubbles

---
 scripts/generate_gn_new.py                    | 253 ++++++++++++++++++
 symfem/elements/_guzman_neilan_tetrahedron.py | 194 --------------
 symfem/elements/_guzman_neilan_triangle.py    |  65 ++---
 symfem/elements/guzman_neilan.py              |  23 +-
 test/test_guzman_neilan.py                    |  12 +
 5 files changed, 303 insertions(+), 244 deletions(-)
 create mode 100644 scripts/generate_gn_new.py
 delete mode 100644 symfem/elements/_guzman_neilan_tetrahedron.py

diff --git a/scripts/generate_gn_new.py b/scripts/generate_gn_new.py
new file mode 100644
index 00000000..1302bd78
--- /dev/null
+++ b/scripts/generate_gn_new.py
@@ -0,0 +1,253 @@
+"""Generate coefficients for Guzman-Neilan element."""
+
+import os
+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.piecewise_functions import PiecewiseFunction
+from symfem.symbols import t, x
+
+TESTING = "test" in sys.argv
+
+
+def poly(reference, k):
+    if k < 2:
+        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)
+        ]
+    if k == 2:
+        assert reference.name == "tetrahedron"
+
+        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)
+        ]
+
+        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:]
+
+        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:]
+
+        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:]
+
+        #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]
+        #poly[20] -= poly[0]
+        #poly[21] -= poly[0]
+        #poly[22] += poly[0]
+        #poly[23] += poly[0]
+        #poly[24] += poly[0]
+        #poly[25] += poly[0]
+        #poly[26] += poly[0]
+
+        ned1 = symfem.create_element("tetrahedron", "Nedelec", 1).get_polynomial_basis()
+
+        for n, p in enumerate(poly):
+            print(n, [(p.dot(f)).integral(reference) for f in ned1])
+
+        return poly
+
+    raise NotImplementedError()
+
+
+def find_solution(mat, aim):
+    if len(mat) == 40:
+        mat = mat[:29] + mat[30:]
+        mat = mat[:18] + mat[19:]
+        mat = mat[:7] + mat[8:]
+
+    s_mat = sympy.Matrix(mat)
+    print(s_mat.shape)
+    if False:
+        def filter(mat, r, c):
+            return [row[:c] + row[c+1:] for row in mat[:r] + mat[r+1:]]
+        coords = {25: 1}
+        #mat = filter(mat, 25, 1)
+        #mat = filter(mat, 1, 8)
+        #mat = filter(mat, 1, 9)
+        #mat = filter(mat, 2, 8)
+        #mat = filter(mat, 2, 8)
+        print()
+        for n, row in enumerate(mat):
+            print((f" {n}")[-2:], "".join([" " if r == 0 else "*" for r in row]))
+
+        print(mat[7])
+        print(mat[18])
+        print(mat[29])
+
+        # from IPython import embed; embed()
+    solution = s_mat[:s_mat.shape[1], :].inv() @ sympy.Matrix(aim[:s_mat.shape[1]])
+
+    assert s_mat @ solution == sympy.Matrix(aim)
+    return list(solution)
+
+
+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")
+
+for ref in ["triangle", "tetrahedron"]:
+    reference = symfem.create_reference(ref)
+    br = symfem.create_element(ref, "Bernardi-Raugel", 1)
+    mid = reference.midpoint()
+
+    sub_cells: typing.List[symfem.geometry.SetOfPoints] = []
+
+    if ref == "triangle":
+        fs = br.get_basis_functions()[-3:]
+        sub_cells = [
+            (reference.vertices[0], reference.vertices[1], mid),
+            (reference.vertices[0], reference.vertices[2], mid),
+            (reference.vertices[1], reference.vertices[2], mid),
+        ]
+        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 = [
+            (reference.vertices[0], reference.vertices[1], reference.vertices[2], mid),
+            (reference.vertices[0], reference.vertices[1], reference.vertices[3], mid),
+            (reference.vertices[0], reference.vertices[2], reference.vertices[3], mid),
+            (reference.vertices[1], reference.vertices[2], reference.vertices[3], mid),
+        ]
+        xx, yy, zz = x
+        terms = [1, xx, yy, zz, xx**2, yy**2, zz**2, xx * yy, xx * zz, yy * zz]
+
+        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'
+    output += "import sympy\n\nbubbles = [\n"
+
+    for f in fs:
+        assert isinstance(f, VectorFunction)
+        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()
+
+        for term in fun:
+            assert term in terms
+        aim = [fun[term] if term in fun else 0 for term in terms] * (br.reference.tdim + 1)
+
+        mat: typing.List[typing.List[symfem.functions.ScalarFunction]] = [
+            [] for t in terms for p in sub_basis[0].pieces
+        ]
+        for b in sub_basis:
+            i = 0
+            for p in b.pieces.values():
+                assert isinstance(p, VectorFunction)
+                d_s = p.div().as_sympy()
+                assert isinstance(d_s, sympy.core.expr.Expr)
+                d = d_s.expand().as_coefficients_dict()
+                for term in d:
+                    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)  # type: ignore
+                    i += 1
+
+        coeffs = find_solution(mat, aim)
+        bubble = f
+        for i, j in zip(coeffs, sub_basis):
+            bubble -= i * j
+
+        output += "    {\n"
+        for cell, f in bubble.pieces.items():
+            output += "        "
+            output += "(" + ", ".join(["(" + ", ".join([
+                f"{c}" if isinstance(c, sympy.Integer) else "sympy.S('" + f"{c}".replace(" ", "") + "')" for c in p
+            ]) + ")" for p in cell]) + ")"
+            output += ": (\n"
+            output += ",\n".join(["            sympy.S('" + f"{c.as_sympy().expand()}".replace(" ", "") + "')" for c in f])
+            output += "),\n"
+        output += "    },\n"
+
+    output += "]\n"
+
+    with open(filename, "w") as ff:
+        ff.write(output)
diff --git a/symfem/elements/_guzman_neilan_tetrahedron.py b/symfem/elements/_guzman_neilan_tetrahedron.py
deleted file mode 100644
index aaa5fa78..00000000
--- a/symfem/elements/_guzman_neilan_tetrahedron.py
+++ /dev/null
@@ -1,194 +0,0 @@
-"""Values for Guzman-Neilan element."""
-
-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),
-    ],
-]
diff --git a/symfem/elements/_guzman_neilan_triangle.py b/symfem/elements/_guzman_neilan_triangle.py
index 58e8da0f..6246c4f1 100644
--- a/symfem/elements/_guzman_neilan_triangle.py
+++ b/symfem/elements/_guzman_neilan_triangle.py
@@ -2,35 +2,38 @@
 
 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/guzman_neilan.py b/symfem/elements/guzman_neilan.py
index 1447e38f..b6e3b7c5 100644
--- a/symfem/elements/guzman_neilan.py
+++ b/symfem/elements/guzman_neilan.py
@@ -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()
 
@@ -75,24 +75,7 @@ def _make_polyset_triangle(
         lagrange = VectorLagrange(reference, order)
         basis: typing.List[FunctionInput] = [
             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))
+        ] + [PiecewiseFunction(b, 2) for b in bubbles]
         return basis
 
     def _make_polyset_tetrahedron(
@@ -165,6 +148,7 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
     continuity = "L2"
     value_type = "vector macro"
     last_updated = "2024.10.2"
+    cache = False
 
 
 class GuzmanNeilanSecondKind(CiarletElement):
@@ -358,6 +342,7 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
     continuity = "L2"
     value_type = "vector macro"
     last_updated = "2024.10"
+    cache = False
 
 
 def make_piecewise_lagrange(
diff --git a/test/test_guzman_neilan.py b/test/test_guzman_neilan.py
index 373d81f0..a90fd70d 100644
--- a/test/test_guzman_neilan.py
+++ b/test/test_guzman_neilan.py
@@ -10,6 +10,18 @@
 from symfem.symbols import x
 
 
+def test_triangle_bubbles():
+    from symfem.elements._guzman_neilan_triangle import bubbles
+
+    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])
 def test_guzman_neilan_triangle(order):
     e = symfem.create_element("triangle", "Guzman-Neilan second kind", order)