From 7d418fa0a372ac6f5e103533ab77ad6a9fac764c Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Mon, 19 Apr 2021 10:04:22 +0100
Subject: [PATCH] GL/GLL tabulation via barycentric interpolation (#65)

Implement numerical stable tabulation of (derivatives) of GL/GLL basis functions
using the second barycentric interpolation formula of Berrut and Trefethen (2004).
---
 FIAT/barycentric_interpolation.py | 46 +++++++++++++++++++++++++++++++
 FIAT/gauss_legendre.py            | 23 ++++++++++++++++
 FIAT/gauss_lobatto_legendre.py    | 23 ++++++++++++++++
 3 files changed, 92 insertions(+)
 create mode 100644 FIAT/barycentric_interpolation.py

diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py
new file mode 100644
index 00000000..ceee2f1a
--- /dev/null
+++ b/FIAT/barycentric_interpolation.py
@@ -0,0 +1,46 @@
+# Copyright (C) 2021 Pablo D. Brubeck
+#
+# 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), 2021
+
+import numpy
+
+
+def barycentric_interpolation(xsrc, xdst, order=0):
+    """Return tabulations of a 1D Lagrange nodal basis via the second barycentric interpolation formula
+
+    See Berrut and Trefethen (2004) https://doi.org/10.1137/S0036144502417715 Eq. (4.2) & (9.4)
+
+    :arg xsrc: a :class:`numpy.array` with the nodes defining the Lagrange polynomial basis
+    :arg xdst: a :class:`numpy.array` with the interpolation points
+    :arg order: the integer order of differentiation
+    :returns: dict of tabulations up to the given order (in the same format as :meth:`~.CiarletElement.tabulate`)
+    """
+
+    # w = barycentric weights
+    # D = spectral differentiation matrix (D.T : u(xsrc) -> u'(xsrc))
+    # I = barycentric interpolation matrix (I.T : u(xsrc) -> u(xdst))
+
+    D = numpy.add.outer(-xsrc, xsrc)
+    numpy.fill_diagonal(D, 1.0E0)
+    w = 1.0E0 / numpy.prod(D, axis=0)
+    D = numpy.divide.outer(w, w) / D
+    numpy.fill_diagonal(D, numpy.diag(D) - numpy.sum(D, axis=0))
+
+    I = numpy.add.outer(-xsrc, xdst)
+    idx = numpy.argwhere(numpy.isclose(I, 0.0E0, 1E-14))
+    I[idx[:, 0], idx[:, 1]] = 1.0E0
+    I = 1.0E0 / I
+    I *= w[:, None]
+    I[:, idx[:, 1]] = 0.0E0
+    I[idx[:, 0], idx[:, 1]] = 1.0E0
+    I = (1.0E0 / numpy.sum(I, axis=0)) * I
+
+    derivs = {(0,): I}
+    for k in range(0, order):
+        derivs[(k+1,)] = numpy.matmul(D, derivs[(k,)])
+
+    return derivs
diff --git a/FIAT/gauss_legendre.py b/FIAT/gauss_legendre.py
index 3ed840e0..18f68cbb 100644
--- a/FIAT/gauss_legendre.py
+++ b/FIAT/gauss_legendre.py
@@ -5,9 +5,14 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 #
 # Written by David A. Ham (david.ham@imperial.ac.uk), 2015
+#
+# Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021
+
+import numpy
 
 from FIAT import finite_element, polynomial_set, dual_set, functional, quadrature
 from FIAT.reference_element import LINE
+from FIAT.barycentric_interpolation import barycentric_interpolation
 
 
 class GaussLegendreDualSet(dual_set.DualSet):
@@ -31,3 +36,21 @@ def __init__(self, ref_el, degree):
         dual = GaussLegendreDualSet(ref_el, degree)
         formdegree = ref_el.get_spatial_dimension()  # n-form
         super(GaussLegendre, self).__init__(poly_set, dual, degree, formdegree)
+
+    def tabulate(self, order, points, entity=None):
+        # This overrides the default with a more numerically stable algorithm
+
+        if entity is None:
+            entity = (self.ref_el.get_dimension(), 0)
+
+        entity_dim, entity_id = entity
+        transform = self.ref_el.get_entity_transform(entity_dim, entity_id)
+
+        xsrc = []
+        for node in self.dual.nodes:
+            # Assert singleton point for each node.
+            (pt,), = node.get_point_dict().keys()
+            xsrc.append(pt)
+        xsrc = numpy.asarray(xsrc)
+        xdst = numpy.array(list(map(transform, points))).flatten()
+        return barycentric_interpolation(xsrc, xdst, order=order)
diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py
index cb469c8b..ee10c180 100644
--- a/FIAT/gauss_lobatto_legendre.py
+++ b/FIAT/gauss_lobatto_legendre.py
@@ -5,9 +5,14 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 #
 # Written by David A. Ham (david.ham@imperial.ac.uk), 2015
+#
+# Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2021
+
+import numpy
 
 from FIAT import finite_element, polynomial_set, dual_set, functional, quadrature
 from FIAT.reference_element import LINE
+from FIAT.barycentric_interpolation import barycentric_interpolation
 
 
 class GaussLobattoLegendreDualSet(dual_set.DualSet):
@@ -31,3 +36,21 @@ def __init__(self, ref_el, degree):
         dual = GaussLobattoLegendreDualSet(ref_el, degree)
         formdegree = 0  # 0-form
         super(GaussLobattoLegendre, self).__init__(poly_set, dual, degree, formdegree)
+
+    def tabulate(self, order, points, entity=None):
+        # This overrides the default with a more numerically stable algorithm
+
+        if entity is None:
+            entity = (self.ref_el.get_dimension(), 0)
+
+        entity_dim, entity_id = entity
+        transform = self.ref_el.get_entity_transform(entity_dim, entity_id)
+
+        xsrc = []
+        for node in self.dual.nodes:
+            # Assert singleton point for each node.
+            (pt,), = node.get_point_dict().keys()
+            xsrc.append(pt)
+        xsrc = numpy.asarray(xsrc)
+        xdst = numpy.array(list(map(transform, points))).flatten()
+        return barycentric_interpolation(xsrc, xdst, order=order)