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