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           0 :   Irregular(
      43             :       const Mesh<Dim>& source_mesh,
      44             :       const tnsr::I<double, Dim, Frame::ElementLogical>& target_point);
      45           0 :   Irregular();
      46             : 
      47             :   /// Serialization for Charm++
      48             :   // NOLINTNEXTLINE(google-runtime-references)
      49           1 :   void pup(PUP::er& p);
      50             : 
      51             :   /// @{
      52             :   /// Performs the interpolation on a `Variables` with grid points corresponding
      53             :   /// to the `Mesh<Dim>` specified in the constructor.
      54             :   /// The result is a `Variables` whose internal `DataVector` goes over the
      55             :   /// list of target_points that were specified in the constructor.
      56             :   /// \note for the void function, `result` will be resized to the proper size.
      57             :   template <typename TagsList>
      58           1 :   void interpolate(gsl::not_null<Variables<TagsList>*> result,
      59             :                    const Variables<TagsList>& vars) const;
      60             :   template <typename TagsList>
      61           1 :   Variables<TagsList> interpolate(const Variables<TagsList>& vars) const;
      62             :   /// @}
      63             : 
      64             :   /// @{
      65             :   /// \brief Interpolate a DataVector or ComplexDataVector onto the target
      66             :   /// points.
      67             :   ///
      68             :   /// \note When interpolating multiple tensors, the Variables interface is more
      69             :   /// efficient. However, this DataVector interface is useful for applications
      70             :   /// where only some components of a Tensor or Variables need to be
      71             :   /// interpolated.
      72           1 :   void interpolate(gsl::not_null<DataVector*> result,
      73             :                    const DataVector& input) const;
      74           1 :   DataVector interpolate(const DataVector& input) const;
      75           1 :   void interpolate(gsl::not_null<ComplexDataVector*> result,
      76             :                    const ComplexDataVector& input) const;
      77           1 :   ComplexDataVector interpolate(const ComplexDataVector& input) const;
      78             :   /// @}
      79             : 
      80             :   /// \brief Interpolate multiple variables on the grid to the target points.
      81             :   /// @{
      82           1 :   void interpolate(gsl::not_null<gsl::span<double>*> result,
      83             :                    const gsl::span<const double>& input) const;
      84           1 :   void interpolate(gsl::not_null<gsl::span<std::complex<double>>*> result,
      85             :                    const gsl::span<const std::complex<double>>& input) const;
      86           1 :   void interpolate(gsl::not_null<gsl::span<float>*> result,
      87             :                    const gsl::span<const float>& input) const;
      88             :   /// @}
      89             : 
      90             :  private:
      91           0 :   friend bool operator==(const Irregular& lhs, const Irregular& rhs) {
      92             :     return lhs.interpolation_matrix_ == rhs.interpolation_matrix_;
      93             :   }
      94           0 :   Matrix interpolation_matrix_;
      95             : };
      96             : 
      97             : template <size_t Dim>
      98             : template <typename TagsList>
      99           0 : void Irregular<Dim>::interpolate(
     100             :     const gsl::not_null<Variables<TagsList>*> result,
     101             :     const Variables<TagsList>& vars) const {
     102             :   if (UNLIKELY(result->number_of_grid_points() !=
     103             :                interpolation_matrix_.rows())) {
     104             :     *result = Variables<TagsList>(interpolation_matrix_.rows(), 0.);
     105             :   }
     106             :   ASSERT(interpolation_matrix_.columns() == vars.number_of_grid_points(),
     107             :          "Number of grid points in source 'vars', "
     108             :              << vars.number_of_grid_points()
     109             :              << ",\n disagrees with the size of the source_mesh, "
     110             :              << interpolation_matrix_.columns()
     111             :              << ", that was passed into the constructor");
     112             :   auto result_span = gsl::make_span(result->data(), result->size());
     113             :   const auto vars_span = gsl::make_span(vars.data(), vars.size());
     114             :   interpolate(make_not_null(&result_span), vars_span);
     115             : }
     116             : 
     117             : template <size_t Dim>
     118             : template <typename TagsList>
     119             : Variables<TagsList> Irregular<Dim>::interpolate(
     120             :     const Variables<TagsList>& vars) const {
     121             :   Variables<TagsList> result{interpolation_matrix_.rows()};
     122             :   interpolate(make_not_null(&result), vars);
     123             :   return result;
     124             : }
     125             : 
     126             : template <size_t Dim>
     127           0 : bool operator!=(const Irregular<Dim>& lhs, const Irregular<Dim>& rhs);
     128             : 
     129             : }  // namespace intrp