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