SpECTRE
v2024.09.16

Functionality associated with a particular choice of basis functions and quadrature for spectral operations. More...
Namespaces  
namespace  filtering 
Matrices for filtering spectral coefficients.  
namespace  Swsh 
Namespace for spinweighted spherical harmonic utilities.  
Classes  
struct  MortarSizeHash 
Typedefs  
using  MortarSize = ChildSize 
The portion of an element covered by a mortar.  
Enumerations  
enum class  Basis : uint8_t { Uninitialized = 0 << basis_shift , Chebyshev = 1 << basis_shift , Legendre = 2 << basis_shift , FiniteDifference = 3 << basis_shift , SphericalHarmonic = 4 << basis_shift } 
Either the basis functions used by a spectral or discontinuous Galerkin (DG) method, or the value FiniteDifference when a finite difference method is used. More...  
enum class  ChildSize : uint8_t { Uninitialized = 0 , Full = 1 , UpperHalf = 2 , LowerHalf = 3 } 
The portion of a mesh covered by a child mesh.  
enum class  Quadrature : uint8_t { Uninitialized = 0 , Gauss , GaussLobatto , CellCentered , FaceCentered , Equiangular } 
Either the choice of quadrature method to compute integration weights for a spectral or discontinuous Galerkin (DG) method, or the locations of grid points when a finite difference method is used. More...  
Functions  
std::array< Basis, 5 >  all_bases () 
All possible values of Basis.  
Basis  to_basis (const std::string &basis) 
Convert a string to a Basis enum.  
std::ostream &  operator<< (std::ostream &os, const Basis &basis) 
Output operator for a Basis.  
template<typename CoefficientDataType , typename DataType >  
DataType  evaluate_clenshaw (const std::vector< CoefficientDataType > &coefficients, const DataType &alpha, const DataType &beta, const DataType &f0_of_x, const DataType &f1_of_x) 
Clenshaw's algorithm for evaluating linear combinations of basis functions obeying a given three term recurrence relation. See numerical recipes 3rd edition 5.4.2 for details. alpha and beta define the recursion \( f_n(x) = \alpha f_{n1}(x) + \beta f_{n2}(x) \) f0_of_x and f1_of_x are \( f_0(x)\) and \(f_1(x)\) respectively.  
std::ostream &  operator<< (std::ostream &os, ChildSize mortar_size) 
template<size_t Dim>  
bool  needs_projection (const Mesh< Dim > &mesh1, const Mesh< Dim > &mesh2, const std::array< ChildSize, Dim > &child_sizes) 
Determine whether data needs to be projected between a child mesh and its parent mesh. If no projection is necessary the data may be used asis. Projection is necessary if the child is either prefined or hrefined relative to its parent, or both. This operation is symmetric, i.e. it is irrelevant in which order the child and the parent mesh are passed in.  
const Matrix &  projection_matrix_child_to_parent (const Mesh< 1 > &child_mesh, const Mesh< 1 > &parent_mesh, ChildSize size, bool operand_is_massive=false) 
The projection matrix from a child mesh to its parent. More...  
template<size_t Dim>  
std::array< std::reference_wrapper< const Matrix >, Dim >  projection_matrix_child_to_parent (const Mesh< Dim > &child_mesh, const Mesh< Dim > &parent_mesh, const std::array< ChildSize, Dim > &child_sizes, bool operand_is_massive=false) 
The projection matrix from a child mesh to its parent, in Dim dimensions.  
const Matrix &  projection_matrix_parent_to_child (const Mesh< 1 > &parent_mesh, const Mesh< 1 > &child_mesh, ChildSize size) 
The projection matrix from a parent mesh to one of its children. More...  
template<size_t Dim>  
std::array< std::reference_wrapper< const Matrix >, Dim >  projection_matrix_parent_to_child (const Mesh< Dim > &parent_mesh, const Mesh< Dim > &child_mesh, const std::array< ChildSize, Dim > &child_sizes) 
The projection matrix from a parent mesh to one of its children, in Dim dimensions.  
template<size_t Dim>  
std::array< std::reference_wrapper< const Matrix >, Dim >  p_projection_matrices (const Mesh< Dim > &source_mesh, const Mesh< Dim > &target_mesh) 
The projection matrices from a source mesh to a target mesh where the meshes cover the same physical volume.  
std::array< Quadrature, 6 >  all_quadratures () 
All possible values of Quadrature.  
Quadrature  to_quadrature (const std::string &quadrature) 
Convert a string to a Quadrature enum.  
std::ostream &  operator<< (std::ostream &os, const Quadrature &quadrature) 
Output operator for a Quadrature.  
template<Basis BasisType, typename T >  
T  compute_basis_function_value (size_t k, const T &x) 
Compute the function values of the basis function \(\Phi_k(x)\) (zeroindexed).  
template<Basis >  
DataVector  compute_inverse_weight_function_values (const DataVector &) 
Compute the inverse of the weight function \(w(x)\) w.r.t. which the basis functions are orthogonal. See the description of quadrature_weights(size_t) for details. This is arbitrarily set to 1 for FiniteDifference basis, to integrate using the midpoint method (see quadrature_weights (size_t) for details).  
template<Basis BasisType>  
double  compute_basis_function_normalization_square (size_t k) 
Compute the normalization square of the basis function \(\Phi_k\) (zeroindexed), i.e. the weighted definite integral over its square.  
template<Basis BasisType, Quadrature QuadratureType>  
std::pair< DataVector, DataVector >  compute_collocation_points_and_weights (size_t num_points) 
Compute the collocation points and weights associated to the basis and quadrature. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const DataVector &  collocation_points (size_t num_points) 
Collocation points. More...  
const DataVector &  collocation_points (const Mesh< 1 > &mesh) 
Collocation points for a onedimensional mesh. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const DataVector &  quadrature_weights (size_t num_points) 
Weights to compute definite integrals. More...  
const DataVector &  quadrature_weights (const Mesh< 1 > &mesh) 
Quadrature weights for a onedimensional mesh. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const Matrix &  weak_flux_differentiation_matrix (size_t num_points) 
Matrix used to compute the divergence of the flux in weak form. More...  
const Matrix &  weak_flux_differentiation_matrix (const Mesh< 1 > &mesh) 
Matrix used to compute the divergence of the flux in weak form. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const Matrix &  integration_matrix (size_t num_points) 
Matrix used to perform an indefinite integral of a function over the logical grid. The left boundary condition is such that the integral is 0 at \(\xi=1\) More...  
const Matrix &  integration_matrix (const Mesh< 1 > &mesh) 
Indefinite integration matrix for a onedimensional mesh. More...  
template<Basis BasisType, Quadrature QuadratureType, typename T >  
Matrix  interpolation_matrix (size_t num_points, const T &target_points) 
Matrix used to interpolate to the target_points . More...  
template<typename T >  
Matrix  interpolation_matrix (const Mesh< 1 > &mesh, const T &target_points) 
Interpolation matrix to the target_points for a onedimensional mesh. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const Matrix &  modal_to_nodal_matrix (size_t num_points) 
Matrix used to transform from the spectral coefficients (modes) of a function to its nodal coefficients. Also referred to as the Vandermonde matrix. More...  
const Matrix &  modal_to_nodal_matrix (const Mesh< 1 > &mesh) 
Transformation matrix from modal to nodal coefficients for a onedimensional mesh. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const Matrix &  nodal_to_modal_matrix (size_t num_points) 
Matrix used to transform from the nodal coefficients of a function to its spectral coefficients (modes). Also referred to as the inverse Vandermonde matrix. More...  
const Matrix &  nodal_to_modal_matrix (const Mesh< 1 > &mesh) 
Transformation matrix from nodal to modal coefficients for a onedimensional mesh. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const Matrix &  linear_filter_matrix (size_t num_points) 
Matrix used to linearize a function. More...  
const Matrix &  linear_filter_matrix (const Mesh< 1 > &mesh) 
Linear filter matrix for a onedimensional mesh. More...  
template<size_t DimMinusOne>  
size_t  hash (const std::array< Spectral::ChildSize, DimMinusOne > &mortar_size) 
Performs a perfect hash of the mortars into \(2^{d1}\) slots on the range \([0, 2^{d1})\). More...  
template<Basis BasisType, Quadrature QuadratureType>  
const Matrix &  differentiation_matrix (size_t num_points) 
Matrix used to compute the derivative of a function. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const Matrix &  differentiation_matrix_transpose (size_t num_points) 
Matrix used to compute the derivative of a function. More...  
const Matrix &  differentiation_matrix (const Mesh< 1 > &mesh) 
Differentiation matrix for a onedimensional mesh. More...  
const Matrix &  differentiation_matrix_transpose (const Mesh< 1 > &mesh) 
Differentiation matrix for a onedimensional mesh. More...  
const std::pair< Matrix, Matrix > &  boundary_interpolation_matrices (const Mesh< 1 > &mesh) 
Matrices that interpolate to the lower and upper boundaries of the element. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const std::pair< Matrix, Matrix > &  boundary_interpolation_matrices (size_t num_points) 
Matrices that interpolate to the lower and upper boundaries of the element. More...  
const std::pair< DataVector, DataVector > &  boundary_interpolation_term (const Mesh< 1 > &mesh) 
Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhustype boundary conditions in a discontinuous Galerkin scheme using Gauss points. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const std::pair< DataVector, DataVector > &  boundary_interpolation_term (size_t num_points) 
Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhustype boundary conditions in a discontinuous Galerkin scheme using Gauss points. More...  
const std::pair< DataVector, DataVector > &  boundary_lifting_term (const Mesh< 1 > &mesh) 
Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points. More...  
template<Basis BasisType, Quadrature QuadratureType>  
const std::pair< DataVector, DataVector > &  boundary_lifting_term (size_t num_points) 
Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points. More...  
Variables  
constexpr uint8_t  basis_shift = 4 
The amount we leftshift the basis integers by. More...  
template<Basis basis, Quadrature quadrature>  
constexpr size_t  minimum_number_of_points 
Minimum number of possible collocation points for a quadrature type. More...  
template<Basis basis>  
constexpr size_t  maximum_number_of_points 
Maximum number of allowed collocation points. More...  
Functionality associated with a particular choice of basis functions and quadrature for spectral operations.
The functions in this namespace provide lowlevel access to collocation points, quadrature weights and associated matrices, such as for differentiation and interpolation. They are available in two versions: either templated directly on the enum cases of the Spectral::Basis and Spectral::Quadrature types, or taking a onedimensional Mesh as their argument.
Most algorithms in this namespace are adapted from [110].

