Test that basis functions are independent (#84)

* add test that the basis functions are linearly independent

* adjust parameters for direct serendipity

* only run independence test for DirectElements, as Ciarlet elements are guaranteed to be independent

Author: mscroggs

symfem/elements/direct_serendipity.py

@@ -23,15 +23,15 @@ def __init__(self, reference, order, variant):
         # Functions on edges
         if order >= 2:
             alpha_h = 1
-            beta_h = 1
+            beta_h = 2
             gamma_h = 1
-            xi_h = 1
+            xi_h = 2
             eta_h = 1
 
             alpha_v = 1
-            beta_v = 1
+            beta_v = 2
             gamma_v = 1
-            xi_v = 1
+            xi_v = 2
             eta_v = 1
 
             lambda_h = alpha_h * (1 - x[1]) + beta_h * x[1] + gamma_h

test/test_elements.py

@@ -2,7 +2,7 @@
 import sympy
 import symfem
 from symfem import create_element
-from symfem.core.finite_element import CiarletElement
+from symfem.core.finite_element import CiarletElement, DirectElement
 from symfem.core.symbolic import subs, x, PiecewiseFunction
 from symfem.core.vectors import vsub
 from utils import test_elements, all_symequal
@@ -20,6 +20,67 @@ def test_all_tested():
                 raise ValueError(f"{e.names[0]} on a {r} is not tested")
 
 
+@pytest.mark.parametrize(
+    ("cell_type", "element_type", "order", "variant"),
+    [[reference, element, order, variant]
+     for reference, i in test_elements.items() for element, j in i.items()
+     for variant, k in j.items() for order in k])
+def test_independence(
+    elements_to_test, cells_to_test, cell_type, element_type, order, variant,
+    speed
+):
+    """Test that DirectElements have independent basis functions."""
+    if elements_to_test != "ALL" and element_type not in elements_to_test:
+        pytest.skip()
+    if cells_to_test != "ALL" and cell_type not in cells_to_test:
+        pytest.skip()
+    if speed == "fast":
+        if order > 2:
+            pytest.skip()
+        if order == 2 and cell_type in ["tetrahedron", "hexahedron", "prism", "pyramid"]:
+            pytest.skip()
+
+    space = create_element(cell_type, element_type, order, variant)
+
+    # Only run this test for DirectElements
+    if not isinstance(space, DirectElement):
+        pytest.skip()
+
+    basis = space.get_basis_functions()
+    all_terms = set()
+
+    try:
+        basis[0].as_coefficients_dict()
+        scalar = True
+    except AttributeError:
+        scalar = False
+
+    if scalar:
+        for f in basis:
+            for term in f.as_coefficients_dict():
+                all_terms.add(term)
+        mat = [[0 for i in all_terms] for j in basis]
+        for i, t in enumerate(all_terms):
+            for j, f in enumerate(basis):
+                fd = f.as_coefficients_dict()
+                if t in fd:
+                    mat[j][i] = fd[t]
+    else:
+        for f in basis:
+            for fi, fpart in enumerate(f):
+                for term in fpart.as_coefficients_dict():
+                    all_terms.add((fi, term))
+        mat = [[0 for i in all_terms] for j in basis]
+        for i, (fi, t) in enumerate(all_terms):
+            for j, f in enumerate(basis):
+                fd = f[fi].as_coefficients_dict()
+                if t in fd:
+                    mat[j][i] = fd[t]
+    mat = sympy.Matrix(mat)
+
+    assert mat.rank() == mat.rows
+
+
 @pytest.mark.parametrize(
     ("cell_type", "element_type", "order", "variant"),
     [[reference, element, order, variant]

test/utils.py

@@ -27,7 +27,7 @@
         "matrix discontinuous Lagrange": {"equispaced": range(3), "lobatto": range(3)},
         "symmetric matrix discontinuous Lagrange": {
             "equispaced": range(3), "lobatto": range(3)},
-        "bubble": {"equispaced": range(2, 5), "lobatto": range(2, 3)},
+        "bubble": {"equispaced": range(3, 5), "lobatto": range(3, 3)},
         "bubble enriched Lagrange": {"equispaced": range(1, 3), "lobatto": range(1, 3)},
         "bubble enriched vector Lagrange":
             {"equispaced": range(1, 3), "lobatto": range(1, 3)},
@@ -60,7 +60,7 @@
         "vP": {"equispaced": range(3), "lobatto": range(3)},
         "matrix discontinuous Lagrange": {"equispaced": range(2), "lobatto": range(3)},
         "symmetric matrix discontinuous Lagrange": {"equispaced": range(2), "lobatto": range(3)},
-        "bubble": {"equispaced": [3], "lobatto": [3]},
+        "bubble": {"equispaced": [4], "lobatto": [4]},
         "CR": {"equispaced": [1]},
         "Regge": {None: range(3)},
         "Nedelec1": {"equispaced": range(1, 3)},