more robust tests for HCT/HCT-red (#94)

Author: pbrubeck

test/unit/test_hct.py

@@ -1,9 +1,10 @@
 import pytest
 import numpy
 
-from FIAT import HsiehCloughTocher as HCT
-from FIAT.reference_element import ufc_simplex, make_lattice
+from FIAT import RestrictedElement, HsiehCloughTocher as HCT
+from FIAT.reference_element import ufc_simplex
 from FIAT.functional import PointEvaluation
+from FIAT.macro import CkPolynomialSet
 
 
 @pytest.fixture
@@ -13,12 +14,28 @@ def cell():
     return K
 
 
+def span_greater_equal(A, B):
+    # span(A) >= span(B)
+    _, residual, *_ = numpy.linalg.lstsq(A.reshape(A.shape[0], -1).T,
+                                         B.reshape(B.shape[0], -1).T)
+    return numpy.allclose(residual, 0)
+
+
+def make_points(K, degree):
+    top = K.get_topology()
+    pts = []
+    for dim in top:
+        for entity in top[dim]:
+            pts.extend(K.make_points(dim, entity, degree))
+    return pts
+
+
 @pytest.mark.parametrize("reduced", (False, True))
 def test_hct_constant(cell, reduced):
     # Test that bfs associated with point evaluation sum up to 1
     fe = HCT(cell, reduced=reduced)
 
-    pts = make_lattice(cell.get_vertices(), 3)
+    pts = make_points(cell, 4)
     tab = fe.tabulate(2, pts)
 
     coefs = numpy.zeros((fe.space_dimension(),))
@@ -33,3 +50,27 @@ def test_hct_constant(cell, reduced):
         expected = 1 if sum(alpha) == 0 else 0
         vals = numpy.dot(coefs, tab[alpha])
         assert numpy.allclose(vals, expected)
+
+
+@pytest.mark.parametrize("reduced", (False, True))
+def test_full_polynomials(cell, reduced):
+    # Test that HCT/HCT-red contains all cubics/quadratics
+    fe = HCT(cell, reduced=reduced)
+    if reduced:
+        fe = RestrictedElement(fe, restriction_domain="vertex")
+
+    ref_complex = fe.get_reference_complex()
+    pts = make_points(ref_complex, 4)
+    tab = fe.tabulate(0, pts)[(0, 0)]
+
+    degree = fe.degree()
+    if reduced:
+        degree -= 1
+
+    P = CkPolynomialSet(cell, degree, variant="bubble")
+    P_tab = P.tabulate(pts)[(0, 0)]
+    assert span_greater_equal(tab, P_tab)
+
+    C1 = CkPolynomialSet(ref_complex, degree, order=1, variant="bubble")
+    C1_tab = C1.tabulate(pts)[(0, 0)]
+    assert span_greater_equal(tab, C1_tab)

test/unit/test_stokes_complex.py

@@ -35,10 +35,8 @@ def rAQ(cell):
 
 def span_greater_equal(A, B):
     # span(A) >= span(B)
-    dimA = A.shape[0]
-    dimB = B.shape[0]
-    _, residual, *_ = numpy.linalg.lstsq(A.reshape(dimA, -1).T,
-                                         B.reshape(dimB, -1).T)
+    _, residual, *_ = numpy.linalg.lstsq(A.reshape(A.shape[0], -1).T,
+                                         B.reshape(B.shape[0], -1).T)
     return numpy.allclose(residual, 0)