Recursive, parameter-free, explicitly defined interpolation nodes for simplices

A rule for constructing interpolation nodes for $n$th degree polynomials on the simplex is presented. These nodes are simple to define recursively from families of 1D node sets, such as the Lobatto-Gauss-Legendre (LGL) nodes. The resulting nodes have attractive properties: they are fully symmetric, they match the 1D family used in construction on the edges of the simplex, and the nodes constructed for the $(d-1)$-simplex are the boundary traces of the nodes constructed for the $d$-simplex. When compared using the Lebesgue constant to other explicit rules for defining interpolation nodes, the nodes recursively constructed from LGL nodes are nearly as good as the"warp&blend"nodes of Warburton in 2D (which, though defined differently, are very similar), and in 3D are better than other known explicit rules by increasing margins for $n>6$. By that same measure, these recursively defined nodes are not as good as implicitly defined nodes found by optimizing the Lebesgue constant or related functions, but such optimal node sets have yet to be computed for the tetrahedron. A reference python implementation has been distributed as the `recursivenodes` package, but the simplicity of the recursive construction makes them easy to implement.


Definition of the recursive rule
The motivating example for this work is the use of Lagrange polynomials as shape functions for the finite element approximation space P n (∆ d ): polynomials of degree at most n on the d-simplex. A Lagrange polynomial basis Φ X = {ϕ i } ⊂ P n (∆ d ) is defined by a set of interpolation nodes X = {x i } ⊂ ∆ d as Φ X = {ϕ i ∈ P n (∆ d ) : ϕ i (x j ) = δ ij }. While some of the properties of an implementation of the finite element method depend only on the approximation space, the basis used can affect the convergence, numerical stability, and computational efficiency. Discussion of each of these aspects follows the definition of the interpolation nodes that are the main contribution of this work.
The nodes are dimensionally recursive, building from points on the interval [0, 1]. A 1D node set is a set of points X n = {x n,i } n i=0 ⊂ [0, 1] that is increasing and symmetric about 1/2, x n,i = 1 − x n,n−i . A 1D node family is a collection X = {X n } n∈N0 . Examples include equispaced nodes, symmetric Gauss-Jacobi quadrature nodes, and symmetric Lobatto-Gauss-Jacobi quadrature nodes. This work uses the standard notation |α| = i α i , and further defines: • #x as the length of a tuple (multi-index or vector), • x \i as the tuple formed by removing the ith element, and • x +i as the augmentation of a tuple by inserting a zero for the ith element.
Given a 1D node family X, the recursive definition of the interpolation node b X (α) ∈ ∆ #α−1 bary is b X (α) =    (1), #α = 1, The full d-simplex node set is and the full d-simplex node family is Unless otherwise specified, the 1D node family X is taken to be the Lobatto-Gauss-Legendre (LGL) family X LGL . Some examples are illustrated in Fig. 1.  ((1, 2, 3)). If X is the equispaced 1D node family (left), the projections meet at an equispaced node in the triangle. If X is the LGL 1D node family (right), the lines do not actually meet at one point. From this comes the idea that, if a good family of interpolation nodes on the (d − 1)-simplex are already known, a heuristic for locating the node b X (α) in the d-simplex is to choose a point whose projection onto each of the facets is one of those good nodes,

Intuition behind the recursive rule
Unfortunately, this is an overdetermined set of requirements.
The system (4) has a solution if X is the family of equispaced nodes, and the solution is an equispaced node in the triangle, as seen in Fig. 2 (left).
Any point in the interior of the triangle is in the convex hull of its projections onto the edges, so if a node location does satisfy (4), then it can be expressed as a barycentric combination of its projections. The equispaced nodes of the triangle not only have projections that are equispaced nodes on the edges, but their barycentric weights have a remarkable property. In Proposition 2.1, the barycentric weights describing equispaced nodes in the triangle are themselves 1D equispaced nodes on the edge (see Fig. 3, left). By analogy, a heuristic for approximating a solution to the overdetermined system (4) is to use the same 1D node family X that was used for the projection points as barycentric weights for combining them (see Fig. 3, right), which restates (1).

Comparison to other node families
The recursive rule (1) generates node families R d X for the d-simplex in each dimension. This section compares them to other node families with respect to several metrics that are relevant to finite element computations.

