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