Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <complex>
       7              : #include <cstddef>
       8              : 
       9              : #include "DataStructures/ComplexDataVector.hpp"
      10              : #include "DataStructures/DataVector.hpp"
      11              : #include "DataStructures/Matrix.hpp"
      12              : #include "DataStructures/Tensor/TypeAliases.hpp"
      13              : #include "DataStructures/Variables.hpp"
      14              : #include "Utilities/Blas.hpp"
      15              : #include "Utilities/ErrorHandling/Assert.hpp"
      16              : #include "Utilities/Gsl.hpp"
      17              : 
      18              : /// \cond
      19              : template <size_t Dim>
      20              : class Mesh;
      21              : namespace PUP {
      22              : class er;
      23              : }  // namespace PUP
      24              : 
      25              : /// \endcond
      26              : 
      27              : namespace intrp {
      28              : 
      29              : /// \ingroup NumericalAlgorithmsGroup
      30              : /// \brief Interpolates a `Variables` onto an arbitrary set of points.
      31              : ///
      32              : /// \details If the `source_mesh` uses Spectral::Basis::FiniteDifference,
      33              : /// linear interpolation is done in each dimension; otherwise it uses the
      34              : /// barycentric interpolation provided by Spectral::interpolation_matrix in each
      35              : /// dimension.
      36              : template <size_t Dim>
      37            1 : class Irregular {
      38              :  public:
      39            0 :   Irregular(
      40              :       const Mesh<Dim>& source_mesh,
      41              :       const tnsr::I<DataVector, Dim, Frame::ElementLogical>& target_points,
      42              :       const std::optional<size_t>& fd_to_fd_interp_order = std::nullopt);
      43            0 :   Irregular(const Mesh<Dim>& source_mesh,
      44              :             const tnsr::I<double, Dim, Frame::ElementLogical>& target_point,
      45              :             const std::optional<size_t>& fd_to_fd_interp_order = std::nullopt);
      46            0 :   Irregular();
      47              : 
      48              :   /// Serialization for Charm++
      49              :   // NOLINTNEXTLINE(google-runtime-references)
      50            1 :   void pup(PUP::er& p);
      51              : 
      52              :   /// @{
      53              :   /// Performs the interpolation on a `Variables` with grid points corresponding
      54              :   /// to the `Mesh<Dim>` specified in the constructor.
      55              :   /// The result is a `Variables` whose internal `DataVector` goes over the
      56              :   /// list of target_points that were specified in the constructor.
      57              :   /// \note for the void function, `result` will be resized to the proper size.
      58              :   template <typename TagsList>
      59            1 :   void interpolate(gsl::not_null<Variables<TagsList>*> result,
      60              :                    const Variables<TagsList>& vars) const;
      61              :   template <typename TagsList>
      62            1 :   Variables<TagsList> interpolate(const Variables<TagsList>& vars) const;
      63              :   /// @}
      64              : 
      65              :   /// @{
      66              :   /// \brief Interpolate a DataVector or ComplexDataVector onto the target
      67              :   /// points.
      68              :   ///
      69              :   /// \note When interpolating multiple tensors, the Variables interface is more
      70              :   /// efficient. However, this DataVector interface is useful for applications
      71              :   /// where only some components of a Tensor or Variables need to be
      72              :   /// interpolated.
      73            1 :   void interpolate(gsl::not_null<DataVector*> result,
      74              :                    const DataVector& input) const;
      75            1 :   DataVector interpolate(const DataVector& input) const;
      76            1 :   void interpolate(gsl::not_null<ComplexDataVector*> result,
      77              :                    const ComplexDataVector& input) const;
      78            1 :   ComplexDataVector interpolate(const ComplexDataVector& input) const;
      79              :   /// @}
      80              : 
      81              :   /// \brief Interpolate multiple variables on the grid to the target points.
      82              :   /// @{
      83            1 :   void interpolate(gsl::not_null<gsl::span<double>*> result,
      84              :                    const gsl::span<const double>& input) const;
      85            1 :   void interpolate(gsl::not_null<gsl::span<std::complex<double>>*> result,
      86              :                    const gsl::span<const std::complex<double>>& input) const;
      87            1 :   void interpolate(gsl::not_null<gsl::span<float>*> result,
      88              :                    const gsl::span<const float>& input) const;
      89              :   /// @}
      90              : 
      91              :  private:
      92            0 :   friend bool operator==(const Irregular& lhs, const Irregular& rhs) {
      93              :     return lhs.interpolation_matrix_ == rhs.interpolation_matrix_;
      94              :   }
      95            0 :   Matrix interpolation_matrix_;
      96              : };
      97              : 
      98              : template <size_t Dim>
      99              : template <typename TagsList>
     100            0 : void Irregular<Dim>::interpolate(
     101              :     const gsl::not_null<Variables<TagsList>*> result,
     102              :     const Variables<TagsList>& vars) const {
     103              :   if (UNLIKELY(result->number_of_grid_points() !=
     104              :                interpolation_matrix_.rows())) {
     105              :     *result = Variables<TagsList>(interpolation_matrix_.rows(), 0.);
     106              :   }
     107              :   ASSERT(interpolation_matrix_.columns() == vars.number_of_grid_points(),
     108              :          "Number of grid points in source 'vars', "
     109              :              << vars.number_of_grid_points()
     110              :              << ",\n disagrees with the size of the source_mesh, "
     111              :              << interpolation_matrix_.columns()
     112              :              << ", that was passed into the constructor");
     113              :   auto result_span = gsl::make_span(result->data(), result->size());
     114              :   const auto vars_span = gsl::make_span(vars.data(), vars.size());
     115              :   interpolate(make_not_null(&result_span), vars_span);
     116              : }
     117              : 
     118              : template <size_t Dim>
     119              : template <typename TagsList>
     120              : Variables<TagsList> Irregular<Dim>::interpolate(
     121              :     const Variables<TagsList>& vars) const {
     122              :   Variables<TagsList> result{interpolation_matrix_.rows()};
     123              :   interpolate(make_not_null(&result), vars);
     124              :   return result;
     125              : }
     126              : 
     127              : template <size_t Dim>
     128            0 : bool operator!=(const Irregular<Dim>& lhs, const Irregular<Dim>& rhs);
     129              : 
     130              : }  // namespace intrp