Boundary and symmetry properties
The nodes R d X,n have three non-numerical properties that make them convenient to use when implementing the finite element method.
I. Symmetry: The symmetry group of the d-simplex is the group S d+1 : for ∆ d bary , each symmetry corresponds to a permutation of the coordinates. It is clear that the recursive rule (1) respects these symmetries, that b X (σ(α)) = σ(b X (α)). This is useful when a d-simplex is viewed from multiple orientations, such as when it is the interface between cells.
In other words, R 1 X,n is the 1D node set X n mapped to the barycentric line ∆ 1 bary . III. Recursive boundary traces: Problems solved by the finite element method can have forms computed over surfaces: data for Neumann boundary conditions or jump terms in discontinuous Galerkin methods, for example. A good node set should induce good shape functions for P n (∆ d ), but also for the trace spaces on the boundary facets, which are embeddings of P n (∆ d−1 ).
The following proposition show that if the 1D node family X has nodes at the endpoints, then R d X,n has nodes on each boundary facet of ∆ d bary that are the R d−1 X,n nodes mapped onto that facet, and so they are appropriate for defining Lagrange polynomials on the trace space.
Now assume the property holds if #α ≤ d, and let #α = d + 1. If i = j, then α \i has a zero at an index ∈ {j, j − 1}, so The order can be switched: So by relabeling and using (1) this implies From this equality (5) simplifies: Properties (II) and (III) together mean that the nodes of R d X,n on an edge are always mappings of the 1D node set X n . This is useful when simplices appear in hybrid meshes with tensor-product cells, which often use tensor products of 1D node sets, because common edges between the two cell types will have the same nodes.

Interpolation properties
A problem discretized by the finite element method may require the approximation of an arbitrary function f in P n (∆ d ). Certain problems have optimal projection operators for this purpose, such as L 2 projection or H 1 projection, but these operators can only be approximated with numerical integration rules, and may be implicit or expensive. When Lagrange polynomials are used as a basis, interpolation at the nodes is an appealing projection onto P n (∆ d ), because it requires the minimum number of function evaluations. Let I X : B(∆ d ) → P n (∆ d ) be the interpolation operator defined by nodes X acting on bounded, measurable functions on the d-simplex. The interpolation error can be bounded by where Λ max (X) is the Lebesgue constant, defined by the shape functions Φ X associated with X, Lebesgue constants for R d X are compared against some other node families on the triangle in Fig. 4 and on the tetrahedron in Fig. 5. These include: • equispaced: Equispaced nodes, defined by b eq,i (α) = α i /|α|. • BLP: The nodes of Blyth, Luo, & Pozrikidis [BP06], [LP06], which like the recursively defined nodes are based on the LGL nodes X LGL . If α > 0, which indicates that the node will be in the interior of the simplex, they are defined by Points on the boundary are mapped from the same rule applied to the The nodes of Warburton [War06], which define the node location b wb (α) as the image of the equispaced node b eq (α) under a smooth bijection of the d-simplex. The bijection sends equispaced nodes to LGL nodes on the edges. The smooth map is nearly isoparametric, but a blending parameter is introduced that controls the distortion in the interior of the element, and optimal values of this blending parameter have been computed for n up to 15 in d = 2 and 3. and Vianello [RSV12] by numerical minimization of Λ max , subject to the constraints that the nodes remain symmetric and that the nodes on the edges be LGL nodes. • CB: Nodes for the tetrahedron computed by Chen and Babuška [CB96] by numerical minimization of the related interpolation metric • HT: Nodes for the tetrahedron computed by Hesthaven and Teng [HT00] as the equilibrium distribution of charged particles.
All of these node families except Roth-leb and RSV-leb are symmetric for all n, and all except equispaced, Roth- 2D: In 2D, the R 2 X node family has Lebesgue constants that are not much worse than those for node families implicitly defined to minimize the Lebesgue Figure 6: The LEBGLS nodes of Rapetti, Sommariva, and Vianello [RSV12], showing the abrupt change in layout between n = 9 (left) and n = 10 (right). constant for n ≤ 9 (Fig. 4, left). For n ≥ 10, the Lebesgue constant grows much faster for R 2 X than for the best implicitly defined nodes. Not coincidentally, at n = 10 the layout of implicitly defined nodes that minimize the Lebesgue constant changes significantly. Until then, the RSV-lebgls nodes look "lattice-like," as though they have been smoothly, symmetrically, and monotonically mapped from the equispaced nodes, the same as the explicit node families. At n = 10, however, this pattern changes (Fig. 6). This suggests that no node family that retains the lattice-like structure, including R d X , can attain a slow growth of Λ max like the implicitly defined families.
In comparison to the other explicitly defined nodes (Fig. 4, right), the R 2 X family is nearly as good as the warp & blend family, which has the best Lebesgue constants: Λ max is never more than 10% different between them for n ≤ 15. In fact, despite the differences in their definitions-warp & blend by continuous bijections, R 2 X by recursion-the node families are remarkably similar for n ≤ 15: b X (α) − b wb (α) ≤ 0.01 for every node in these node sets.