strong 
Either the basis functions used by a spectral or discontinuous Galerkin (DG) method, or the value FiniteDifference
when a finite difference method is used.
The particular choices of Basis and Quadrature determine where the collocation points of a Mesh are located in an Element. For a spectral or DG method, the Basis also represents the choice of basis functions used to represent a function on an Element, which then provides a convenient choice for the operators used for differentiation, interpolation, etc. For a finite difference method, one needs to choose the order of the scheme (and hence the weights, differentiation matrix, integration weights, and interpolant) locally in space and time to handle discontinuous solutions.
Legendre
for a generalpurpose spectral or DG mesh, unless you have a particular reason for choosing Chebyshev
.SphericalHarmonic
to choose a spherical harmonic basis. By convention, the first dimension represents the polar/zentith angle (or colatitude), while the second dimension represents the azimuthal angle (or longitude)Uninitialized
value. The number of bits to shift is encoded in the variable Spectral::detail::basis_shift
.

strong 
Either the choice of quadrature method to compute integration weights for a spectral or discontinuous Galerkin (DG) method, or the locations of grid points when a finite difference method is used.
The particular choices of Basis and Quadrature determine where the collocation points of a Mesh are located in an Element. For a spectral or DG method, integrals using \(N\) collocation points with Gauss quadrature are exact to polynomial order \(p=2N1\). GaussLobatto quadrature is exact only to polynomial order \(p=2N3\), but includes collocation points at the domain boundary. For a finite difference method, one needs to choose the order of the scheme (and hence the weights, differentiation matrix, integration weights, and interpolant) locally in space and time to handle discontinuous solutions.
Gauss
or GaussLobatto
when using Basis::Legendre or Basis::Chebyshev.CellCentered
or FaceCentered
when using Basis::FiniteDifference.Gauss
for the first dimension and Equiangular
in the second dimension.Uninitialized
value. const std::pair< Matrix, Matrix > & Spectral::boundary_interpolation_matrices  (  const Mesh< 1 > &  mesh  ) 
Matrices that interpolate to the lower and upper boundaries of the element.
Assumes that the logical coordinates are \([1, 1]\). The first element of the pair interpolates to \(\xi=1\) and the second to \(\xi=1\). These are just the Lagrange interpolating polynomials evaluated at \(\xi=\pm1\). For GaussLobatto points the only nonzero element is at the boundaries and is one and so is not implemented.
const std::pair< Matrix, Matrix > & Spectral::boundary_interpolation_matrices  (  size_t  num_points  ) 
Matrices that interpolate to the lower and upper boundaries of the element.
Assumes that the logical coordinates are \([1, 1]\). The first element of the pair interpolates to \(\xi=1\) and the second to \(\xi=1\). These are just the Lagrange interpolating polynomials evaluated at \(\xi=\pm1\). For GaussLobatto points the only nonzero element is at the boundaries and is one and so is not implemented.
const std::pair< DataVector, DataVector > & Spectral::boundary_interpolation_term  (  const Mesh< 1 > &  mesh  ) 
Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhustype boundary conditions in a discontinuous Galerkin scheme using Gauss points.
Assumes that the logical coordinates are \([1, 1]\). The interpolation is done by assuming the time derivative correction is zero on interior nodes. With a nodal Lagrange polynomial basis this means that only the \(\ell^{\mathrm{GaussLobatto}}_{0}\) and \(\ell^{\mathrm{GaussLobatto}}_{N}\) polynomials/basis functions contribute to the correction. In order to interpolate the correction from the nodes at the boundary, the GaussLobatto Lagrange polynomials must be evaluated at the Gauss grid points. The returned pair of DataVector
s stores
\begin{align*} &\ell^{\mathrm{GaussLobatto}}_{0}(\xi_j^{\mathrm{Gauss}}), \\ &\ell^{\mathrm{GaussLobatto}}_{N}(\xi_j^{\mathrm{Gauss}}). \end{align*}
This is a different correction from lifting. Lifting is done using the mass matrix, which is an integral over the basis functions, while here we use interpolation.
const std::pair< DataVector, DataVector > & Spectral::boundary_interpolation_term  (  size_t  num_points  ) 
Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhustype boundary conditions in a discontinuous Galerkin scheme using Gauss points.
Assumes that the logical coordinates are \([1, 1]\). The interpolation is done by assuming the time derivative correction is zero on interior nodes. With a nodal Lagrange polynomial basis this means that only the \(\ell^{\mathrm{GaussLobatto}}_{0}\) and \(\ell^{\mathrm{GaussLobatto}}_{N}\) polynomials/basis functions contribute to the correction. In order to interpolate the correction from the nodes at the boundary, the GaussLobatto Lagrange polynomials must be evaluated at the Gauss grid points. The returned pair of DataVector
s stores
\begin{align*} &\ell^{\mathrm{GaussLobatto}}_{0}(\xi_j^{\mathrm{Gauss}}), \\ &\ell^{\mathrm{GaussLobatto}}_{N}(\xi_j^{\mathrm{Gauss}}). \end{align*}
This is a different correction from lifting. Lifting is done using the mass matrix, which is an integral over the basis functions, while here we use interpolation.
const std::pair< DataVector, DataVector > & Spectral::boundary_lifting_term  (  const Mesh< 1 > &  mesh  ) 
Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points.
Assumes that the logical coordinates are \([1, 1]\). The first element of the pair is the Lagrange polyonmials evaluated at \(\xi=1\) divided by the weights and the second at \(\xi=1\). Specifically,
\begin{align*} \frac{\ell_j(\xi=\pm1)}{w_j} \end{align*}
const std::pair< DataVector, DataVector > & Spectral::boundary_lifting_term  (  size_t  num_points  ) 
Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points.
Assumes that the logical coordinates are \([1, 1]\). The first element of the pair is the Lagrange polyonmials evaluated at \(\xi=1\) divided by the weights and the second at \(\xi=1\). Specifically,
\begin{align*} \frac{\ell_j(\xi=\pm1)}{w_j} \end{align*}
const DataVector & Spectral::collocation_points  (  const Mesh< 1 > &  mesh  ) 
Collocation points for a onedimensional mesh.
const DataVector & Spectral::collocation_points  (  size_t  num_points  ) 
Collocation points.
num_points  The number of collocation points 
std::pair< DataVector, DataVector > Spectral::compute_collocation_points_and_weights  (  size_t  num_points  ) 
Compute the collocation points and weights associated to the basis and quadrature.
This function is expected to return the tuple \((\xi_k,w_k)\) where the \(\xi_k\) are the collocation points and the \(w_k\) are defined in the description of quadrature_weights(size_t)
.
FiniteDifference
basis or CellCentered
and FaceCentered
quadratures, the weights are defined to integrate with the midpoint method Differentiation matrix for a onedimensional mesh.
const Matrix & Spectral::differentiation_matrix  (  size_t  num_points  ) 
Matrix used to compute the derivative of a function.
For a function represented by the nodal coefficients \(u_j\) a matrix multiplication with the differentiation matrix \(D_{ij}\) gives the coefficients of the function's derivative. Since \(u(x)\) is expanded in Lagrange polynomials \(u(x)=\sum_j u_j l_j(x)\) the differentiation matrix is computed as \(D_{ij}=l_j^\prime(\xi_i)\) where the \(\xi_i\) are the collocation points.
The finite difference matrix uses summation by parts operators, \(D_{21}, D_{42}, D_{43}\), and \(D_{65}\) from [53].
num_points  The number of collocation points 
Differentiation matrix for a onedimensional mesh.
const Matrix & Spectral::differentiation_matrix_transpose  (  size_t  num_points  ) 
Matrix used to compute the derivative of a function.
For a function represented by the nodal coefficients \(u_j\) a matrix multiplication with the differentiation matrix \(D_{ij}\) gives the coefficients of the function's derivative. Since \(u(x)\) is expanded in Lagrange polynomials \(u(x)=\sum_j u_j l_j(x)\) the differentiation matrix is computed as \(D_{ij}=l_j^\prime(\xi_i)\) where the \(\xi_i\) are the collocation points.
The finite difference matrix uses summation by parts operators, \(D_{21}, D_{42}, D_{43}\), and \(D_{65}\) from [53].
num_points  The number of collocation points 
size_t Spectral::hash  (  const std::array< Spectral::ChildSize, DimMinusOne > &  mortar_size  ) 
Performs a perfect hash of the mortars into \(2^{d1}\) slots on the range \([0, 2^{d1})\).
This is particularly useful when hashing into staticallysized maps based on the number of dimensions.
Indefinite integration matrix for a onedimensional mesh.
const Matrix & Spectral::integration_matrix  (  size_t  num_points  ) 
Matrix used to perform an indefinite integral of a function over the logical grid. The left boundary condition is such that the integral is 0 at \(\xi=1\)
Currently only Legendre and Chebyshev polynomials are implemented, but we provide a derivation for how to compute the indefinite integration matrix for general Jacobi polynomials.
The Legendre polynomials have the identity:
\begin{align*} P_n(x) = \frac{1}{2n+1}\frac{d}{dx}\left(P_{n+1}(x)  P_{n1}(x)\right) \end{align*}
The goal is to evaluate the integral of a function \(u\) expanded in terms of Legendre polynomials as
\begin{align*} u(x) = \sum_{i=0}^{N} c_i P_i(x) \end{align*}
We similarly expand the indefinite integral of \(u\) as
\begin{align*} \left.\int u(y) dy\right\rvert_{x}=&\sum_{i=0}^N \tilde{c}_i P_i(x) \\ =&\left.\int\sum_{i=1}^{N}\frac{c_i}{2i+1} \left(P_{i+1}(y)P_{i1}(y)\right)dy\right\rvert_{x} + \tilde{c}_0 P_0(x) \\ =&\sum_{i=1}^{N}\left(\frac{c_{i1}}{2i1}  \frac{c_{i+1}}{2i+3}\right) P_i(x) + \tilde{c}_0 P_0(x) \end{align*}
Thus we get that for \(i>0\)
\begin{align*} \tilde{c}_i=\frac{c_{i1}}{2i1}\frac{c_{i+1}}{2i+3} \end{align*}
and \(\tilde{c}_0\) is a constant of integration, which we choose such that the integral is 0 at the left boundary of the domain ( \(x=1\)). The condition for this is:
\begin{align*} \tilde{c}_0=\sum_{i=1}^{N}(1)^{i+1}\tilde{c}_i \end{align*}
The matrix returned by this function is the product of the tridiagonal matrix for the \(\tilde{c}_i\) and the matrix for the boundary condition.
A similar derivation leads to the relations:
\begin{align*} \tilde{c}_i=&\frac{c_{i1}c_{i+1}}{2i},&\mathrm{if}\;i>1 \\ \tilde{c}_1=&c_0  \frac{c_2}{2},&\mathrm{if}\;i=1 \\ \end{align*}
We again have:
\begin{align*} \tilde{c}_0=\sum_{i=1}^N(1)^{i+1}\tilde{c}_i \end{align*}
These are then used to define the indefinite integration matrix.
For general Jacobi polynomials \(P^{(\alpha,\beta)}_n(x)\) given by
\begin{align*} (1x)^\alpha(1+x)^\beta P^{(\alpha,\beta)}_n(x)=\frac{(1)^n}{2^n n!} \frac{d^n}{dx^n}\left[(1x)^{\alpha+n}(1+x)^{\beta+n}\right] \end{align*}
we have that
\begin{align*} P^{(\alpha,\beta)}_n(x)=\frac{d}{dx}\left( b^{(\alpha,\beta)}_{n1,n}P^{(\alpha,\beta)}_{n1}(x) + b^{(\alpha,\beta)}_{n,n}P^{(\alpha,\beta)}_n(x) + b^{(\alpha,\beta)}_{n+1,n}P^{(\alpha,\beta)}_{n+1}(x) \right) \end{align*}
where
\begin{align*} b^{(\alpha,\beta)}_{n1,n}=&\frac{1}{n+\alpha+\beta} a^{(\alpha,\beta)}_{n1,n} \\ b^{(\alpha,\beta)}_{n,n}=&\frac{2}{\alpha+\beta} a^{(\alpha,\beta)}_{n,n} \\ b^{(\alpha,\beta)}_{n+1,n}=&\frac{1}{n+1} a^{(\alpha,\beta)}_{n+1,n} \\ a^{(\alpha,\beta)}_{n1,n}=&\frac{2(n+\alpha)(n+\beta)} {(2n+\alpha+\beta+1)(2n+\alpha+\beta)} \\ a^{(\alpha,\beta)}_{n,n}=&\frac{\alpha^2\beta^2} {(2n+\alpha+\beta+2)(2n+\alpha+\beta)} \\ a^{(\alpha,\beta)}_{n1,n}=&\frac{2(n+1)(n+\alpha+\beta+1)} {(2n+\alpha+\beta+2)(2n+\alpha+\beta+1)} \end{align*}
Following the same derivation we get that
\begin{align*} \tilde{c}_i=c_{i+1}b^{(\alpha,\beta)}_{i,i+1} +c_i b^{(\alpha,\beta)}_{i,i} +c_{i1}b^{(\alpha,\beta)}_{i,i1} \end{align*}
and the boundary condition is
\begin{align*} \tilde{c}_0=\sum_{i=1}^N(1)^{i+1} \frac{\Gamma(i+\alpha+1)}{i!\Gamma(\alpha+1)} \tilde{c}_i \end{align*}
where \(\Gamma(x)\) is the Gamma function.
Matrix Spectral::interpolation_matrix  (  const Mesh< 1 > &  mesh, 
const T &  target_points  
) 
Interpolation matrix to the target_points
for a onedimensional mesh.
Matrix Spectral::interpolation_matrix  (  size_t  num_points, 
const T &  target_points  
) 
Matrix used to interpolate to the target_points
.
BasisType
(often \([1,1]\)), the resulting matrix performs polynomial extrapolation instead of interpolation. The extrapolation will be correct but may suffer from reduced accuracy, especially for higherorder polynomials (i.e., larger values of num_points
).num_points  The number of collocation points 
target_points  The points to interpolate to 
Linear filter matrix for a onedimensional mesh.
const Matrix & Spectral::linear_filter_matrix  (  size_t  num_points  ) 
Matrix used to linearize a function.
Filters out all except the lowest two modes by applying \(\mathcal{V}^{1}\cdot\mathrm{diag}(1,1,0,0,...)\cdot\mathcal{V}\) to the nodal coefficients, where \(\mathcal{V}\) is the Vandermonde matrix computed in modal_to_nodal_matrix(size_t)
.
num_points  The number of collocation points 
Transformation matrix from modal to nodal coefficients for a onedimensional mesh.
const Matrix & Spectral::modal_to_nodal_matrix  (  size_t  num_points  ) 
Matrix used to transform from the spectral coefficients (modes) of a function to its nodal coefficients. Also referred to as the Vandermonde matrix.
The Vandermonde matrix is computed as \(\mathcal{V}_{ij}=\Phi_j(\xi_i)\) where the \(\Phi_j(x)\) are the spectral basis functions used in the modal expansion \(u(x)=\sum_j \widetilde{u}_j\Phi_j(x)\), e.g. normalized Legendre polynomials, and the \(\xi_i\) are the collocation points used to construct the interpolating Lagrange polynomials in the nodal expansion \(u(x)=\sum_j u_j l_j(x)\). Then the Vandermonde matrix arises as \(u(\xi_i)=u_i=\sum_j \widetilde{u}_j\Phi_j(\xi_i)=\sum_j \mathcal{V}_{ij}\widetilde{u}_j\).
num_points  The number of collocation points 
Transformation matrix from nodal to modal coefficients for a onedimensional mesh.
const Matrix & Spectral::nodal_to_modal_matrix  (  size_t  num_points  ) 
Matrix used to transform from the nodal coefficients of a function to its spectral coefficients (modes). Also referred to as the inverse Vandermonde matrix.
This is the inverse to the Vandermonde matrix \(\mathcal{V}\) computed in modal_to_nodal_matrix(size_t). It can be computed analytically for Gauss quadrature by evaluating \(\sum_j\mathcal{V}^{1}_{ij}u_j=\widetilde{u}_i= \frac{(u,\Phi_i)}{\gamma_i}\) for a Lagrange basis function \(u(x)=l_k(x)\) to find \(\mathcal{V}^{1}_{ij}=\mathcal{V}_{ji}\frac{w_j}{\gamma_i}\) where the \(w_j\) are the Gauss quadrature weights and \(\gamma_i\) is the norm square of the spectral basis function \(\Phi_i\).
num_points  The number of collocation points 
const Matrix & Spectral::projection_matrix_child_to_parent  (  const Mesh< 1 > &  child_mesh, 
const Mesh< 1 > &  parent_mesh,  
ChildSize  size,  
bool  operand_is_massive = false 

) 
The projection matrix from a child mesh to its parent.
The projection matrices returned by this function (and by projection_matrix_parent_to_child()) define orthogonal projection operators between the spaces of functions on a parent mesh and its children. These projections are usually the correct way to transfer data between meshes in a meshrefinement hierarchy, as well as between an element face and its adjacent mortars.
These functions assume that the child_mesh
is at least as fine as the parent_mesh
, i.e. functions on the parent_mesh
can be represented exactly on the child_mesh
. In practice this means that functions can be projected to a mortar (the child_mesh
) from both adjacent element faces (the parent_mesh
) without losing accuracy. Similarly, functions in a meshrefinement hierarchy don't lose accuracy when an element is split (hrefined). For this reason, the projection_matrix_child_to_parent
is sometimes referred to as a "restriction operator" and the projection_matrix_parent_to_child
as a "prolongation operator".
true
for the parameter operand_is_massive
(default is false
). The restriction operator for this case is just the transpose of the prolongation operator, i.e. just an interpolation matrix transpose. Note that the "massive" residual already takes the difference in element size between parent and children into account by including a Jacobian in the volume element of the integral.\begin{equation*} T_{jk} = \frac{2 j + 1}{2} 2^j \sum_{n=0}^{jk} \binom{j}{k+n} \binom{(j + k + n  1)/2}{j} \frac{(k + n)!^2}{(2 k + n + 1)! n!} \end{equation*}
const Matrix & Spectral::projection_matrix_parent_to_child  (  const Mesh< 1 > &  parent_mesh, 
const Mesh< 1 > &  child_mesh,  
ChildSize  size  
) 
The projection matrix from a parent mesh to one of its children.
const DataVector & Spectral::quadrature_weights  (  const Mesh< 1 > &  mesh  ) 
Quadrature weights for a onedimensional mesh.
const DataVector & Spectral::quadrature_weights  (  size_t  num_points  ) 
Weights to compute definite integrals.
These are the coefficients to contract with the nodal function values \(f_k\) to approximate the definite integral \(I[f]=\int f(x)\mathrm{d}x\).
Note that the term quadrature also often refers to the quantity \(Q[f]=\int f(x)w(x)\mathrm{d}x\approx \sum_k f_k w_k\). Here, \(w(x)\) denotes the basisspecific weight function w.r.t. to which the basis functions \(\Phi_k\) are orthogonal, i.e \(\int\Phi_i(x)\Phi_j(x)w(x)=0\) for \(i\neq j\). The weights \(w_k\) approximate this inner product. To approximate the definite integral \(I[f]\) we must employ the coefficients \(\frac{w_k}{w(\xi_k)}\) instead, where the \(\xi_k\) are the collocation points. These are the coefficients this function returns. Only for a unit weight function \(w(x)=1\), i.e. a Legendre basis, is \(I[f]=Q[f]\) so this function returns the \(w_k\) identically.
For a FiniteDifference
basis or CellCentered
and FaceCentered
quadratures, the interpretation of the quadrature weights in term of an approximation to \(I(q)\) remains correct, but its explanation in terms of orthonormal basis is not, i.e. we set \(w_k\) to the grid spacing at each point, and the inverse weight \(\frac{1}{w(\xi_k)}=1\) to recover the midpoint method for definite integrals.
num_points  The number of collocation points 
Matrix used to compute the divergence of the flux in weak form.
const Matrix & Spectral::weak_flux_differentiation_matrix  (  size_t  num_points  ) 
Matrix used to compute the divergence of the flux in weak form.
This is the transpose of the differentiation matrix multiplied by quadrature weights that appear in DG integrals:
\begin{equation} \frac{D^T_{ij}} \frac{w_j}{w_i} \end{equation}
num_points  The number of collocation points 

constexpr 
The amount we leftshift the basis integers by.
We do this to be able to represent the combined Basis and Quadrature as one 8bit integer.

constexpr 
Maximum number of allowed collocation points.
We choose a limit of 24 FD grid points because for DGsubcell the number of points in an element is 2 * number_dg_points  1
for cell centered, and 2 * number_dg_points
for facecentered. Because there is no way of generically retrieving the maximum number of grid points for a nonFD basis, we need to hardcode both values here. If the number of grid points is increased for the nonFD bases, it should also be increased for the FD basis. Note that for good taskbased parallelization 24 grid points is already a fairly large number.

constexpr 
Minimum number of possible collocation points for a quadrature type.
Since GaussLobatto quadrature has points on the domain boundaries it must have at least two collocation points. Gauss quadrature can have only one collocation point.
For CellCentered
the minimum number of points is 1, while for FaceCentered
it is 2.