Test GLL/GL edge dofs

FIAT/gauss_legendre.py

@@ -21,7 +21,7 @@ class GaussLegendrePointSet(recursive_points.RecursivePointSet):
     def __init__(self):
         ref_el = UFCInterval()
         lr = quadrature.GaussLegendreQuadratureLineRule
-        f = lambda n: lr(ref_el, n + 1).pts
+        f = lambda n: lr(ref_el, n + 1).get_points()
         super(GaussLegendrePointSet, self).__init__(f)
 
 

FIAT/gauss_lobatto_legendre.py

@@ -21,7 +21,7 @@ class GaussLobattoLegendrePointSet(recursive_points.RecursivePointSet):
     def __init__(self):
         ref_el = UFCInterval()
         lr = quadrature.GaussLobattoLegendreQuadratureLineRule
-        f = lambda n: lr(ref_el, n + 1).pts if n else None
+        f = lambda n: lr(ref_el, n + 1).get_points() if n else None
         super(GaussLobattoLegendrePointSet, self).__init__(f)
 
 

FIAT/recursive_points.py

@@ -37,25 +37,26 @@ def multiindex_equal(d, isum, imin=0):
 
 
 class RecursivePointSet(object):
-    """Family of points on the unit interval.  This class essentially is a
+    """Class to construct recursive points on simplices based on a family of
+    points on the unit interval.  This class essentially is a
     lazy-evaluate-and-cache dictionary: the user passes a routine to evaluate
     entries for unknown keys """
 
-    def __init__(self, f):
-        self._f = f
+    def __init__(self, family):
+        self._family = family
         self._cache = {}
 
     def interval_points(self, degree):
         try:
             return self._cache[degree]
         except KeyError:
-            x = self._f(degree)
-            if x is None:
-                x_ro = x
-            else:
-                x_ro = numpy.array(x).flatten()
-                x_ro.setflags(write=False)
-            return self._cache.setdefault(degree, x_ro)
+            x = self._family(degree)
+            if x is not None:
+                if not isinstance(x, numpy.ndarray):
+                    x = numpy.array(x)
+                x = x.reshape((-1,))
+                x.setflags(write=False)
+            return self._cache.setdefault(degree, x)
 
     def _recursive(self, alpha):
         """The barycentric (d-1)-simplex coordinates for a
@@ -102,51 +103,34 @@ def make_points(self, ref_el, dim, entity_id, order):
             raise ValueError("illegal dimension")
 
 
-def make_recursive_point_set(family):
-    from FIAT import quadrature, reference_element
-    ref_el = reference_element.UFCInterval()
-    if family == "equispaced":
-        f = lambda n: numpy.linspace(0.0, 1.0, n + 1)
-    elif family == "dg_equispaced":
-        f = lambda n: numpy.linspace(1.0/(n+2.0), (n+1.0)/(n+2.0), n + 1)
-    elif family == "gl":
-        lr = quadrature.GaussLegendreQuadratureLineRule
-        f = lambda n: lr(ref_el, n + 1).pts
-    elif family == "gll":
-        lr = quadrature.GaussLobattoLegendreQuadratureLineRule
-        f = lambda n: lr(ref_el, n + 1).pts if n else None
-    else:
-        raise ValueError("Invalid node family %s" % family)
-    return RecursivePointSet(f)
-
-
 if __name__ == "__main__":
-    from FIAT import reference_element
+    from FIAT import reference_element, quadrature
     from matplotlib import pyplot as plt
+
+    interval = reference_element.UFCInterval()
+    cg = lambda n: numpy.linspace(0.0, 1.0, n + 1)
+    dg = lambda n: numpy.linspace(1.0/(n+2.0), (n+1.0)/(n+2.0), n + 1)
+    gll = lambda n: quadrature.GaussLegendreQuadratureLineRule(interval, n + 1).get_points()
+    gl = lambda n: quadrature.GaussLobattoLegendreQuadratureLineRule(interval, n + 1).get_points() if n else None
+
+    order = 11
+    cg_points = RecursivePointSet(gll)
+    dg_points = RecursivePointSet(gl)
+
     ref_el = reference_element.ufc_simplex(2)
     h = numpy.sqrt(3)
     s = 2*h/3
     ref_el.vertices = [(0, s), (-1.0, s-h), (1.0, s-h)]
     x = numpy.array(ref_el.vertices + ref_el.vertices[:1])
     plt.plot(x[:, 0], x[:, 1], "k")
 
