diff --git a/.github/workflows/pythonapp.yml b/.github/workflows/pythonapp.yml index 53c7aac61..39dd4bcc9 100644 --- a/.github/workflows/pythonapp.yml +++ b/.github/workflows/pythonapp.yml @@ -12,7 +12,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, macos-latest] - python-version: [3.7, 3.8] + python-version: [3.8, 3.9] steps: - uses: actions/checkout@v2 diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py index 64c6913ef..2bc33de42 100644 --- a/FIAT/barycentric_interpolation.py +++ b/FIAT/barycentric_interpolation.py @@ -95,4 +95,4 @@ def get_expansion_set(ref_el, pts): if ref_el.get_shape() == reference_element.LINE: return LagrangeLineExpansionSet(ref_el, pts) else: - raise Exception("Unknown reference element type.") + raise ValueError("Invalid reference element type.") diff --git a/FIAT/gauss_legendre.py b/FIAT/gauss_legendre.py index dc17c4183..4d3e1c0e8 100644 --- a/FIAT/gauss_legendre.py +++ b/FIAT/gauss_legendre.py @@ -8,38 +8,52 @@ # # Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021 -from FIAT import finite_element, dual_set, functional, quadrature -from FIAT.reference_element import LINE +from FIAT import finite_element, polynomial_set, dual_set, functional +from FIAT.reference_element import POINT, LINE, TRIANGLE, TETRAHEDRON from FIAT.orientation_utils import make_entity_permutations_simplex from FIAT.barycentric_interpolation import LagrangePolynomialSet +from FIAT.reference_element import make_lattice class GaussLegendreDualSet(dual_set.DualSet): - """The dual basis for 1D discontinuous elements with nodes at the - Gauss-Legendre points.""" + """The dual basis for discontinuous elements with nodes at the + (recursive) Gauss-Legendre points.""" + def __init__(self, ref_el, degree): - entity_ids = {0: {0: [], 1: []}, - 1: {0: list(range(0, degree+1))}} - lr = quadrature.GaussLegendreQuadratureLineRule(ref_el, degree+1) - nodes = [functional.PointEvaluation(ref_el, x) for x in lr.pts] + entity_ids = {} entity_permutations = {} - entity_permutations[0] = {0: {0: []}, 1: {0: []}} - entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree + 1)} + top = ref_el.get_topology() + for dim in sorted(top): + entity_ids[dim] = {} + entity_permutations[dim] = {} + perms = make_entity_permutations_simplex(dim, degree + 1 if dim == len(top)-1 else -1) + for entity in sorted(top[dim]): + entity_ids[dim][entity] = [] + entity_permutations[dim][entity] = perms + # make nodes by getting points + pts = make_lattice(ref_el.get_vertices(), degree, variant="gl") + nodes = [functional.PointEvaluation(ref_el, x) for x in pts] + entity_ids[dim][0] = list(range(len(nodes))) super(GaussLegendreDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) class GaussLegendre(finite_element.CiarletElement): - """1D discontinuous element with nodes at the Gauss-Legendre points.""" + """Simplicial discontinuous element with nodes at the (recursive) 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.") + if ref_el.shape not in {POINT, LINE, TRIANGLE, TETRAHEDRON}: + raise ValueError("Gauss-Legendre elements are only defined on simplices.") dual = GaussLegendreDualSet(ref_el, degree) - points = [] - for node in dual.nodes: - # Assert singleton point for each node. - pt, = node.get_point_dict().keys() - points.append(pt) - poly_set = LagrangePolynomialSet(ref_el, points) + if ref_el.shape == LINE: + # In 1D we can use the primal basis as the expansion set, + # avoiding any round-off coming from a basis transformation + points = [] + for node in dual.nodes: + # Assert singleton point for each node. + pt, = node.get_point_dict().keys() + points.append(pt) + poly_set = LagrangePolynomialSet(ref_el, points) + else: + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form super(GaussLegendre, self).__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py index b40557a1b..d47f758b0 100644 --- a/FIAT/gauss_lobatto_legendre.py +++ b/FIAT/gauss_lobatto_legendre.py @@ -8,38 +8,27 @@ # # Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021 -from FIAT import finite_element, dual_set, functional, quadrature -from FIAT.reference_element import LINE -from FIAT.orientation_utils import make_entity_permutations_simplex +from FIAT import finite_element, polynomial_set, lagrange +from FIAT.reference_element import LINE, TRIANGLE, TETRAHEDRON from FIAT.barycentric_interpolation import LagrangePolynomialSet -class GaussLobattoLegendreDualSet(dual_set.DualSet): - """The dual basis for 1D continuous elements with nodes at the - Gauss-Lobatto points.""" - def __init__(self, ref_el, degree): - entity_ids = {0: {0: [0], 1: [degree]}, - 1: {0: list(range(1, degree))}} - lr = quadrature.GaussLobattoLegendreQuadratureLineRule(ref_el, degree+1) - nodes = [functional.PointEvaluation(ref_el, x) for x in lr.pts] - entity_permutations = {} - entity_permutations[0] = {0: {0: [0]}, 1: {0: [0]}} - entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree - 1)} - - super(GaussLobattoLegendreDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) - - class GaussLobattoLegendre(finite_element.CiarletElement): - """1D continuous element with nodes at the Gauss-Lobatto points.""" + """Simplicial continuous element with nodes at the (recursive) Gauss-Lobatto points.""" def __init__(self, ref_el, degree): - if ref_el.shape != LINE: - raise ValueError("Gauss-Lobatto-Legendre elements are only defined in one dimension.") - dual = GaussLobattoLegendreDualSet(ref_el, degree) - points = [] - for node in dual.nodes: - # Assert singleton point for each node. - pt, = node.get_point_dict().keys() - points.append(pt) - poly_set = LagrangePolynomialSet(ref_el, points) + if ref_el.shape not in {LINE, TRIANGLE, TETRAHEDRON}: + raise ValueError("Gauss-Lobatto-Legendre elements are only defined on simplices.") + dual = lagrange.LagrangeDualSet(ref_el, degree, variant="gll") + if ref_el.shape == LINE: + # In 1D we can use the primal basis as the expansion set, + # avoiding any round-off coming from a basis transformation + points = [] + for node in dual.nodes: + # Assert singleton point for each node. + pt, = node.get_point_dict().keys() + points.append(pt) + poly_set = LagrangePolynomialSet(ref_el, points) + else: + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) formdegree = 0 # 0-form super(GaussLobattoLegendre, self).__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index 7852bd2bb..12816ad0a 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -14,7 +14,7 @@ class LagrangeDualSet(dual_set.DualSet): simplices of any dimension. Nodes are point evaluation at equispaced points.""" - def __init__(self, ref_el, degree): + def __init__(self, ref_el, degree, variant=None): entity_ids = {} nodes = [] entity_permutations = {} @@ -29,7 +29,7 @@ def __init__(self, ref_el, degree): entity_permutations[dim] = {} perms = {0: [0]} if dim == 0 else make_entity_permutations_simplex(dim, degree - dim) for entity in sorted(top[dim]): - pts_cur = ref_el.make_points(dim, entity, degree) + pts_cur = ref_el.make_points(dim, entity, degree, variant=variant) nodes_cur = [functional.PointEvaluation(ref_el, x) for x in pts_cur] nnodes_cur = len(nodes_cur) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 9bf0b71dd..679c0a462 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -18,6 +18,7 @@ import numpy from FIAT import expansions from FIAT.functional import index_iterator +from FIAT.reference_element import make_lattice def mis(m, n): @@ -155,23 +156,20 @@ def __init__(self, ref_el, degree, shape=tuple()): cur_idx = tuple([cur_bf] + list(idx) + [exp_bf]) coeffs[cur_idx] = 1.0 cur_bf += 1 - # construct dmats if degree == 0: dmats = [numpy.array([[0.0]], "d") for i in range(sd)] else: - pts = ref_el.make_points(sd, 0, degree + sd + 1) + pts = make_lattice(ref_el.get_vertices(), degree, variant="gl") v = numpy.transpose(expansion_set.tabulate(degree, pts)) - vinv = numpy.linalg.inv(v) dv = expansion_set.tabulate_derivatives(degree, pts) dtildes = [[[a[1][i] for a in dvrow] for dvrow in dv] for i in range(sd)] - dmats = [numpy.dot(vinv, numpy.transpose(dtilde)) + dmats = [numpy.linalg.solve(v, numpy.transpose(dtilde)) for dtilde in dtildes] - PolynomialSet.__init__(self, ref_el, degree, embedded_degree, expansion_set, coeffs, dmats) diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index a93b73c1a..5e0787be9 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -8,10 +8,10 @@ # Modified by David A. Ham (david.ham@imperial.ac.uk), 2015 import itertools -import math import numpy +from recursivenodes.quadrature import gaussjacobi -from FIAT import reference_element, expansions, jacobi, orthopoly +from FIAT import reference_element, expansions, orthopoly class QuadratureRule(object): @@ -33,7 +33,7 @@ def get_weights(self): return numpy.array(self.wts) def integrate(self, f): - return sum([w * f(x) for (x, w) in zip(self.pts, self.wts)]) + return sum(w * f(x) for x, w in zip(self.pts, self.wts)) class GaussJacobiQuadratureLineRule(QuadratureRule): @@ -42,8 +42,8 @@ class GaussJacobiQuadratureLineRule(QuadratureRule): def __init__(self, ref_el, m): # this gives roots on the default (-1,1) reference element - # (xs_ref, ws_ref) = compute_gauss_jacobi_rule(a, b, m) - (xs_ref, ws_ref) = compute_gauss_jacobi_rule(0., 0., m) + # (xs_ref, ws_ref) = gaussjacobi(m, a, b) + (xs_ref, ws_ref) = gaussjacobi(m, 0., 0.) Ref1 = reference_element.DefaultLine() A, b = reference_element.make_affine_mapping(Ref1.get_vertices(), @@ -186,8 +186,8 @@ class CollapsedQuadratureTriangleRule(QuadratureRule): from the square to the triangle.""" def __init__(self, ref_el, m): - ptx, wx = compute_gauss_jacobi_rule(0., 0., m) - pty, wy = compute_gauss_jacobi_rule(1., 0., m) + ptx, wx = gaussjacobi(m, 0., 0.) + pty, wy = gaussjacobi(m, 1., 0.) # map ptx , pty pts_ref = [expansions.xi_triangle((x, y)) @@ -213,9 +213,9 @@ class CollapsedQuadratureTetrahedronRule(QuadratureRule): from the cube to the tetrahedron.""" def __init__(self, ref_el, m): - ptx, wx = compute_gauss_jacobi_rule(0., 0., m) - pty, wy = compute_gauss_jacobi_rule(1., 0., m) - ptz, wz = compute_gauss_jacobi_rule(2., 0., m) + ptx, wx = gaussjacobi(m, 0., 0.) + pty, wy = gaussjacobi(m, 1., 0.) + ptz, wz = gaussjacobi(m, 2., 0.) # map ptx , pty pts_ref = [expansions.xi_tetrahedron((x, y, z)) @@ -319,53 +319,3 @@ def make_tensor_product_quadrature(*quad_rules): wts = [numpy.prod(wt_tuple) for wt_tuple in itertools.product(*[q.wts for q in quad_rules])] return QuadratureRule(ref_el, pts, wts) - - -# rule to get Gauss-Jacobi points -def compute_gauss_jacobi_points(a, b, m): - """Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method. - The initial guesses are the Chebyshev points. Algorithm - implemented in Python from the pseudocode given by Karniadakis and - Sherwin""" - x = [] - eps = 1.e-8 - max_iter = 100 - for k in range(0, m): - r = -math.cos((2.0 * k + 1.0) * math.pi / (2.0 * m)) - if k > 0: - r = 0.5 * (r + x[k - 1]) - j = 0 - delta = 2 * eps - while j < max_iter: - s = 0 - for i in range(0, k): - s = s + 1.0 / (r - x[i]) - f = jacobi.eval_jacobi(a, b, m, r) - fp = jacobi.eval_jacobi_deriv(a, b, m, r) - delta = f / (fp - f * s) - - r = r - delta - - if math.fabs(delta) < eps: - break - else: - j = j + 1 - - x.append(r) - return x - - -def compute_gauss_jacobi_rule(a, b, m): - xs = compute_gauss_jacobi_points(a, b, m) - - a1 = math.pow(2, a + b + 1) - a2 = math.gamma(a + m + 1) - a3 = math.gamma(b + m + 1) - a4 = math.gamma(a + b + m + 1) - a5 = math.factorial(m) - a6 = a1 * a2 * a3 / a4 / a5 - - ws = [a6 / (1.0 - x**2.0) / jacobi.eval_jacobi_deriv(a, b, m, x)**2.0 - for x in xs] - - return xs, ws diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 7a4f63e6e..66225a674 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -24,6 +24,7 @@ from collections import defaultdict import operator from math import factorial +from recursivenodes.nodes import _recursive, _decode_family from FIAT.orientation_utils import make_cell_orientation_reflection_map_simplex, make_cell_orientation_reflection_map_tensorproduct @@ -39,22 +40,36 @@ TENSORPRODUCT = 99 +def multiindex_equal(d, isum, imin=0): + """A generator for d-tuple multi-indices whose sum is isum and minimum is imin. + """ + if d <= 0: + return + imax = isum - (d - 1) * imin + if imax < imin: + return + for i in range(imin, imax): + for a in multiindex_equal(d - 1, isum - i, imin=imin): + yield a + (i,) + yield (imin,) * (d - 1) + (imax,) + + def lattice_iter(start, finish, depth): """Generator iterating over the depth-dimensional lattice of integers between start and (finish-1). This works on simplices in - 1d, 2d, 3d, and beyond""" + 0d, 1d, 2d, 3d, and beyond""" if depth == 0: - return + yield tuple() elif depth == 1: for ii in range(start, finish): - yield [ii] + yield (ii,) else: for ii in range(start, finish): for jj in lattice_iter(start, finish - ii, depth - 1): - yield jj + [ii] + yield jj + (ii,) -def make_lattice(verts, n, interior=0): +def make_lattice(verts, n, interior=0, variant=None): """Constructs a lattice of points on the simplex defined by verts. For example, the 1:st order lattice will be just the vertices. The optional argument interior specifies how many points from @@ -62,15 +77,15 @@ def make_lattice(verts, n, interior=0): and interior = 0, this function will return the vertices and midpoint, but with interior = 1, it will only return the midpoint.""" - - vs = numpy.array(verts) - hs = (vs - vs[0])[1:, :] / n - - m = hs.shape[0] - result = [tuple(vs[0] + numpy.array(indices).dot(hs)) - for indices in lattice_iter(interior, n + 1 - interior, m)] - - return result + if variant is None or variant == "equispaced": + variant = "equi" + elif variant == "gll": + variant = "lgl" + family = _decode_family(variant) + D = len(verts) + X = numpy.array(verts) + get_point = lambda alpha: tuple(numpy.dot(_recursive(D - 1, n, alpha, family), X)) + return list(map(get_point, multiindex_equal(D, n, interior))) def linalg_subspace_intersection(A, B): @@ -393,7 +408,7 @@ def compute_face_edge_tangents(self, dim, entity_id): edge_ts.append(vert_coords[dest] - vert_coords[source]) return edge_ts - def make_points(self, dim, entity_id, order): + def make_points(self, dim, entity_id, order, variant=None): """Constructs a lattice of points on the entity_id:th facet of dimension dim. Order indicates how many points to include in each direction.""" @@ -403,9 +418,9 @@ def make_points(self, dim, entity_id, order): entity_verts = \ self.get_vertices_of_subcomplex( self.get_topology()[dim][entity_id]) - return make_lattice(entity_verts, order, 1) + return make_lattice(entity_verts, order, 1, variant=variant) elif dim == self.get_spatial_dimension(): - return make_lattice(self.get_vertices(), order, 1) + return make_lattice(self.get_vertices(), order, 1, variant=variant) else: raise ValueError("illegal dimension") @@ -1266,7 +1281,9 @@ def make_affine_mapping(xs, ys): def default_simplex(spatial_dim): """Factory function that maps spatial dimension to an instance of the default reference simplex of that dimension.""" - if spatial_dim == 1: + if spatial_dim == 0: + return Point() + elif spatial_dim == 1: return DefaultLine() elif spatial_dim == 2: return DefaultTriangle() diff --git a/setup.py b/setup.py index b4cb69674..fa5a47dc4 100755 --- a/setup.py +++ b/setup.py @@ -27,4 +27,4 @@ download_url=tarball, license="LGPL v3 or later", packages=["FIAT"], - install_requires=["numpy", "sympy"]) + install_requires=["numpy", "sympy", "recursivenodes"]) diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 463c5fc3c..05efd3776 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -247,9 +247,21 @@ def __init__(self, a, b): "GaussLegendre(I, 0)", "GaussLegendre(I, 1)", "GaussLegendre(I, 2)", + "GaussLegendre(T, 0)", + "GaussLegendre(T, 1)", + "GaussLegendre(T, 2)", + "GaussLegendre(S, 0)", + "GaussLegendre(S, 1)", + "GaussLegendre(S, 2)", "GaussLobattoLegendre(I, 1)", "GaussLobattoLegendre(I, 2)", "GaussLobattoLegendre(I, 3)", + "GaussLobattoLegendre(T, 1)", + "GaussLobattoLegendre(T, 2)", + "GaussLobattoLegendre(T, 3)", + "GaussLobattoLegendre(S, 1)", + "GaussLobattoLegendre(S, 2)", + "GaussLobattoLegendre(S, 3)", "Bubble(I, 2)", "Bubble(T, 3)", "Bubble(S, 4)", diff --git a/test/unit/test_gauss_legendre.py b/test/unit/test_gauss_legendre.py index 4fc27863f..849f4c58c 100644 --- a/test/unit/test_gauss_legendre.py +++ b/test/unit/test_gauss_legendre.py @@ -23,25 +23,66 @@ import numpy as np -@pytest.mark.parametrize("degree", range(1, 7)) -def test_gl_basis_values(degree): +def symmetric_simplex(dim): + from FIAT.reference_element import ufc_simplex + s = ufc_simplex(dim) + if dim == 1: + s.vertices = [(-1.,), (1.,)] + elif dim == 2: + h = 3.**0.5 / dim + s.vertices = [(0., 1.), (-h, -0.5), (h, -0.5)] + elif dim == 3: + h = 3.**0.5 / dim + s.vertices = [(-h, h, h), (h, -h, h), (h, h, -h), (h, h, h)] + return s + + +@pytest.mark.parametrize("degree", range(0, 8)) +@pytest.mark.parametrize("dim", (1, 2, 3)) +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 + from FIAT import GaussLegendre, make_quadrature - s = ufc_simplex(1) + s = symmetric_simplex(dim) q = make_quadrature(s, degree + 1) - fe = GaussLegendre(s, degree) - tab = fe.tabulate(0, q.pts)[(0,)] + tab = fe.tabulate(0, q.pts)[(0,)*dim] for test_degree in range(degree + 1): - coefs = [n(lambda x: x[0]**test_degree) for n in fe.dual.nodes] + v = lambda x: sum(x)**test_degree + coefs = [n(v) 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) + reference = np.dot([v(x) for x in q.pts], q.wts) assert np.allclose(integral, reference, rtol=1e-14) +@pytest.mark.parametrize("dim, degree", [(1, 4), (2, 4), (3, 4)]) +def test_edge_dofs(dim, degree): + """ Ensure edge DOFs are point evaluations at GL points.""" + from FIAT import GaussLegendre, quadrature, expansions + + 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 + line = s if dim == 1 else s.construct_subelement(1) + lr = quadrature.GaussLegendreQuadratureLineRule(line, 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)))) + + 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 a2beb454c..adf933d3e 100644 --- a/test/unit/test_gauss_lobatto_legendre.py +++ b/test/unit/test_gauss_lobatto_legendre.py @@ -23,25 +23,67 @@ import numpy as np -@pytest.mark.parametrize("degree", range(1, 7)) -def test_gll_basis_values(degree): +def symmetric_simplex(dim): + from FIAT.reference_element import ufc_simplex + s = ufc_simplex(dim) + if dim == 1: + s.vertices = [(-1.,), (1.,)] + elif dim == 2: + h = 3.**0.5 / dim + s.vertices = [(0., 1.), (-h, -0.5), (h, -0.5)] + elif dim == 3: + h = 3.**0.5 / dim + s.vertices = [(-h, h, h), (h, -h, h), (h, h, -h), (h, h, h)] + return s + + +@pytest.mark.parametrize("degree", range(1, 8)) +@pytest.mark.parametrize("dim", (1, 2, 3)) +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 + from FIAT import GaussLobattoLegendre, make_quadrature - s = ufc_simplex(1) + s = symmetric_simplex(dim) q = make_quadrature(s, degree + 1) - fe = GaussLobattoLegendre(s, degree) - tab = fe.tabulate(0, q.pts)[(0,)] + tab = fe.tabulate(0, q.pts)[(0,)*dim] for test_degree in range(degree + 1): - coefs = [n(lambda x: x[0]**test_degree) for n in fe.dual.nodes] + v = lambda x: sum(x)**test_degree + coefs = [n(v) 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) + reference = np.dot([v(x) for x in q.pts], q.wts) assert np.allclose(integral, reference, rtol=1e-14) +@pytest.mark.parametrize("dim, degree", [(1, 4), (2, 4), (3, 4)]) +def test_edge_dofs(dim, degree): + """ Ensure edge DOFs are point evaluations at GL points.""" + from FIAT import GaussLobattoLegendre, quadrature, expansions + + 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 + line = s if dim == 1 else s.construct_subelement(1) + lr = quadrature.GaussLobattoLegendreQuadratureLineRule(line, 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)))) + + if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__))