GLL/GL on simplices with recursive points

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,61 @@
 #
 # 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,
+                  quadrature, recursive_points)
+from FIAT.reference_element import POINT, LINE, TRIANGLE, TETRAHEDRON, UFCInterval
 from FIAT.orientation_utils import make_entity_permutations_simplex
 from FIAT.barycentric_interpolation import LagrangePolynomialSet
 
 
+class GaussLegendrePointSet(recursive_points.RecursivePointSet):
+    """Recursive point set on simplices based on the Gauss-Legendre points on
+    the interval"""
+    def __init__(self):
+        ref_el = UFCInterval()
+        lr = quadrature.GaussLegendreQuadratureLineRule
+        f = lambda n: lr(ref_el, n + 1).get_points()
+        super(GaussLegendrePointSet, self).__init__(f)
+
+
 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."""
+    point_set = GaussLegendrePointSet()
+
     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 = self.point_set.recursive_points(ref_el.get_vertices(), degree)
+        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:
+            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,69 @@
 #
 # 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,
+                  quadrature, recursive_points)
+from FIAT.reference_element import LINE, TRIANGLE, TETRAHEDRON, UFCInterval
 from FIAT.orientation_utils import make_entity_permutations_simplex
 from FIAT.barycentric_interpolation import LagrangePolynomialSet
 
 
+class GaussLobattoLegendrePointSet(recursive_points.RecursivePointSet):
+    """Recursive point set on simplices based on the Gauss-Lobatto points on
+    the interval"""
+    def __init__(self):
+        ref_el = UFCInterval()
+        lr = quadrature.GaussLobattoLegendreQuadratureLineRule
+        f = lambda n: lr(ref_el, n + 1).get_points() if n else None
+        super(GaussLobattoLegendrePointSet, self).__init__(f)
+
+
 class GaussLobattoLegendreDualSet(dual_set.DualSet):
-    """The dual basis for 1D continuous elements with nodes at the
-    Gauss-Lobatto points."""
+    """The dual basis for continuous elements with nodes at the
+    (recursive) Gauss-Lobatto points."""
+    point_set = GaussLobattoLegendrePointSet()
+
     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_ids = {}
+        nodes = []
         entity_permutations = {}
-        entity_permutations[0] = {0: {0: [0]}, 1: {0: [0]}}
-        entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree - 1)}
+
+        # make nodes by getting points
+        # need to do this dimension-by-dimension, facet-by-facet
+        top = ref_el.get_topology()
+
+        cur = 0
+        for dim in sorted(top):
+            entity_ids[dim] = {}
+            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 = self.point_set.make_points(ref_el, dim, entity, degree)
+                nodes_cur = [functional.PointEvaluation(ref_el, x)
+                             for x in pts_cur]
+                nnodes_cur = len(nodes_cur)
+                nodes += nodes_cur
+                entity_ids[dim][entity] = list(range(cur, cur + nnodes_cur))
+                cur += nnodes_cur
+                entity_permutations[dim][entity] = perms
 
         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.")
+        if ref_el.shape not in {LINE, TRIANGLE, TETRAHEDRON}:
+            raise ValueError("Gauss-Lobatto-Legendre elements are only defined on simplices.")
         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 == LINE:
+            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/polynomial_set.py

@@ -18,6 +18,9 @@
 import numpy
 from FIAT import expansions
 from FIAT.functional import index_iterator
+from FIAT.reference_element import UFCInterval
+from FIAT.quadrature import GaussLegendreQuadratureLineRule
+from FIAT.recursive_points import RecursivePointSet
 
 
 def mis(m, n):
@@ -125,6 +128,7 @@ class ONPolynomialSet(PolynomialSet):
     for vector- and tensor-valued sets as well.
 
     """
+    point_set = RecursivePointSet(lambda n: GaussLegendreQuadratureLineRule(UFCInterval(), n + 1).get_points())
 
     def __init__(self, ref_el, degree, shape=tuple()):
 
@@ -155,23 +159,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 = self.point_set.recursive_points(ref_el.get_vertices(), degree)
 
             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/recursive_points.py

