Add Gauss-Radau points and a corresponding finite element.

FIAT/__init__.py

@@ -21,6 +21,7 @@
 from FIAT.lagrange import Lagrange
 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre
 from FIAT.gauss_legendre import GaussLegendre
+from FIAT.gauss_radau import GaussRadau
 from FIAT.morley import Morley
 from FIAT.nedelec import Nedelec
 from FIAT.nedelec_second_kind import NedelecSecondKind

FIAT/gauss_radau.py

@@ -0,0 +1,36 @@
+# 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 Robert C. Kirby ([email protected]), 2020
+
+
+from FIAT import finite_element, polynomial_set, dual_set, functional, quadrature
+from FIAT.reference_element import LINE
+
+
+class GaussRadauDualSet(dual_set.DualSet):
+    """The dual basis for 1D discontinuous elements with nodes at the
+    Gauss-Radau points."""
+    def __init__(self, ref_el, degree, right=True):
+        # Do DG connectivity because it's bonkers to do one-sided assembly even
+        # though we have an endpoint in the point set!
+        entity_ids = {0: {0: [], 1: []},
+                      1: {0: list(range(0, degree+1))}}
+        lr = quadrature.RadauQuadratureLineRule(ref_el, degree+1, right)
+        nodes = [functional.PointEvaluation(ref_el, x) for x in lr.pts]
+
+        super(GaussRadauDualSet, self).__init__(nodes, ref_el, entity_ids)
+
+
+class GaussRadau(finite_element.CiarletElement):
+    """1D discontinuous element with nodes at the Gauss-Radau points."""
+    def __init__(self, ref_el, degree):
+        if ref_el.shape != LINE:
+            raise ValueError("Gauss-Radau elements are only defined in one dimension.")
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
+        dual = GaussRadauDualSet(ref_el, degree)
+        formdegree = ref_el.get_spatial_dimension()  # n-form
+        super(GaussRadau, self).__init__(poly_set, dual, degree, formdegree)

FIAT/quadrature.py

@@ -123,6 +123,63 @@ def __init__(self, ref_el, m):
         QuadratureRule.__init__(self, ref_el, xs, ws)
 
 
+class RadauQuadratureLineRule(QuadratureRule):
+    """Produce the Gauss--Radau quadrature rules on the interval using
+    an adaptation of Winkel's Matlab code. 
+
+    The quadrature rule uses m points for a degree of precision of 2m-1.
+    """
+    def __init__(self, ref_el, m, right=True):
+        assert m >= 1
+        N = m - 1
+        # Use Chebyshev-Gauss-Radau nodes as initial guess for LGR nodes
+        x=-numpy.cos(2*numpy.pi*numpy.linspace(0, N, m)/(2*N+1))
+
+        P = numpy.zeros((N+1,N+2))
+
+        xold = 2
+
+        free=numpy.arange(1, N+1, dtype='int')
+
+        while numpy.max(numpy.abs(x-xold))>5e-16:
+            xold=x.copy()
+
+            P[0,:]=(-1)**numpy.arange(0, N+2)
+            P[free,0]=1
+            P[free,1]=x[free]
+
+            for k in range(2, N+2):
+                P[free,k]=( (2*k-1)*x[free]*P[free,k-1]-(k-1)*P[free,k-2] )/k
+
+            x[free]=xold[free]-((1-xold[free])/(N+1))*(P[free,N]+P[free,N+1])/(P[free,N]-P[free,N+1])
+
+        # The Legendre-Gauss-Radau Vandermonde
+        P=P[:,:-1]
+        # Compute the weights
+        w=numpy.zeros(N+1)
+        w[0]=2/(N+1)**2
+        w[free]=(1-x[free])/((N+1)*P[free,-1])**2
+
+        if right:
+            x = numpy.flip(-x)
+            w = numpy.flip(w)
+
+        xs_ref = x
+        ws_ref = w
+
+        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)
+
+        xs = tuple([tuple(mapping(x_ref)[0]) for x_ref in xs_ref])
+        ws = tuple([scale * w for w in ws_ref])
+
+        QuadratureRule.__init__(self, ref_el, xs, ws)
+
+
 class CollapsedQuadratureTriangleRule(QuadratureRule):
     """Implements the collapsed quadrature rules defined in
     Karniadakis & Sherwin by mapping products of Gauss-Jacobi rules