IrregularInterpolant.hpp
1 // 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"
12 #include "ErrorHandling/Assert.hpp"
13 #include "Utilities/Blas.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 /// Interpolates a `Variables` onto an arbitrary set of points.
29 template <size_t Dim>
30 class Irregular {
31  public:
32  Irregular(
33  const Mesh<Dim>& source_mesh,
34  const tnsr::I<DataVector, Dim, Frame::Logical>& target_points) noexcept;
35  Irregular();
36 
37  // clang-tidy: no runtime references
38  /// Serialization for Charm++
39  void pup(PUP::er& p) noexcept; // NOLINT
40 
41  //@{
42  /// Performs the interpolation on a `Variables` with grid points corresponding
43  /// to the `Mesh<Dim>` specified in the constructor.
44  /// The result is a `Variables` whose internal `DataVector` goes over the
45  /// list of target_points that were specified in the constructor.
46  /// \note for the void function, `result` will be resized to the proper size.
47  template <typename TagsList>
48  void interpolate(gsl::not_null<Variables<TagsList>*> result,
49  const Variables<TagsList>& vars) const noexcept;
50  template <typename TagsList>
51  Variables<TagsList> interpolate(const Variables<TagsList>& vars) const
52  noexcept;
53  //@}
54 
55  private:
56  friend bool operator==(const Irregular& lhs, const Irregular& rhs) noexcept {
57  return lhs.interpolation_matrix_ == rhs.interpolation_matrix_;
58  }
59  Matrix interpolation_matrix_;
60 };
61 
62 template <size_t Dim>
63 template <typename TagsList>
65  const gsl::not_null<Variables<TagsList>*> result,
66  const Variables<TagsList>& vars) const noexcept {
67  // For matrix multiplication of Interp . Source = Result:
68  // matrix Interp is m rows by k columns
69  // matrix Source is k rows by n columns
70  // matrix Result is m rows by n columns
71  const size_t m = interpolation_matrix_.rows();
72  const size_t k = interpolation_matrix_.columns();
73  const size_t n = vars.number_of_independent_components;
74  ASSERT(k == vars.number_of_grid_points(),
75  "Number of grid points in source 'vars', "
76  << vars.number_of_grid_points()
77  << ",\n disagrees with the size of the source_mesh, " << k
78  << ", that was passed into the constructor");
79  if (result->number_of_grid_points() != m) {
80  *result = Variables<TagsList>(m, 0.);
81  }
82  dgemm_('n', 'n', m, n, k, 1.0, interpolation_matrix_.data(), m, vars.data(),
83  k, 0.0, result->data(), m);
84 }
85 
86 template <size_t Dim>
87 template <typename TagsList>
88 Variables<TagsList> Irregular<Dim>::interpolate(
89  const Variables<TagsList>& vars) const noexcept {
90  Variables<TagsList> result;
91  interpolate(make_not_null(&result), vars);
92  return result;
93 }
94 
95 template <size_t Dim>
96 bool operator!=(const Irregular<Dim>& lhs,
97  const Irregular<Dim>& rhs) noexcept;
98 
99 } // namespace intrp
100 
Defines class Matrix.
Definition: Strahlkorper.hpp:14
Definition: AddTemporalIdsToInterpolationTarget.hpp:17
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Interpolates a Variables onto an arbitrary set of points.
Definition: IrregularInterpolant.hpp:30
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
Defines class Variables.
Declares the interfaces for the BLAS used.
Defines a list of useful type aliases for tensors.
void dgemm_(const char &TRANSA, const char &TRANSB, const size_t &M, const size_t &N, const size_t &K, const double &ALPHA, const double *A, const size_t &LDA, const double *B, const size_t &LDB, const double &BETA, double *C, const size_t &LDC)
Perform a matrix-matrix multiplication.
Definition: Blas.hpp:89
Defines macro ASSERT.
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