These html pages are based on the PhD thesis "Cluster-Based Parallelization of Simulations on Dynamically Adaptive Grids and Dynamic Resource Management" by Martin Schreiber.
There is also more information and a PDF version available.

2.3 Basis functions

To approximate the conserved quantities in each reference element, a representation of our state U(ξ,η,t) at coordinate ξ and η at time stamp t has to be chosen. We use a basis which is generated by i polynomial functions φi(ξ,η) for each reference element. A linear combination of these polynomials and their associated weights qik(t) for each conserved quantity k then approximates the solution

U k(ξ,η,t) ≈ qk(ξ,η,t) =    qk(t)φi(ξ,η ) = ˆU k(ξ,η,t)
                        i  i
and in vector form
            (∑    1         )
            |∑  iqi2(t)φi(ξ,η)|
U (ξ,η,t) ≈ |   iqi(t)φi(ξ,η)| = Uˆ(ξ,η,t).
            (∑     ...      )

Next, we choose a set of basis functions for the approximation of our solution and introduce the nodal and modal basis functions. These different sets of basis functions have different properties such as their direct applicability to flux computations and their computational intensity, see e.g.  [LNT07]. We introduce nodal and modal basis functions to show the applicability of the developed interfaces for both of them.

2.3.1 Nodal basis functions

Nodal polynomial basis functions are generated for a particular set of nodal points pi T with T the triangle reference space. Evaluating the corresponding basis function at each nodal point holds φi(pj) = δi,j with δ being the Kronecker delta. Thus, the related coefficient qi directly represents the approximated solution at a nodal point pi. Based on a given set of nodal points, the basis functions can then e.g. be computed using a Lagrange interpolation method. For triangles, a manifold of choices for nodal points exists and the locations of nodal points are chosen very carefully. Otherwise, this leads to ill-conditioned or even singular mass matrices [SSD03]. Integrals which are evaluated during or in advance of DG simulations can be computed with numerical integration. Therefore, choosing quadrature points from numerical integration as nodal points would be an obvious choice.

Extending the well known one-dimensional Gaussian quadrature formula to two-dimensions using tensor product would lead to quadrature points outside the reference triangle and more quadrature points would have to be evaluated than necessary. Dunavant [Dun85] developed a distribution of quadrature points with symmetric layout also minimizing the number of quadrature points, however with quadrature points still outside the triangle for higher-order quadrature.

An alternative to these quadrature points is given by Gauss-Lobatto (GL) points [HW08]. In contrast to Dunavant quadrature points, Gauss-Lobatto quadrature points are also placed at the boundaries of the reference triangle and have the advantageous property of creating sparse matrices for flux computations (see Section 2.7). An additionally required property of our nodal basis is a unique representation of the polynomials used as basis, also denoted as unisolvency [SSD03]. This unisolvency can be assured by testing for a non-singular Vandermonde matrix e.g. by inverting it to compute the coefficients for the basis polynomials. This property is also fulfilled by the (GL) points. With the properties of creating sparse matrices for flux computations and non-singular Vandermonde matrices, we focus on using these (GL) nodal points in our simulation.

Independent of the nodal points, we further assume a given degree d of the polynomials which we use for our discretization. Then, the maximum number of basis functions is given by

N :=       2      .

2.3.2 Modal basis functions

Using orthogonal basis functions, this would lead to diagonal mass matrices (derived in Section 2.5). This is stated to lead to a computational efficient way to compute the mandatorily required inverted mass matrix [KDDLPI07]. Whereas one-dimensional Jacobi polynomials provide such an orthogonality on the unit interval, Dubiner presented a required extension to the projection of those polynomials to conserve the orthogonality also for integrals over our reference triangle [Dub91]. Here, n-th Jacobi polynomial Pn(a,b)(ξ) is parameterized with the coefficients a and b and the spatially related component ξ.

Using PQ(n,a,b,ξ) := Pn(a,b)(2ξ - 1), the m and n-th orthogonal Jacobi polynomial on a triangle is then given by

                     (             )
P o  (m,n, ξ,η) := PQ  m, 0,0,--ξ--  (1- η)m ⋅P Q (n,2m + 1,0,η).
   tri                         1-  η

These orthogonal Jacobi polynomials can be further transformed into orthonormal polynomials for a triangle by normalization:

                      P otri(m, n,ξ,η)
Ptri(m, n,ξ,η) :=  ∘-∫----------------------.
                    T P otri(m, n,ξ,η)2d(ξ,η)

Examples are given in Appendix A.1.2.