RegularGridInterpolant.hpp
1 // 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/DataVector.hpp"
10 #include "DataStructures/Index.hpp"
11 #include "DataStructures/Matrix.hpp" // IWYU pragma: keep
13 #include "NumericalAlgorithms/LinearOperators/ApplyMatrices.hpp"
14 #include "Utilities/Gsl.hpp"
15 #include "Utilities/MakeArray.hpp"
16 
17 // IWYU pragma: no_forward_declare Variables
18 
19 /// \cond
20 template <size_t Dim>
21 class Mesh;
22 namespace PUP {
23 class er;
24 } // namespace PUP
25 /// \endcond
26 
27 namespace intrp {
28 
29 /// \ingroup NumericalAlgorithmsGroup
30 /// \brief Interpolate a Variables from a Mesh onto a regular grid of points.
31 ///
32 /// The target points must lie on a tensor-product grid that is aligned with the
33 /// source Mesh. In any direction where the source and target points share an
34 /// identical Mesh (i.e., where the underlying 1-dimensional meshes share the
35 /// same extent, basis, and quadrature), the code is optimized to avoid
36 /// performing identity interpolations.
37 ///
38 /// Note, however, that in each dimension of the target grid, the points can
39 /// be freely distributed; in particular, the grid points need not be the
40 /// collocation points corresponding to a particular basis and quadrature. Note
41 /// also that the target grid need not overlap the source grid. In this case,
42 /// polynomial extrapolation is performed, with order set by the order of the
43 /// basis in the source grid. The extrapolation will be correct but may suffer
44 /// from reduced accuracy, especially for higher-order polynomials.
45 template <size_t Dim>
46 class RegularGrid {
47  public:
48  /// \brief An interpolator between two regular grids.
49  ///
50  /// When the optional third argument is NOT passed, creates an interpolator
51  /// between two regular meshes.
52  ///
53  /// The optional third argument allows the caller to override the distribution
54  /// of grid points in any dimension(s) of the target grid. Each non-empty
55  /// element of `override_target_mesh_with_1d_logical_coords` gives the logical
56  /// coordinates which will override the default coordinates of `target_mesh`.
57  RegularGrid(const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
59  override_target_mesh_with_1d_logical_coords =
60  make_array<Dim>(DataVector())) noexcept;
61 
62  RegularGrid();
63 
64  // clang-tidy: no runtime references
65  void pup(PUP::er& p) noexcept; // NOLINT
66 
67  /// \brief Interpolate Variables onto new mesh.
68  //@{
69  template <typename TagsList>
70  void interpolate(gsl::not_null<Variables<TagsList>*> result,
71  const Variables<TagsList>& vars) const noexcept;
72  template <typename TagsList>
73  Variables<TagsList> interpolate(const Variables<TagsList>& vars) const
74  noexcept;
75  //@}
76 
77  private:
78  template <size_t LocalDim>
79  // NOLINTNEXTLINE(readability-redundant-declaration)
80  friend bool operator==(const RegularGrid<LocalDim>& lhs,
81  const RegularGrid<LocalDim>& rhs) noexcept;
82 
83  size_t number_of_target_points_{};
84  Index<Dim> source_extents_;
85  std::array<Matrix, Dim> interpolation_matrices_;
86 };
87 
88 template <size_t Dim>
89 template <typename TagsList>
91  const gsl::not_null<Variables<TagsList>*> result,
92  const Variables<TagsList>& vars) const noexcept {
93  if (result->number_of_grid_points() != number_of_target_points_) {
94  result->initialize(number_of_target_points_);
95  }
96  apply_matrices(result, interpolation_matrices_, vars, source_extents_);
97 }
98 
99 template <size_t Dim>
100 template <typename TagsList>
101 Variables<TagsList> RegularGrid<Dim>::interpolate(
102  const Variables<TagsList>& vars) const noexcept {
103  Variables<TagsList> result;
104  interpolate(make_not_null(&result), vars);
105  return result;
106 }
107 
108 template <size_t Dim>
109 bool operator!=(const RegularGrid<Dim>& lhs,
110  const RegularGrid<Dim>& rhs) noexcept;
111 
112 } // namespace intrp
Defines class Matrix.
Definition: Strahlkorper.hpp:14
Definition: AddTemporalIdsToInterpolationTarget.hpp:17
Defines function make_array.
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
void apply_matrices(const gsl::not_null< Variables< VariableTags > *> result, const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension.
Definition: ApplyMatrices.hpp:62
Defines class template Index.
Defines class Variables.
Interpolate a Variables from a Mesh onto a regular grid of points.
Definition: RegularGridInterpolant.hpp:46
Stores a collection of function values.
Definition: DataVector.hpp:46
An integer multi-index.
Definition: Index.hpp:28
Defines functions and classes from the GSL.
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12