Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <array>
       7              : #include <cstddef>
       8              : #include <optional>
       9              : 
      10              : #include "DataStructures/Matrix.hpp"
      11              : #include "DataStructures/Tensor/IndexType.hpp"
      12              : #include "DataStructures/Tensor/TypeAliases.hpp"
      13              : #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
      14              : #include "Utilities/Gsl.hpp"
      15              : 
      16              : /// \cond
      17              : class DataVector;
      18              : template <size_t Dim>
      19              : class Mesh;
      20              : namespace PUP {
      21              : class er;
      22              : }  // namespace PUP
      23              : /// \endcond
      24              : 
      25              : namespace intrp {
      26              : /*!
      27              :  * \brief Interpolates by doing partial summation in each dimension using
      28              :  * one-dimensional interpolation
      29              :  *
      30              :  * \details The one-dimensional matrices used to do the interpolation depend
      31              :  * upon the Spectral::Basis used in each dimension:
      32              :  * - For a Chebyshev or Legendre basis, the matrices are given by
      33              :  * Spectral::fornberg_interpolation_matrix at the quadrature points of the
      34              :  * source_mesh.  (These are equivalent to those returned by
      35              :  * Spectral::interpolation_matrix.)
      36              :  * - For a Fourier basis, the matrix is given by
      37              :  * Spectral::fourier_interpolation_matrix at the quadrature points of the
      38              :  * source_mesh
      39              :  *
      40              :  * For multidimensional bases such as SphericalHarmonic, ZernikeB2, or
      41              :  * ZernikeB3, the matrices used for interpolating cannot be applied per
      42              :  * dimension but must be handled specially.
      43              :  *
      44              :  */
      45              : template <size_t Dim>
      46            1 : class Cardinal {
      47              :  public:
      48            0 :   Cardinal(
      49              :       const Mesh<Dim>& source_mesh,
      50              :       const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points);
      51            0 :   Cardinal(const Mesh<Dim>& source_mesh,
      52              :            const tnsr::I<double, Dim, Frame::ElementLogical>& target_point);
      53              : 
      54            0 :   Cardinal();
      55              : 
      56              :   /// Interpolates the function `f` provided on the `source_mesh` to the
      57              :   /// `target_points` with which the interpolator was constructed.
      58            1 :   DataVector interpolate(const DataVector& f) const;
      59              : 
      60              :   /// The one-dimensional interpolation matrices used to do the interpolation
      61            1 :   const std::array<Matrix, Dim>& interpolation_matrices() const;
      62              : 
      63              :   // NOLINTNEXTLINE(google-runtime-references)
      64            0 :   void pup(PUP::er& p);
      65              : 
      66              :   /// Logic for `set_zernike_b2_weights()`. Also used by IrregularInterpolant.
      67            1 :   static void compute_zernike_b2_weights(
      68              :       gsl::not_null<Matrix*> weights, const Mesh<Dim>& source_mesh,
      69              :       const std::array<Matrix, Dim>& interpolation_matrices,
      70              :       size_t n_target_points);
      71              : 
      72              :   /// Logic for `set_zernike_b3_weights()`. Also used by IrregularInterpolant.
      73            1 :   static void compute_zernike_b3_weights(
      74              :       gsl::not_null<Matrix*> weights, const Mesh<Dim>& source_mesh,
      75              :       const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points,
      76              :       const ylm::Spherepack& b3_ylm, size_t n_target_points);
      77              : 
      78              :  private:
      79              :   /// Precomputes `zernike_b2_weights_`, which is all work independent of
      80              :   /// `f_source`, to avoid redundant computations. This is only needed when
      81              :   /// the source mesh has a B2 basis
      82            1 :   void set_zernike_b2_weights();
      83              : 
      84              :   /// General routine called by `interpolate()` for a mesh using B2 bases
      85            1 :   DataVector interpolate_zernike_b2(const DataVector& f_source) const;
      86              : 
      87              :   /// Precomputes `zernike_b3_weights_`, which is all work independent of
      88              :   /// `f_source`, to avoid redundant computations. This is only needed when
      89              :   /// the source mesh has a B3 basis
      90            1 :   void set_zernike_b3_weights();
      91              : 
      92              :   /// General routine called by `interpolate()` for a mesh using B3 bases
      93            1 :   DataVector interpolate_zernike_b3(const DataVector& f_source) const;
      94              : 
      95              :   template <size_t LocalDim>
      96              :   // NOLINTNEXTILNE(readability-redundant-declaration)
      97            0 :   friend bool operator==(const Cardinal<LocalDim>& lhs,
      98              :                          const Cardinal<LocalDim>& rhs);
      99              : 
     100            0 :   size_t n_target_points_ = 0;
     101              :   // Only required by B3 weights
     102              :   std::optional<tnsr::I<DataVector, Dim, Frame::ElementLogical>>
     103            0 :       target_points_{};
     104            0 :   Mesh<Dim> source_mesh_{};
     105            0 :   std::array<Matrix, Dim> interpolation_matrices_{};
     106            0 :   bool using_spherical_harmonics_{false};
     107            0 :   bool using_zernike_b2_{false};
     108            0 :   Matrix zernike_b2_weights_{};
     109            0 :   bool using_zernike_b3_{false};
     110            0 :   Matrix zernike_b3_weights_{};
     111            0 :   std::optional<ylm::Spherepack> b3_ylm_{};
     112              : };
     113              : 
     114              : template <size_t Dim>
     115            0 : bool operator==(const Cardinal<Dim>& lhs, const Cardinal<Dim>& rhs);
     116              : 
     117              : template <size_t Dim>
     118            0 : bool operator!=(const Cardinal<Dim>& lhs, const Cardinal<Dim>& rhs);
     119              : }  // namespace intrp