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