@@ -0,0 +1,166 @@
+# Copyright (C) 2015 Imperial College London and others.
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+#
+# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2023
+
+import numpy
+
+"""
+@article{isaac2020recursive,
+  title={Recursive, parameter-free, explicitly defined interpolation nodes for simplices},
+  author={Isaac, Tobin},
+  journal={SIAM Journal on Scientific Computing},
+  volume={42},
+  number={6},
+  pages={A4046--A4062},
+  year={2020},
+  publisher={SIAM}
+}
+"""
+
+
+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,)
+
+
+class RecursivePointSet(object):
+    """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, family):
+        self._family = family
+        self._cache = {}
+
+    def interval_points(self, degree):
+        try:
+            return self._cache[degree]
+        except KeyError:
+            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
+        multiindex alpha of length d and sum n, based on a 1D node family."""
+        d = len(alpha)
+        n = sum(alpha)
+        b = numpy.zeros((d,), dtype="d")
+        xn = self.interval_points(n)
+        if xn is None:
+            return b
+        if d == 2:
+            b[:] = xn[list(alpha)]
+            return b
+        weight = 0.0
+        for i in range(d):
+            w = xn[n - alpha[i]]
+            alpha_noti = alpha[:i] + alpha[i+1:]
+            br = self._recursive(alpha_noti)
+            b[:i] += w * br[:i]
+            b[i+1:] += w * br[i:]
+            weight += w
+        b /= weight
+        return b
+
+    def recursive_points(self, vertices, order, interior=0):
+        X = numpy.array(vertices)
+        get_point = lambda alpha: tuple(numpy.dot(self._recursive(alpha), X))
+        return list(map(get_point, multiindex_equal(len(vertices), order, interior)))
+
+    def make_points(self, ref_el, dim, entity_id, order):
+        """Constructs a lattice of points on the entity_id:th
+        facet of dimension dim.  Order indicates how many points to
+        include in each direction."""
+        if dim == 0:
+            return (ref_el.get_vertices()[entity_id], )
+        elif 0 < dim < ref_el.get_spatial_dimension():
+            entity_verts = \
+                ref_el.get_vertices_of_subcomplex(
+                    ref_el.get_topology()[dim][entity_id])
+            return self.recursive_points(entity_verts, order, 1)
+        elif dim == ref_el.get_spatial_dimension():
+            return self.recursive_points(ref_el.get_vertices(), order, 1)
+        else:
+            raise ValueError("illegal dimension")
+
+
+if __name__ == "__main__":
+    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
+    ref_el.vertices = [(0, s), (-1.0, s-h), (1.0, s-h)]
+    x = numpy.array(ref_el.vertices + ref_el.vertices[:1])
+    plt.plot(x[:, 0], x[:, 1], "k")
+
+    for d in range(1, 4):
+        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 = cg_points.make_points(ref_el, dim, entity, order)
+            if len(pts):
+                x = numpy.array(pts)
+                for r in range(1, 3):
+                    th = r * (2*numpy.pi)/3
+                    ct = numpy.cos(th)
+                    st = numpy.sin(th)
+                    Q = numpy.array([[ct, -st], [st, ct]])
+                    x = numpy.concatenate([x, numpy.dot(x, Q)])
+                plt.scatter(x[:, 0], x[:, 1])
+
+    x0 = 2.0
+    h = -h
+    s = 2*h/3
+    ref_el = reference_element.ufc_simplex(2)
+    ref_el.vertices = [(x0, s), (x0-1.0, s-h), (x0+1.0, s-h)]
+
+    x = numpy.array(ref_el.vertices + ref_el.vertices[:1])
+    d = len(ref_el.vertices)
+    x0 = sum(x[:d])/d
+    plt.plot(x[:, 0], x[:, 1], "k")
+
+    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
+        ct = numpy.cos(th)
+        st = numpy.sin(th)
+        Q = numpy.array([[ct, -st], [st, ct]])
+        x = numpy.concatenate([x, numpy.dot(x-x0, Q)+x0])
+    plt.scatter(x[:, 0], x[:, 1])
+
+    plt.gca().set_aspect('equal', 'box')
+    plt.show()