From ffeef6cd0f2d12df6da2c3e9255758d975f3e0f3 Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Thu, 27 Jun 2019 18:47:14 +0100 Subject: [PATCH 1/2] proper testing and EGL/Edge elements --- FIAT/__init__.py | 5 + FIAT/edge.py | 176 ++++++++++++++++++++++ FIAT/extended_gauss_legendre.py | 44 ++++++ FIAT/quadrature.py | 35 ++++- test/unit/test_edge.py | 72 +++++++++ test/unit/test_extended_gauss_legendre.py | 60 ++++++++ 6 files changed, 391 insertions(+), 1 deletion(-) create mode 100644 FIAT/edge.py create mode 100644 FIAT/extended_gauss_legendre.py create mode 100644 test/unit/test_edge.py create mode 100644 test/unit/test_extended_gauss_legendre.py diff --git a/FIAT/__init__.py b/FIAT/__init__.py index f5c1677ce..27921df68 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -38,6 +38,8 @@ from FIAT.mixed import MixedElement # noqa: F401 from FIAT.restricted import RestrictedElement # noqa: F401 from FIAT.quadrature_element import QuadratureElement # noqa: F401 +from FIAT.extended_gauss_legendre import ExtendedGaussLegendre +from FIAT.edge import EdgeGaussLobattoLegendre, EdgeExtendedGaussLegendre # Important functionality from FIAT.quadrature import make_quadrature # noqa: F401 @@ -65,6 +67,9 @@ "Lagrange": Lagrange, "Gauss-Lobatto-Legendre": GaussLobattoLegendre, "Gauss-Legendre": GaussLegendre, + "Extended-Gauss-Legendre": ExtendedGaussLegendre, + "Gauss-Lobatto-Legendre Edge": EdgeGaussLobattoLegendre, + "Extended-Gauss-Legendre Edge": EdgeExtendedGaussLegendre, "Morley": Morley, "Nedelec 1st kind H(curl)": Nedelec, "Nedelec 2nd kind H(curl)": NedelecSecondKind, diff --git a/FIAT/edge.py b/FIAT/edge.py new file mode 100644 index 000000000..4c518da75 --- /dev/null +++ b/FIAT/edge.py @@ -0,0 +1,176 @@ +# Copyright (C) 2019 Chris Eldred (Inria Grenoble Rhone-Alpes) and Werner Bauer (Inria Rennes) +# +# 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 . +# + +import sympy +import numpy as np +from FIAT.finite_element import FiniteElement +from FIAT.dual_set import make_entity_closure_ids +from FIAT.polynomial_set import mis +from FIAT import quadrature + +x = sympy.symbols('x') + + +def _lagrange_poly(i, xi): + '''Returns the Langrange polynomial P_i(x) of pts xi + which interpolates the values (0,0,..1,..,0) where 1 is at the ith support point + :param i: non-zero location + :param xi: set of N support points + ''' + index = list(range(len(xi))) + index.pop(i) + return sympy.prod([(x-xi[j][0])/(xi[i][0]-xi[j][0]) for j in index]) + + +def _lagrange_basis(spts): + symbas = [] + for i in range(len(spts)): + symbas.append(_lagrange_poly(i, spts)) + return symbas + + +def _create_compatible_l2_basis(cg_symbas): + ncg_basis = len(cg_symbas) + symbas = [] + for i in range(1, ncg_basis): + basis = 0 + for j in range(i): + basis = basis + sympy.diff(-cg_symbas[j]) + symbas.append(basis) + return symbas + + +class _EdgeElement(FiniteElement): + + def __new__(cls, ref_el, degree): + dim = ref_el.get_spatial_dimension() + if dim == 1: + self = super().__new__(cls) + return self + else: + raise IndexError("only intervals supported for _IntervalElement") + + def tabulate(self, order, points, entity=None): + + if entity is None: + entity = (self.ref_el.get_spatial_dimension(), 0) + entity_dim, entity_id = entity + dim = self.ref_el.get_spatial_dimension() + + transform = self.ref_el.get_entity_transform(entity_dim, entity_id) + points = list(map(transform, points)) + + phivals = {} + for o in range(order + 1): + alphas = mis(dim, o) + for alpha in alphas: + try: + basis = self.basis[alpha] + except KeyError: + basis = sympy.diff(self.basis[(0,)], x) + self.basis[alpha] = basis + T = np.zeros((len(basis), len(points))) + for i in range(len(points)): + subs = {x: points[i][0]} + for j, f in enumerate(basis): + T[j, i] = f.evalf(subs=subs) + phivals[alpha] = T + + return phivals + + def degree(self): + return self._degree + + def get_nodal_basis(self): + raise NotImplementedError("get_nodal_basis not implemented for _IntervalElement") + + def get_dual_set(self): + raise NotImplementedError("get_dual_set is not implemented for _IntervalElement") + + def get_coeffs(self): + raise NotImplementedError("get_coeffs not implemented for _IntervalElement") + + def entity_dofs(self): + """Return the map of topological entities to degrees of + freedom for the finite element.""" + return self.entity_ids + + def entity_closure_dofs(self): + """Return the map of topological entities to degrees of + freedom on the closure of those entities for the finite element.""" + return self.entity_closure_ids + + def value_shape(self): + return () + + def dmats(self): + raise NotImplementedError + + def get_num_members(self, arg): + raise NotImplementedError + + def space_dimension(self): + return len(self.basis[(0,)]) + + def __init__(self, ref_el, degree): + + dim = ref_el.get_spatial_dimension() + topology = ref_el.get_topology() + + if not dim == 1: + raise IndexError("only intervals supported for DMSE") + + formdegree = 1 + + # This is general code to build empty entity_ids + entity_ids = {} + for dim in range(dim+1): + entity_ids[dim] = {} + for entity in sorted(topology[dim]): + entity_ids[dim][entity] = [] + + # The only dofs for DMSE are with the interval! + entity_ids[dim][0] = list(range(degree + 1)) + + # Build the basis + # This is a dictionary basis[o] that gives the "basis" functions for derivative order o, where o=0 is just the basis itself + # This is filled as needed for o > 0 by tabulate + + self.entity_ids = entity_ids + self.entity_closure_ids = make_entity_closure_ids(ref_el, entity_ids) + self._degree = degree + + super(_EdgeElement, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree) + + +class EdgeGaussLobattoLegendre(_EdgeElement): + def __init__(self, ref_el, degree): + super(EdgeGaussLobattoLegendre, self).__init__(ref_el, degree) + cont_pts = quadrature.GaussLobattoLegendreQuadratureLineRule(ref_el, degree+2).pts + cont_basis = _lagrange_basis(cont_pts) + basis = _create_compatible_l2_basis(cont_basis) + self.basis = {(0,): sympy.Array(basis)} + + +class EdgeExtendedGaussLegendre(_EdgeElement): + def __init__(self, ref_el, degree): + super(EdgeExtendedGaussLegendre, self).__init__(ref_el, degree) + cont_pts = quadrature.ExtendedGaussLegendreQuadratureLineRule(ref_el, degree+2).pts + cont_basis = _lagrange_basis(cont_pts) + basis = _create_compatible_l2_basis(cont_basis) + self.basis = {(0,): sympy.Array(basis)} diff --git a/FIAT/extended_gauss_legendre.py b/FIAT/extended_gauss_legendre.py new file mode 100644 index 000000000..aba06a7ed --- /dev/null +++ b/FIAT/extended_gauss_legendre.py @@ -0,0 +1,44 @@ +# Copyright (C) 2019 Chris Eldred (Inria Grenoble Rhone-Alpes) and Werner Bauer (Inria Rennes) +# +# 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 . +# +# Written by Chris Eldred (Inria Grenoble Rhone-Alpes) and Werner Bauer (Inria Rennes) 2019 + +from FIAT import finite_element, polynomial_set, dual_set, functional, quadrature +from FIAT.reference_element import LINE + + +class ExtendedGaussLegendreDualSet(dual_set.DualSet): + """The dual basis for 1D continuous elements with nodes at the + Extended-Gauss-Legendre points.""" + def __init__(self, ref_el, degree): + entity_ids = {0: {0: [0], 1: [degree]}, + 1: {0: list(range(1, degree))}} + lr = quadrature.ExtendedGaussLegendreQuadratureLineRule(ref_el, degree+1) + nodes = [functional.PointEvaluation(ref_el, x) for x in lr.pts] + + super(ExtendedGaussLegendreDualSet, self).__init__(nodes, ref_el, entity_ids) + + +class ExtendedGaussLegendre(finite_element.CiarletElement): + """1D continuous element with nodes at the Extended-Gauss-Legendre points.""" + def __init__(self, ref_el, degree): + if ref_el.shape != LINE: + raise ValueError("Extended-Gauss-Legendre elements are only defined in one dimension.") + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) + dual = ExtendedGaussLegendreDualSet(ref_el, degree) + formdegree = 0 # 0-form + super(ExtendedGaussLegendre, self).__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index 95292bb0b..e82b98189 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -108,7 +108,7 @@ def __init__(self, ref_el, m): class GaussLegendreQuadratureLineRule(QuadratureRule): - """Produce the Gauss--Legendre quadrature rules on the interval using + """Produce the Gauss-Legendre quadrature rules on the interval using the implementation in numpy. This facilitates implementing discontinuous spectral elements. @@ -134,6 +134,39 @@ def __init__(self, ref_el, m): QuadratureRule.__init__(self, ref_el, xs, ws) +class ExtendedGaussLegendreQuadratureLineRule(QuadratureRule): + """Produce the Extended-Gauss-Legendre quadrature rules on the interval using + the implementation in numpy. This facilitates implementing + dual continuous mimetic spectral elements. + + The quadrature rule uses m points for a degree of precision of 2(m-2)-3 = 2m - 7. + Thus, not recommend for actual use (but is useful to define the Extended-Gauss-Legendre elements) + """ + def __init__(self, ref_el, m): + if m < 3: + raise ValueError( + "Extended-Gauss-Legendre quadrature invalid for fewer than 3 points") + + xs_ref, ws_ref = numpy.polynomial.legendre.leggauss(m-2) + 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) + + ws = list([scale * w for w in ws_ref]) + xs = [tuple(mapping(x_ref)[0]) for x_ref in xs_ref] + + xs.insert(0, (0.,)) + xs.append((1.,)) + # The integral of these basis functions over the interval is 0! + ws.insert(0, 0.) + ws.append(0.) + + QuadratureRule.__init__(self, ref_el, tuple(xs), tuple(ws)) + + class CollapsedQuadratureTriangleRule(QuadratureRule): """Implements the collapsed quadrature rules defined in Karniadakis & Sherwin by mapping products of Gauss-Jacobi rules diff --git a/test/unit/test_edge.py b/test/unit/test_edge.py new file mode 100644 index 000000000..bc7230c84 --- /dev/null +++ b/test/unit/test_edge.py @@ -0,0 +1,72 @@ +# 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 . +# +# Authors: +# +# Christopher Eldred + +import pytest +import math +import sympy +import FIAT + +x = sympy.symbols('x') + + +@pytest.mark.parametrize("degree", range(0, 7)) +def test_gll_edge_basis_values(degree): + """Ensure that edge elements are histopolatory""" + + s = FIAT.ufc_simplex(1) + + fe = FIAT.EdgeGaussLobattoLegendre(s, degree) + + basis = fe.basis[(0,)] + cont_pts = FIAT.quadrature.GaussLobattoLegendreQuadratureLineRule(s, degree + 2).pts + + for i in range(len(cont_pts)-1): + for j in range(basis.shape[0]): + int_sub = sympy.integrate(basis[j], (x, cont_pts[i][0], cont_pts[i+1][0])) + if i == j: + assert(math.isclose(int_sub, 1.)) + else: + assert(math.isclose(int_sub, 0., abs_tol=1e-9)) + + +@pytest.mark.parametrize("degree", range(2, 7)) +def test_egl_edge_basis_values(degree): + """Ensure that edge elements are histopolatory""" + + s = FIAT.ufc_simplex(1) + + fe = FIAT.EdgeExtendedGaussLegendre(s, degree) + + basis = fe.basis[(0,)] + cont_pts = FIAT.quadrature.ExtendedGaussLegendreQuadratureLineRule(s, degree + 2).pts + + for i in range(len(cont_pts)-1): + for j in range(basis.shape[0]): + int_sub = sympy.integrate(basis[j], (x, cont_pts[i][0], cont_pts[i+1][0])) + if i == j: + assert(math.isclose(int_sub, 1.)) + else: + assert(math.isclose(int_sub, 0., abs_tol=1e-9)) + + +if __name__ == '__main__': + import os + pytest.main(os.path.abspath(__file__)) diff --git a/test/unit/test_extended_gauss_legendre.py b/test/unit/test_extended_gauss_legendre.py new file mode 100644 index 000000000..25d4ad949 --- /dev/null +++ b/test/unit/test_extended_gauss_legendre.py @@ -0,0 +1,60 @@ +# Copyright (C) 2019 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 . +# +# Authors: +# +# Christopher Eldred + +import pytest +import numpy as np +import FIAT +from FIAT.reference_element import UFCInterval + + +@pytest.mark.parametrize("degree", range(2, 7)) +def test_egl_basis_values(degree): + """Ensure that integrating a simple monomial produces the expected results.""" + from FIAT import ufc_simplex, ExtendedGaussLegendre, make_quadrature + + s = ufc_simplex(1) + q = make_quadrature(s, degree + 1) + + fe = ExtendedGaussLegendre(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) + + +@pytest.mark.parametrize(("points, degree"), ((p, d) + for p in range(3, 10) + for d in range(2*p - 7))) +def test_extended_gauss_legendre_quadrature(points, degree): + """Check that the quadrature rules correctly integrate all the right + polynomial degrees.""" + + q = FIAT.quadrature.ExtendedGaussLegendreQuadratureLineRule(UFCInterval(), points) + assert np.round(q.integrate(lambda x: x[0]**degree) - 1./(degree+1), 14) == 0. + + +if __name__ == '__main__': + import os + pytest.main(os.path.abspath(__file__)) From d73a5c3d37c3639b24bc5ffd02da3869806426bd Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Tue, 10 Sep 2019 16:20:45 +0200 Subject: [PATCH 2/2] tests: Attempt to pacify pytest --- test/unit/test_edge.py | 4 ++-- test/unit/test_extended_gauss_legendre.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/test/unit/test_edge.py b/test/unit/test_edge.py index bc7230c84..ee859e07b 100644 --- a/test/unit/test_edge.py +++ b/test/unit/test_edge.py @@ -27,7 +27,7 @@ x = sympy.symbols('x') -@pytest.mark.parametrize("degree", range(0, 7)) +@pytest.mark.parametrize("degree", list(range(0, 7))) def test_gll_edge_basis_values(degree): """Ensure that edge elements are histopolatory""" @@ -47,7 +47,7 @@ def test_gll_edge_basis_values(degree): assert(math.isclose(int_sub, 0., abs_tol=1e-9)) -@pytest.mark.parametrize("degree", range(2, 7)) +@pytest.mark.parametrize("degree", list(range(2, 7))) def test_egl_edge_basis_values(degree): """Ensure that edge elements are histopolatory""" diff --git a/test/unit/test_extended_gauss_legendre.py b/test/unit/test_extended_gauss_legendre.py index 25d4ad949..055494ff5 100644 --- a/test/unit/test_extended_gauss_legendre.py +++ b/test/unit/test_extended_gauss_legendre.py @@ -25,7 +25,7 @@ from FIAT.reference_element import UFCInterval -@pytest.mark.parametrize("degree", range(2, 7)) +@pytest.mark.parametrize("degree", list(range(2, 7))) def test_egl_basis_values(degree): """Ensure that integrating a simple monomial produces the expected results.""" from FIAT import ufc_simplex, ExtendedGaussLegendre, make_quadrature @@ -44,9 +44,9 @@ def test_egl_basis_values(degree): assert np.allclose(integral, reference, rtol=1e-14) -@pytest.mark.parametrize(("points, degree"), ((p, d) - for p in range(3, 10) - for d in range(2*p - 7))) +@pytest.mark.parametrize(("points, degree"), list((p, d) + for p in range(3, 10) + for d in range(2*p - 7))) def test_extended_gauss_legendre_quadrature(points, degree): """Check that the quadrature rules correctly integrate all the right polynomial degrees."""