-    order = 7
-    rule = "gll"
-    dg_rule = "gl"
-
-    # rule = "equispaced"
-    # dg_rule = "dg_equispaced"
-
-    family = make_recursive_point_set(rule)
-    dg_family = make_recursive_point_set(dg_rule)
-
     for d in range(1, 4):
-        print(family.make_points(reference_element.ufc_simplex(d), d, 0, d))
+        assert cg_points.make_points(reference_element.ufc_simplex(d), d, 0, d) == []
 
     topology = ref_el.get_topology()
     for dim in topology:
         for entity in topology[dim]:
-            pts = family.make_points(ref_el, dim, entity, order)
+            pts = cg_points.make_points(ref_el, dim, entity, order)
             if len(pts):
                 x = numpy.array(pts)
                 for r in range(1, 3):
@@ -168,7 +152,7 @@ def make_recursive_point_set(family):
     x0 = sum(x[:d])/d
     plt.plot(x[:, 0], x[:, 1], "k")
 
-    pts = dg_family.recursive_points(ref_el.vertices, order)
+    pts = dg_points.recursive_points(ref_el.vertices, order)
     x = numpy.array(pts)
     for r in range(1, 3):
         th = r * (2*numpy.pi)/3

test/unit/test_gauss_legendre.py

@@ -23,20 +23,24 @@
 import numpy as np
 
 
-@pytest.mark.parametrize("dim, degree", sum(([pytest.param(d, k) for k in range(0, 8-d)] for d in range(1, 4)), []))
-def test_gl_basis_values(dim, degree):
-    """Ensure that integrating a simple monomial produces the expected results."""
-    from FIAT import ufc_simplex, GaussLegendre, make_quadrature
-
+def symmetric_simplex(dim):
+    from FIAT import ufc_simplex
     s = ufc_simplex(dim)
-    r3 = 3.0 ** 0.5
-    r6 = 6.0 ** 0.5
+    r = lambda x: x ** 0.5
     if dim == 2:
-        s.vertices = [(0.0, 0.0), (-1.0, -r3), (1.0, -r3)]
+        s.vertices = [(0.0, 0.0), (-1.0, -r(3.0)), (1.0, -r(3.0))]
     elif dim == 3:
-        s.vertices = [(r3/3, 0.0, 0.0), (-r3/6, 0.5, 0.0),
-                      (-r3/6, -0.5, 0.0), (0.0, 0.0, r6/3)]
+        s.vertices = [(r(3.0)/3, 0.0, 0.0), (-r(3.0)/6, 0.5, 0.0),
+                      (-r(3.0)/6, -0.5, 0.0), (0.0, 0.0, r(6.0)/3)]
+    return s
+
 
+@pytest.mark.parametrize("dim, degree", sum(([(d, p) for p in range(0, 8-d)] for d in range(1, 4)), []))
+def test_gl_basis_values(dim, degree):
+    """Ensure that integrating a simple monomial produces the expected results."""
+    from FIAT import GaussLegendre, make_quadrature
+
+    s = symmetric_simplex(dim)
     q = make_quadrature(s, degree + 1)
     fe = GaussLegendre(s, degree)
     tab = fe.tabulate(0, q.pts)[(0,)*dim]
@@ -49,6 +53,34 @@ def test_gl_basis_values(dim, degree):
         assert np.allclose(integral, reference, rtol=1e-14)
 
 
+@pytest.mark.parametrize("dim, degree", [(1, 4), (2, 4), (3, 4)])
+def test_symmetry(dim, degree):
+    """ Ensure the dual basis has the right symmetry."""
+    from FIAT import GaussLegendre, quadrature, expansions, ufc_simplex
+
+    s = symmetric_simplex(dim)
+    fe = GaussLegendre(s, degree)
+    ndof = fe.space_dimension()
+    assert ndof == expansions.polynomial_dimension(s, degree)
+
+    points = np.zeros((ndof, dim), "d")
+    for i, node in enumerate(fe.dual_basis()):
+        points[i, :], = node.get_point_dict().keys()
+
+    # Test that edge DOFs are located at the GL quadrature points
+    lr = quadrature.GaussLegendreQuadratureLineRule(ufc_simplex(1), degree + 1)
+    quadrature_points = lr.pts
+
+    entity_dofs = fe.entity_dofs()
+    edge_dofs = entity_dofs[1]
+    for entity in edge_dofs:
+        if len(edge_dofs[entity]) > 0:
+            transform = s.get_entity_transform(1, entity)
+            assert np.allclose(points[edge_dofs[entity]], np.array(list(map(transform, quadrature_points))))
+
+    # TODO add rotational symmetry tests
+
+
 if __name__ == '__main__':
     import os
     pytest.main(os.path.abspath(__file__))

