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 : /// \brief Interpolate multiple variables on the grid to the target points. 72 1 : void interpolate(gsl::not_null<gsl::span<double>*> result, 73 : const gsl::span<const double>& input) const; 74 : 75 : private: 76 0 : friend bool operator==(const Irregular& lhs, const Irregular& rhs) { 77 : return lhs.interpolation_matrix_ == rhs.interpolation_matrix_; 78 : } 79 0 : Matrix interpolation_matrix_; 80 : }; 81 : 82 : template <size_t Dim> 83 : template <typename TagsList> 84 0 : void Irregular<Dim>::interpolate( 85 : const gsl::not_null<Variables<TagsList>*> result, 86 : const Variables<TagsList>& vars) const { 87 : if (UNLIKELY(result->number_of_grid_points() != 88 : interpolation_matrix_.rows())) { 89 : *result = Variables<TagsList>(interpolation_matrix_.rows(), 0.); 90 : } 91 : ASSERT(interpolation_matrix_.columns() == vars.number_of_grid_points(), 92 : "Number of grid points in source 'vars', " 93 : << vars.number_of_grid_points() 94 : << ",\n disagrees with the size of the source_mesh, " 95 : << interpolation_matrix_.columns() 96 : << ", that was passed into the constructor"); 97 : gsl::span<double> result_span{result->data(), result->size()}; 98 : const gsl::span<const double> vars_span{vars.data(), vars.size()}; 99 : interpolate(make_not_null(&result_span), vars_span); 100 : } 101 : 102 : template <size_t Dim> 103 : template <typename TagsList> 104 : Variables<TagsList> Irregular<Dim>::interpolate( 105 : const Variables<TagsList>& vars) const { 106 : Variables<TagsList> result{interpolation_matrix_.rows()}; 107 : interpolate(make_not_null(&result), vars); 108 : return result; 109 : } 110 : 111 : template <size_t Dim> 112 0 : bool operator!=(const Irregular<Dim>& lhs, const Irregular<Dim>& rhs); 113 : 114 : } // namespace intrp