Merged in dham/gauss_legendre (pull request #17)

gauss_legendre quadrature and elements

FIAT/__init__.py

@@ -17,6 +17,7 @@
 from FIAT.hermite import CubicHermite
 from FIAT.lagrange import Lagrange
 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre
+from FIAT.gauss_legendre import GaussLegendre
 from FIAT.morley import Morley
 from FIAT.nedelec import Nedelec
 from FIAT.nedelec_second_kind import NedelecSecondKind
@@ -52,6 +53,7 @@
                       "Hermite": CubicHermite,
                       "Lagrange": Lagrange,
                       "Gauss-Lobatto-Legendre": GaussLobattoLegendre,
+                      "Gauss-Legendre": GaussLegendre,
                       "Morley": Morley,
                       "Nedelec 1st kind H(curl)": Nedelec,
                       "Nedelec 2nd kind H(curl)": NedelecSecondKind,

FIAT/gauss_legendre.py

@@ -0,0 +1,44 @@
+# Copyright (C) 2015 Imperial College London and others.
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+#
+# Written by David A. Ham (david.ham@imperial.ac.uk), 2015
+
+from __future__ import absolute_import, print_function, division
+from . import finite_element, polynomial_set, dual_set, functional, quadrature
+from .reference_element import LINE
+
+
+class GaussLegendreDualSet(dual_set.DualSet):
+    """The dual basis for 1D discontinuous elements with nodes at the
+    Gauss-Legendre points."""
+    def __init__(self, ref_el, degree):
+        entity_ids = {0: {0: [], 1: []},
+                      1: {0: list(range(0, degree+1))}}
+        l = quadrature.GaussLegendreQuadratureLineRule(ref_el, degree+1)
+        nodes = [functional.PointEvaluation(ref_el, x) for x in l.pts]
+
+        dual_set.DualSet.__init__(self, nodes, ref_el, entity_ids)
+
+
+class GaussLegendre(finite_element.FiniteElement):
+    """1D discontinuous element with nodes at the Gauss-Legendre points."""
+    def __init__(self, ref_el, degree):
+        if ref_el.shape != LINE:
+            raise ValueError("Gauss-Legendre elements are only defined in one dimension.")
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
+        dual = GaussLegendreDualSet(ref_el, degree)
+        finite_element.FiniteElement.__init__(self, poly_set, dual, degree)

FIAT/quadrature.py

@@ -111,6 +111,33 @@ def __init__(self, ref_el, m):
         QuadratureRule.__init__(self, ref_el, xs, ws)
 
 
+class GaussLegendreQuadratureLineRule(QuadratureRule):
+    """Produce the Gauss--Legendre quadrature rules on the interval using
+    the implementation in numpy. This facilitates implementing
+    discontinuous spectral elements.
+
+    The quadrature rule uses m points for a degree of precision of 2m-1.
+    """
+    def __init__(self, ref_el, m):
+        if m < 1:
+            raise ValueError(
+                "Gauss-Legendre quadrature invalid for fewer than 2 points")
+
+        xs_ref, ws_ref = numpy.polynomial.legendre.leggauss(m)
+
+        A, b = reference_element.make_affine_mapping(((-1.,), (1.)),
+                                                     ref_el.get_vertices())
+
+        mapping = lambda x: numpy.dot(A, x) + b
+
+        scale = numpy.linalg.det(A)
+
+        xs = tuple([tuple(mapping(x_ref)[0]) for x_ref in xs_ref])
+        ws = tuple([scale * w for w in ws_ref])
+
+        QuadratureRule.__init__(self, ref_el, xs, ws)
+
+
 class CollapsedQuadratureTriangleRule(QuadratureRule):
     """Implements the collapsed quadrature rules defined in
     Karniadakis & Sherwin by mapping products of Gauss-Jacobi rules

test/unit/test_gauss_legendre.py

@@ -0,0 +1,49 @@
+# Copyright (C) 2016 Imperial College London and others
+#
+# This file is part of FIAT.
+#
+# FIAT is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FIAT is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FIAT. If not, see <http://www.gnu.org/licenses/>.
+#
+# Authors:
+#
+# David Ham
+
+from __future__ import absolute_import, print_function, division
+
+import pytest
+import numpy as np
+
+
+@pytest.mark.parametrize("degree", range(1, 7))
+def test_gl_basis_values(degree):
+    """Ensure that integrating a simple monomial produces the expected results."""
+    from FIAT import ufc_simplex, GaussLegendre, make_quadrature
+
+    s = ufc_simplex(1)
+    q = make_quadrature(s, degree + 1)
+
+    fe = GaussLegendre(s, degree)
+    tab = fe.tabulate(0, q.pts)[(0,)]
+
+    for test_degree in range(degree + 1):
+        coefs = [n(lambda x: x[0]**test_degree) for n in fe.dual.nodes]
+        integral = np.dot(coefs, np.dot(tab, q.wts))
+        reference = np.dot([x[0]**test_degree
+                            for x in q.pts], q.wts)
+        assert np.allclose(integral, reference, rtol=1e-14)
+
+
+if __name__ == '__main__':
+    import os
+    pytest.main(os.path.abspath(__file__))

test/unit/test_gauss_lobatto_legendre.py

@@ -26,7 +26,7 @@
 
 
 @pytest.mark.parametrize("degree", range(1, 7))
-def test_basis_values(degree):
+def test_gll_basis_values(degree):
     """Ensure that integrating a simple monomial produces the expected results."""
     from FIAT import ufc_simplex, GaussLobattoLegendre, make_quadrature
 

test/unit/test_quadrature.py

@@ -178,6 +178,18 @@ def test_gauss_lobatto_legendre_quadrature(interval, points, degree):
     assert numpy.round(q.integrate(lambda x: x[0]**degree) - 1./(degree+1), 14) == 0.
 
 
+@pytest.mark.parametrize(("points, degree"), ((p, d)
+                                              for p in range(2, 10)
+                                              for d in range(2*p)))
+def test_gauss_legendre_quadrature(interval, points, degree):
+    """Check that the quadrature rules correctly integrate all the right
+    polynomial degrees."""
+
+    q = FIAT.quadrature.GaussLegendreQuadratureLineRule(interval, points)
+
+    assert numpy.round(q.integrate(lambda x: x[0]**degree) - 1./(degree+1), 14) == 0.
+
+
 if __name__ == '__main__':
     import os
     pytest.main(os.path.abspath(__file__))