ApplyMatrices.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 #include <ostream>
9 
10 #include "DataStructures/DataVector.hpp"
12 #include "ErrorHandling/Assert.hpp"
14 #include "Utilities/Gsl.hpp"
15 
16 /// \cond
17 template <size_t Dim>
18 class Index;
19 // IWYU pragma: no_forward_declare Variables
20 /// \endcond
21 
22 namespace apply_matrices_detail {
23 template <size_t Dim, bool... DimensionIsIdentity>
24 struct Impl {
25  template <typename MatrixType>
26  static void apply(gsl::not_null<double*> result,
27  const std::array<MatrixType, Dim>& matrices,
28  const double* data, const Index<Dim>& extents,
29  size_t number_of_independent_components) noexcept;
30 };
31 
32 template <typename MatrixType, size_t Dim>
33 size_t result_size(const std::array<MatrixType, Dim>& matrices,
34  const Index<Dim>& extents) noexcept {
35  size_t num_points_result = 1;
36  for (size_t d = 0; d < Dim; ++d) {
37  const size_t cols = dereference_wrapper(gsl::at(matrices, d)).columns();
38  if (cols == 0) {
39  // An empty matrix is treated as the identity.
40  num_points_result *= extents[d];
41  } else {
42  ASSERT(cols == extents[d],
43  "Matrix " << d << " has wrong number of columns: "
44  << cols << " (expected " << extents[d] << ")");
45  num_points_result *= dereference_wrapper(gsl::at(matrices, d)).rows();
46  }
47  }
48  return num_points_result;
49 }
50 } // namespace apply_matrices_detail
51 
52 /// \ingroup NumericalAlgorithmsGroup
53 /// \brief Multiply by matrices in each dimension
54 ///
55 /// Multiplies each stripe in the first dimension of `u` by
56 /// `matrices[0]`, each stripe in the second dimension of `u` by
57 /// `matrices[1]`, and so on. If any of the matrices are empty they
58 /// will be treated as the identity, but the matrix multiplications
59 /// will be skipped for increased efficiency.
60 //@{
61 template <typename VariableTags, typename MatrixType, size_t Dim>
62 void apply_matrices(const gsl::not_null<Variables<VariableTags>*> result,
63  const std::array<MatrixType, Dim>& matrices,
64  const Variables<VariableTags>& u,
65  const Index<Dim>& extents) noexcept {
66  ASSERT(u.number_of_grid_points() == extents.product(),
67  "Mismatch between extents (" << extents.product()
68  << ") and variables (" << u.number_of_grid_points() << ").");
69  ASSERT(result->number_of_grid_points() ==
70  apply_matrices_detail::result_size(matrices, extents),
71  "result has wrong size. Expected "
72  << apply_matrices_detail::result_size(matrices, extents)
73  << ", received " << result->number_of_grid_points());
74  apply_matrices_detail::Impl<Dim>::apply(result->data(), matrices, u.data(),
75  extents,
76  u.number_of_independent_components);
77 }
78 
79 template <typename VariableTags, typename MatrixType, size_t Dim>
80 Variables<VariableTags> apply_matrices(
81  const std::array<MatrixType, Dim>& matrices,
82  const Variables<VariableTags>& u, const Index<Dim>& extents) noexcept {
83  Variables<VariableTags> result(
84  apply_matrices_detail::result_size(matrices, extents));
85  apply_matrices(make_not_null(&result), matrices, u, extents);
86  return result;
87 }
88 
89 template <typename MatrixType, size_t Dim>
91  const std::array<MatrixType, Dim>& matrices,
92  const DataVector& u, const Index<Dim>& extents) noexcept {
93  ASSERT(u.size() == extents.product(),
94  "Mismatch between extents (" << extents.product() << ") and size ("
95  << u.size() << ").");
96  ASSERT(result->size() ==
97  apply_matrices_detail::result_size(matrices, extents),
98  "result has wrong size. Expected "
99  << apply_matrices_detail::result_size(matrices, extents)
100  << ", received " << result->size());
101  apply_matrices_detail::Impl<Dim>::apply(result->data(), matrices, u.data(),
102  extents, 1);
103 }
104 
105 template <typename MatrixType, size_t Dim>
107  const DataVector& u,
108  const Index<Dim>& extents) noexcept {
109  DataVector result(apply_matrices_detail::result_size(matrices, extents));
110  apply_matrices(make_not_null(&result), matrices, u, extents);
111  return result;
112 }
113 //@}
Defines function dereference_wrapper.
Definition: ApplyMatrices.cpp:91
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
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
Defines class Variables.
Defines macro ASSERT.
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
decltype(auto) dereference_wrapper(T &&t)
Returns the reference object held by a reference wrapper, if a non-reference_wrapper type is passed i...
Definition: DereferenceWrapper.hpp:22
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid...
Definition: Gsl.hpp:124