3D:
In 3D there are no published examples of Λ max -optimal node sets that have been numerically computed in the same way as in 2D. Λ max (X) is a nonconvex function of the node coordinates in X, and the number of coordinates grows cubically with n, so this is a challenging optimization problem. Instead, the implicitly defined node families CB and HT optimize simpler objectives: the Λ 2 interpolation metric and the electrostatic potential, respectively, and these have only been computed to n ≤ 9. There is little difference in Λ max between R 3 X and these two families (Fig. 5, left), though it is slightly smaller than both for n ≥ 6.
In comparison to the explicitly defined node families BLP and warp & blend (Fig. 5, right), there is little difference for n ≤ 6 (all are within 7% of each other), but R 3 X is increasingly superior for n ≥ 7. For n = 15, the largest for which the warp & blend nodes' blending parameter has been optimized, the Lebesgue constant of R 3 X is 40% smaller.

Asymptotic interpolation properties
A node family is good for approximation by interpolation if the interpolants are known to converge for a large class of functions. In particular, if f is analytic in the neighborhood of ∆ d , then there is a sequence of polynomials p n ⇒ f (converging uniformly on ∆ d ), so it is possible that given the right node family X = {X n } that I Xn f ⇒ f for all analytic f as well.
The weakest known sufficient condition that guarantees this for d > 1 is subexponential growth of the Lebesgue constant: if Λ max (X n ) 1/n → 1, then I Xn f ⇒ f for f analytic in a neighborhood of ∆ d [Blo+92]. The values of Λ max (R d X,n ) that appeared in Figs. 4 & 5 are tabulated in Table 1. They show no evidence of sub-exponential growth.

Finite element matrix conditioning
Matrices that show up repeatedly in applications of the finite element method include the mass matrix M ij = ∆ d ϕ i ϕ j dx and the stiffness matrix K ij = ∆ d ∇ϕ n,i · ∇ϕ n,j dx. For more general problems-nonlinear or with variable coefficients-there may be a secondary set of points Q = {q i } ⊂ ∆ d , and shape functions and their derivatives must be evaluated at these points. Let . Q could be a set of quadrature points or, in nodal methods, the interpolation nodes themselves, in which case the superscript is dropped. The condition numbers of these matrices (using the definition κ 2 (A) = A 2 A † 2 ) for R d X are compared against the condition numbers for the equispaced, BLP, warp & blend, and RSV-lebgls nodes in Fig. 7. The condition numbers of K, G, and L depend on the choice of reference simplex: in this work, they are computed with respect to the biunit simplex ∆ The rankings of the node families by these metrics are essentially the same as by the Lebesgue constant in Section 3.2.   Table 3: Finite element matrix condition numbers R 3 X n κ 2 (M ) κ 2 (M ) 1/n κ 2 (K) κ 2 (K) 1/n κ 2 (G) κ 2 (G) 1/n κ 2 (L) κ 2 (L) In Tables 2 and 3 the growth rates of these condition numbers can be assessed. The values of κ 2 (M ) 1/n , κ 2 (K) 1/n and κ 2 (L) 1/n are not monotonically decreasing in both dimensions for values of n that have been calculated, which suggests super-exponential growth. κ 2 (G) 1/n appears to be monotonically decreasing towards some limit γ d > 1 for d = 2, and d = 3, which suggests exponential growth, but there is no proof of this fact.

Finite element matrix efficiency
To evaluate the basis functions of R d X,n at a set of nodes Q, one can compute the Vandermonde matrices V R d X,n and V Q with respect to a stable basis for R n (∆ d ), such as the Proriol-Koornwinder-Dubiner basis [Pro57], [Koo75], [Dub91], and . Assuming |Q| ∈ Θ(n d ), the cost of constructing this matrix is Θ(n 3d ). There is no structure in R d X,n that would allow for fast application to a vector of nodal coefficients, so the cost of a matrix vector product is Θ(n 2d ). The same costs hold for each directional derivative of the basis functions.
There appears to be no Lagrange polynomial basis for R n (∆ d ) that improves on this for d > 1, so all of the node families discussed above are equal with respect to this metric. It must be noted, however, that the Bernstein-Bézier basis has fast algorithms that allow for optimal construction in Θ(n 2d ) and application in Θ(n d ) for these matrices, for constant coefficients matrices without quadrature [Kir10], for evaluation at the Stroud quadrature points [AAD11], and for the inverse of the mass matrix [Kir16] (albeit with condition numbers that are worse even than equispaced nodes).