test/unit/test_gauss_lobatto_legendre.py

@@ -23,20 +23,24 @@
 import numpy as np
 
 
-@pytest.mark.parametrize("dim, degree", sum(([pytest.param(d, k) for k in range(1, 8-d)] for d in range(1, 4)), []))
-def test_gll_basis_values(dim, degree):
-    """Ensure that integrating a simple monomial produces the expected results."""
-    from FIAT import ufc_simplex, GaussLobattoLegendre, make_quadrature
-
+def symmetric_simplex(dim):
+    from FIAT import ufc_simplex
     s = ufc_simplex(dim)
-    r3 = 3.0 ** 0.5
-    r6 = 6.0 ** 0.5
+    r = lambda x: x ** 0.5
     if dim == 2:
-        s.vertices = [(0.0, 0.0), (-1.0, -r3), (1.0, -r3)]
+        s.vertices = [(0.0, 0.0), (-1.0, -r(3.0)), (1.0, -r(3.0))]
     elif dim == 3:
-        s.vertices = [(r3/3, 0.0, 0.0), (-r3/6, 0.5, 0.0),
-                      (-r3/6, -0.5, 0.0), (0.0, 0.0, r6/3)]
+        s.vertices = [(r(3.0)/3, 0.0, 0.0), (-r(3.0)/6, 0.5, 0.0),
+                      (-r(3.0)/6, -0.5, 0.0), (0.0, 0.0, r(6.0)/3)]
+    return s
+
 
+@pytest.mark.parametrize("dim, degree", sum(([(d, p) for p in range(1, 8-d)] for d in range(1, 4)), []))
+def test_gll_basis_values(dim, degree):
+    """Ensure that integrating a simple monomial produces the expected results."""
+    from FIAT import GaussLobattoLegendre, make_quadrature
+
+    s = symmetric_simplex(dim)
     q = make_quadrature(s, degree + 1)
     fe = GaussLobattoLegendre(s, degree)
     tab = fe.tabulate(0, q.pts)[(0,)*dim]
@@ -49,6 +53,35 @@ def test_gll_basis_values(dim, degree):
         assert np.allclose(integral, reference, rtol=1e-14)
 
 
+@pytest.mark.parametrize("dim, degree", [(1, 4), (2, 4), (3, 4)])
+def test_symmetry(dim, degree):
+    """ Ensure the dual basis has the right symmetry."""
+    from FIAT import GaussLobattoLegendre, quadrature, expansions, ufc_simplex
+
+    s = symmetric_simplex(dim)
+    fe = GaussLobattoLegendre(s, degree)
+    ndof = fe.space_dimension()
+    assert ndof == expansions.polynomial_dimension(s, degree)
+
+    points = np.zeros((ndof, dim), "d")
+    for i, node in enumerate(fe.dual_basis()):
+        points[i, :], = node.get_point_dict().keys()
+
+    # Test that edge DOFs are located at the GLL quadrature points
+    lr = quadrature.GaussLobattoLegendreQuadratureLineRule(ufc_simplex(1), degree + 1)
+    # Edge DOFs are ordered with the two vertex DOFs followed by the interior DOFs
+    quadrature_points = lr.pts[::degree] + lr.pts[1:-1]
+
+    entity_dofs = fe.entity_closure_dofs()
+    edge_dofs = entity_dofs[1]
+    for entity in edge_dofs:
+        if len(edge_dofs[entity]) > 0:
+            transform = s.get_entity_transform(1, entity)
+            assert np.allclose(points[edge_dofs[entity]], np.array(list(map(transform, quadrature_points))))
+
+    # TODO add rotational symmetry tests on each facet
+
+
 if __name__ == '__main__':
     import os
     pytest.main(os.path.abspath(__file__))