Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : /// \cond
       9             : class DataVector;
      10             : class Matrix;
      11             : /// \endcond
      12             : 
      13             : namespace Spectral {
      14             : 
      15             : /*!
      16             :  * \ingroup SpectralGroup
      17             :  *
      18             :  * \brief A collection of helper functions for using a Fourier series
      19             :  *
      20             :  * \details The general real-valued Fourier series is given by:
      21             :  * \f[
      22             :  * f(x) = a_0 + \sum_{n=1}^{\infty} a_n \cos(nx) + b_n \sin(nx)
      23             :  * \f]
      24             :  * where
      25             :  * \f{align*}
      26             :  * a_0 &= \frac{1}{2\pi} \int_0^{2\pi} f(x) dx \\
      27             :  * a_n &= \frac{1}{\pi} \int_0^{2\pi} f(x) \cos(nx) dx \\
      28             :  * b_n &= \frac{1}{\pi} \int_0^{2\pi} f(x) \sin(nx) dx
      29             :  * \f}
      30             :  *
      31             :  * If a function is discretized at \f$N\f$ collocation points, the modal
      32             :  * representation will have \f$N\f$ spectral coefficients consisting of
      33             :  * \f{align*}
      34             :  * a_n & \qquad \text{for } n = 0, \ldots, \frac{N}{2} \\
      35             :  * b_n & \qquad \text{for } n = 1, \ldots, \frac{N-1}{2}
      36             :  * \f}
      37             :  * Thus to represent all terms through the harmonic \f$M\f$ requires \f$N = 2M
      38             :  * +1\f$ collocation points.
      39             :  *
      40             :  * For more details about using Fourier series see e.g. \cite Boyd2001 and
      41             :  * \cite Fornberg1996.
      42             :  */
      43           1 : class Fourier {
      44             :  public:
      45             :   /*!
      46             :    * \brief Value of the basis function \f$\Phi_k(x)\f$
      47             :    *
      48             :    * \details We define the basis functions as
      49             :    * \f{align*}{
      50             :    *   \Phi_k(x)
      51             :    *   &=\begin{cases}
      52             :    *     \hfil \cos(kx) \hfil & \text{if } k \geq 0 \\
      53             :    *     \hfil \sin(-kx) \hfil & \text{if } k < 0
      54             :    *   \end{cases}
      55             :    * \f}
      56             :    */
      57             :   template <typename T>
      58           1 :   static T basis_function_value(int k, const T& x);
      59             : 
      60             :   /*!
      61             :    * \brief The normalization square \f$c_k\f$ of the basis function
      62             :    * \f$\Phi_k(x)\f$, i.e. the definite integral of its square.
      63             :    *
      64             :    * \details
      65             :    * In particular,
      66             :    * \f{align*}{
      67             :    *   c_k
      68             :    *   &=\begin{cases}
      69             :    *     \hfil 2 \pi \hfil & \text{if } k = 0 \\
      70             :    *     \hfil \pi \hfil & \text{otherwise}
      71             :    *   \end{cases}
      72             :    * \f}
      73             :    */
      74           1 :   static double basis_function_normalization_square(int k);
      75             : 
      76             :   /*!
      77             :    * \brief Collocation points \f$\{x_i\}\f$
      78             :    *
      79             :    * \details The collocation points on the interval \f$0 \leq x < 2 \pi\f$ are
      80             :    * given by
      81             :    * \f[
      82             :    * x_i = \frac{2 \pi i}{N}
      83             :    * \f]
      84             :    */
      85           1 :   static DataVector collocation_points(size_t num_points);
      86             : 
      87             :   /*!
      88             :    * \brief Quadrature weights \f$\{w_i\}\f$
      89             :    *
      90             :    * \details The quadrature weights are given by
      91             :    * \f[
      92             :    * w_i = \frac{2 \pi}{N}
      93             :    * \f]
      94             :    */
      95           1 :   static DataVector quadrature_weights(size_t num_points);
      96             : 
      97             :   /*!
      98             :    * \brief Storage index (offset) \f$i\f$ into a ModalVector (representing the
      99             :    * coefficients) of the given mode \f$k\f$
     100             :    *
     101             :    * \details The modal coefficients are stored in a ModalVector as
     102             :    * \f$\{u_0, u_1, u_{-1}, u_2, u_{-2}, \ldots, u_M, u_{-M}\}\f$,
     103             :    * where \f$u_{-M}\f$ is omitted if the number of coefficients is even.
     104             :    * Therefore the storage index \f$i\f$ for the mode \f$u_k\f$ is:
     105             :    * \f{align*}{
     106             :    *   i
     107             :    *   &=\begin{cases}
     108             :    *     \hfil 0 \hfil & \text{if } k = 0 \\
     109             :    *     \hfil 2k-1 \hfil & \text{if } k > 0 \\
     110             :    *     \hfil -2k \hfil & \text{if } k < 0
     111             :    *   \end{cases}
     112             :    * \f}
     113             :    */
     114           1 :   static size_t modal_storage_index(int k);
     115             : 
     116             :   /*!
     117             :    * \brief Mode \f$k\f$ corresponding to given storage index (offset) \f$i\f$
     118             :    * into a ModalVector (representing the coefficients)
     119             :    *
     120             :    * \details The modal coefficients are stored in a ModalVector as
     121             :    * \f$\{u_0, u_1, u_{-1}, u_2, u_{-2}, \ldots, u_M, u_{-M}\}\f$,
     122             :    * where \f$u_{-M}\f$ is omitted if the number of coefficients is even.
     123             :    * Therefore the storage index \f$i\f$ for the mode \f$u_k\f$ is:
     124             :    * \f{align*}{
     125             :    *   k
     126             :    *   &=\begin{cases}
     127             :    *     \hfil 0 \hfil & \text{if } i = 0 \\
     128             :    *     \hfil 1 + \frac{i}{2} \hfil & \text{if } i \text{ is odd} \\
     129             :    *     \hfil -\frac{i}{2} \hfil & \text{if } i \text{ is even}
     130             :    *   \end{cases}
     131             :    * \f}
     132             :    */
     133           1 :   static int mode_at_storage_index(size_t storage_index);
     134             : 
     135             :   /*!
     136             :    * \brief Matrix \f$D_{i,j}\f$ used to obtain the first derivative
     137             :    *
     138             :    * \details The differentiation matrix is given by:
     139             :    * \f{align*}{
     140             :    *   D_{i,j}
     141             :    *   &=\begin{cases}
     142             :    *     \hfil 0 \hfil & \text{if } i = j \\
     143             :    *     \hfil \frac{1}{2} (-1)^{i-j} \csc \left[0.5 (x_i - x_j)\right] \hfil &
     144             :    *     \text{if } i \neq j,\; N \text{ is odd} \\
     145             :    *     \hfil \frac{1}{2} (-1)^{i-j} \cot \left[0.5 (x_i - x_j)\right] \hfil &
     146             :    *     \text{if } i \neq j,\; N \text{ is even}
     147             :    *   \end{cases}
     148             :    * \f}
     149             :    */
     150           1 :   static Matrix differentiation_matrix(size_t num_points);
     151             : 
     152             :   /*!
     153             :    * \brief %Matrix used to interpolate to the \p target_points.
     154             :    *
     155             :    * \details Each row of the matrix is given by the interpolation weights for
     156             :    * interpolating to a particular target point \f$x\f$.  At each target point,
     157             :    * \f$x\f$, the interpolation weights are given by:
     158             :    * \f{align*}{
     159             :    *   C_j(x)
     160             :    *   &=\begin{cases}
     161             :    *     \hfil \frac{1}{N} \sin \left[0.5 N(x - x_j)\right]
     162             :    *     \csc \left[0.5 (x - x_j)\right] \hfil & \text{if } N \text{ is odd} \\
     163             :    *     \hfil \frac{1}{N} \sin \left[0.5 N (x - x_j)\right]
     164             :    *     \cot \left[0.5 (x - x_j)\right] \hfil & \text{if } N \text{ is even}
     165             :    *   \end{cases}
     166             :    * \f}
     167             :    * where \f$x_j\f$ are the collocation_points.
     168             :    */
     169             :   template <typename T>
     170           1 :   static Matrix interpolation_matrix(size_t num_points, const T& target_points);
     171             : };
     172             : }  // namespace Spectral