GLL/GL on simplices with recursive points

.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

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.")

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)

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)

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)

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)
 

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