diff --git a/FIAT/gauss_legendre.py b/FIAT/gauss_legendre.py index 6183bf610..7fe845034 100644 --- a/FIAT/gauss_legendre.py +++ b/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) diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py index a602ae08d..b7ccf7c30 100644 --- a/FIAT/gauss_lobatto_legendre.py +++ b/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) diff --git a/FIAT/recursive_points.py b/FIAT/recursive_points.py index 3ef47e24f..1461ba0e9 100644 --- a/FIAT/recursive_points.py +++ b/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,27 +103,20 @@ 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 @@ -130,23 +124,13 @@ def make_recursive_point_set(family): 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 diff --git a/test/unit/test_gauss_legendre.py b/test/unit/test_gauss_legendre.py index b7ae95d09..1037581ef 100644 --- a/test/unit/test_gauss_legendre.py +++ b/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__)) diff --git a/test/unit/test_gauss_lobatto_legendre.py b/test/unit/test_gauss_lobatto_legendre.py index b72a692c5..3cd5b9a2a 100644 --- a/test/unit/test_gauss_lobatto_legendre.py +++ b/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__))