Ease of computation and implementation
Implicitly defined node families, including Roth-leb, RSV-leb, RSV-lebgls, CB, and HT from Section 3.2 and others not discussed, require the solution of an optimization problem over the choice of node coordinates, a problem size that, even with symmetries enforced, is Θ(n d ) in n. Objective functions like Λ max (X) are quite nonconvex, so care must be taken to avoid local minima. It is fair to characterize these node sets as relatively expensive to compute from scratch.
The ease of implementing node sets from a node family is distinct from the computational complexity of computing the node sets from scratch. Most of the implicitly defined node sets discussed in this work have published node sets for moderate values of n for d = 2 [CB95][Hes98] [RSV12], and a few for d = 3 [CB96] [HT00]. If those values encompass the needs of an application, then implementation is as simple as copying and formatting the published node sets.
Of the explicitly defined node families discussed in this work, the equispaced and BLP nodes are the cheapest to compute: the former requires Θ(d) operations and Θ(1) workspace, the latter requires Θ(d 2 ) operations and Θ(1) workspace per node. The warp & blend nodes additionally require d + 1 evaluations of 1D Jacobi polynomials up to degree n at each node (one per facet of the simplex) and one Θ(n 3 ) inversion of a 1D Vandermonde matrix of size n + 1 per node set.
If the nodes must be computed for higher dimensions, the cost of computing a full node set can be reduced by caching the lower-dimensional nodes. Then the cost of computing the node sets {R d X,i } n i=0 with caching satisfies the recursion , leading to an amortized cost per node that is O((n+d)d 2 /n). The workspace with caching satisfies the recursion S(d, which is an amortized space per node that is O((n + d)d/n). In terms of implementation, once the 1D node family X is available, the code to compute b X (α) is very short. Here is an example implementation in python: import numpy as np def recursive(alpha, family): The barycentric d-simplex coordinates for a multi-index alpha with length d+1 and sum n, based on a 1D node family.

Using 1D families other than LGL nodes
The analysis and comparison in Section 3 was conducted under the assumption that the 1D node family X was X LGL , the Lobatto-Gauss-Legendre nodes on the interval [0, 1]. The recursive rule (1) allows for an arbitrary 1D node family. While X LGL appears to be the best choice according to the metrics in Section 3, for completeness a few alternate choices of X are presented here.
X eq (equispaced): It is not surprising, given the discussion in Section 2 where equispaced nodes provided the intuition behind the recursive rule, that using X = X eq reproduces the equispaced nodes, b Xeq (α) = b eq (α).

X LGC (Lobatto-Gauss-Chebyshev):
The 1D Lobatto-Gauss-Chebyshev (LGC) node family has good interpolation properties in 1D while being a nested family, with X n,LGC ⊂ X 2n,LGC . The recursive nodes R d X LGC inherit the nested property, as demonstrated in Fig. 8, left.
X GL (Gauss-Legendre): The 1D Gauss-Legendre (GC) node family does not include the endpoints 0 and 1, but can still be used to construct nodes. Property (III) in Section 3.1 does not hold, and in fact all nodes will be in the interior of the d-simplex as demonstrated in Fig. 8, right. The metrics from Section 3 are used to compare R d X for X LGL , X LGC , and X GL in Fig. 9. 2 The results for d = 2 and d = 3 resemble the results of the 1D node families, with GL nodes having worse interpolation properties but better conditioned mass matrices than either set of Lobatto nodes, and with LGC nodes having interpolation properties and matrix condition numbers that similar but slightly worse to LGL nodes. The most interesting trend to be observed in Fig. 9 is that, while the growth rate of the Lebesgue constant for the GL nodes is worse than the Lobatto nodes in 2D, it is much closer in 3D. LGL LGC GL Figure 9: Comparing R d X LGL , R d X LGC , and R d X GL according to metrics from Section 3. 16

Conclusion
How and by whom should the nodes R d X defined by the recursive rule (1) be used? The comparisons in this paper have made the case that is the best explicit construction rule thus far, because of its simplicity, its symmetry, and its performance in the metrics that matter to finite element construction (but not for producing asymptotically convergent interpolants). It does not outperform the Warburton's warp & blend node family in 2D, so software already using those would not benefit from switching, but its performance is superior to all other explicit node families in 3D, particularly for n ≥ 7. Likewise, where implicitly defined node families-such as Rapetti, Sommariva, and Vianello's LEBGLS nodes-have been computed and published, they are superior to the R d X node family, especially in 2D for n ≥ 10. But at the time of this writing the tetrahedron has not received nearly as much attention as the triangle, and so this new node family is the best available in 3D.