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              : 
       9              : #include "DataStructures/ApplyMatrices.hpp"
      10              : #include "DataStructures/DataVector.hpp"
      11              : #include "DataStructures/Index.hpp"
      12              : #include "DataStructures/Matrix.hpp"
      13              : #include "DataStructures/Variables.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              : /// \endcond
      23              : 
      24              : namespace intrp {
      25              : 
      26              : /// \ingroup NumericalAlgorithmsGroup
      27              : /// \brief Interpolate data from a Mesh onto a regular grid of points.
      28              : ///
      29              : /// The target points must lie on a tensor-product grid that is aligned with the
      30              : /// source Mesh. In any direction where the source and target points share an
      31              : /// identical Mesh (i.e., where the underlying 1-dimensional meshes share the
      32              : /// same extent, basis, and quadrature), the code is optimized to avoid
      33              : /// performing identity interpolations.
      34              : ///
      35              : /// Note, however, that in each dimension of the target grid, the points can
      36              : /// be freely distributed; in particular, the grid points need not be the
      37              : /// collocation points corresponding to a particular basis and quadrature. Note
      38              : /// also that the target grid need not overlap the source grid. In this case,
      39              : /// polynomial extrapolation is performed, with order set by the order of the
      40              : /// basis in the source grid. The extrapolation will be correct but may suffer
      41              : /// from reduced accuracy, especially for higher-order polynomials.
      42              : template <size_t Dim>
      43            1 : class RegularGrid {
      44              :  public:
      45              :   /// \brief An interpolator between two regular grids.
      46              :   ///
      47              :   /// When the optional third argument is NOT passed, creates an interpolator
      48              :   /// between two regular meshes.
      49              :   ///
      50              :   /// The optional third argument allows the caller to override the distribution
      51              :   /// of grid points in any dimension(s) of the target grid. Each non-empty
      52              :   /// element of `override_target_mesh_with_1d_logical_coords` gives the logical
      53              :   /// coordinates which will override the default coordinates of `target_mesh`.
      54            1 :   RegularGrid(const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
      55              :               const std::array<DataVector, Dim>&
      56              :                   override_target_mesh_with_1d_logical_coords = {});
      57              : 
      58            0 :   RegularGrid();
      59              : 
      60              :   // NOLINTNEXTLINE(google-runtime-references)
      61            0 :   void pup(PUP::er& p);
      62              : 
      63              :   /// @{
      64              :   /// \brief Interpolate a Variables onto the target points.
      65              :   template <typename TagsList>
      66            1 :   void interpolate(gsl::not_null<Variables<TagsList>*> result,
      67              :                    const Variables<TagsList>& vars) const;
      68              :   template <typename TagsList>
      69            1 :   Variables<TagsList> interpolate(const Variables<TagsList>& vars) const;
      70              :   /// @}
      71              : 
      72              :   /// @{
      73              :   /// \brief Interpolate a DataVector onto the target points.
      74              :   ///
      75              :   /// \note When interpolating multiple tensors, the Variables interface is more
      76              :   /// efficient. However, this DataVector interface is useful for applications
      77              :   /// where only some components of a Tensor or Variables need to be
      78              :   /// interpolated.
      79            1 :   void interpolate(gsl::not_null<DataVector*> result,
      80              :                    const DataVector& input) const;
      81            1 :   DataVector interpolate(const DataVector& input) const;
      82            1 :   void interpolate(gsl::not_null<ComplexDataVector*> result,
      83              :                    const ComplexDataVector& input) const;
      84            1 :   ComplexDataVector interpolate(const ComplexDataVector& input) const;
      85              :   /// @}
      86              : 
      87              :   /// \brief Return the internally-stored matrices that interpolate from the
      88              :   /// source grid to the target grid.
      89              :   ///
      90              :   /// Each matrix interpolates in one logical direction, and can be empty if the
      91              :   /// interpolation in that direction is the identity (i.e., if the source
      92              :   /// and target grids have the same logical coordinates in this direction).
      93            1 :   const std::array<Matrix, Dim>& interpolation_matrices() const;
      94              : 
      95              :  private:
      96              :   template <size_t LocalDim>
      97              :   // NOLINTNEXTLINE(readability-redundant-declaration)
      98            0 :   friend bool operator==(const RegularGrid<LocalDim>& lhs,
      99              :                          const RegularGrid<LocalDim>& rhs);
     100              : 
     101            0 :   size_t number_of_target_points_{};
     102            0 :   Index<Dim> source_extents_;
     103            0 :   std::array<Matrix, Dim> interpolation_matrices_;
     104              : };
     105              : 
     106              : template <size_t Dim>
     107              : template <typename TagsList>
     108            0 : void RegularGrid<Dim>::interpolate(
     109              :     const gsl::not_null<Variables<TagsList>*> result,
     110              :     const Variables<TagsList>& vars) const {
     111              :   if (result->number_of_grid_points() != number_of_target_points_) {
     112              :     result->initialize(number_of_target_points_);
     113              :   }
     114              :   apply_matrices(result, interpolation_matrices_, vars, source_extents_);
     115              : }
     116              : 
     117              : template <size_t Dim>
     118              : template <typename TagsList>
     119              : Variables<TagsList> RegularGrid<Dim>::interpolate(
     120              :     const Variables<TagsList>& vars) const {
     121              :   Variables<TagsList> result;
     122              :   interpolate(make_not_null(&result), vars);
     123              :   return result;
     124              : }
     125              : 
     126              : template <size_t Dim>
     127            0 : bool operator!=(const RegularGrid<Dim>& lhs, const RegularGrid<Dim>& rhs);
     128              : 
     129              : }  // namespace intrp