CoordinateMapHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <memory>
8 #include <string>
9 #include <type_traits>
10 #include <unordered_map>
11 
12 #include "DataStructures/Tensor/Identity.hpp"
14 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
17 #include "Utilities/Gsl.hpp"
18 #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
19 
20 namespace domain {
21 namespace CoordinateMap_detail {
22 // @{
23 /// Call the map passing in the time and FunctionsOfTime if the map is
24 /// time-dependent
25 template <typename T, size_t Dim, typename Map>
26 void apply_map(
27  const gsl::not_null<std::array<T, Dim>*> t_map_point, const Map& the_map,
28  const double /*t*/,
29  const std::unordered_map<
31  /*functions_of_time*/,
32  const std::false_type /*is_time_independent*/) {
33  if (LIKELY(not the_map.is_identity())) {
34  *t_map_point = the_map(*t_map_point);
35  }
36 }
37 
38 template <typename T, size_t Dim, typename Map>
39 void apply_map(
40  const gsl::not_null<std::array<T, Dim>*> t_map_point, const Map& the_map,
41  const double t,
42  const std::unordered_map<
44  functions_of_time,
45  const std::true_type
46  /*is_time_dependent*/) {
47  ASSERT(not functions_of_time.empty(),
48  "A function of time must be present if the maps are time-dependent.");
49  ASSERT(
50  [t]() noexcept {
52  const bool isnan = std::isnan(t);
54  return not isnan;
55  }(),
56  "The time must not be NaN for time-dependent maps.");
57  *t_map_point = the_map(*t_map_point, t, functions_of_time);
58 }
59 
60 template <typename T, size_t Dim, typename Map>
61 auto apply_map(
62  const Map& the_map, const std::array<T, Dim>& source_points,
63  const double /*t*/,
64  const std::unordered_map<
66  /*functions_of_time*/,
67  const std::false_type /*is_time_independent*/) {
68  if (LIKELY(not the_map.is_identity())) {
69  return the_map(source_points);
70  }
71  std::decay_t<decltype(the_map(source_points))> result{};
72  for (size_t i = 0; i < result.size(); ++i) {
73  gsl::at(result, i) = gsl::at(source_points, i);
74  }
75  return result;
76 }
77 
78 template <typename T, size_t Dim, typename Map>
79 auto apply_map(
80  const Map& the_map, const std::array<T, Dim>& source_points, const double t,
81  const std::unordered_map<
83  functions_of_time,
84  const std::true_type
85  /*is_time_dependent*/) {
86  // Note: We don't forward to the return-by-not-null version to avoid
87  // additional allocations of the target points array. That is, we would
88  // allocate the target points array once here, and then again inside the call
89  // to the coordinate map.
90  ASSERT(not functions_of_time.empty(),
91  "A function of time must be present if the maps are time-dependent.");
92  ASSERT(
93  [t]() noexcept {
95  const bool isnan = std::isnan(t);
97  return not isnan;
98  }(),
99  "The time must not be NaN for time-dependent maps.");
100  return the_map(source_points, t, functions_of_time);
101 }
102 // @}
103 
104 // @{
105 template <typename T, size_t Dim, typename Map>
106 auto apply_inverse_map(
107  const Map& the_map, const std::array<T, Dim>& target_points,
108  const double /*t*/,
109  const std::unordered_map<
111  /*functions_of_time*/,
112  const std::false_type /*is_time_independent*/) {
113  if (LIKELY(not the_map.is_identity())) {
114  return the_map.inverse(target_points);
115  }
116  using UnwrappedT = tt::remove_cvref_wrap_t<T>;
117  std::decay_t<decltype(the_map.inverse(target_points))> result{
119  for (size_t i = 0; i < target_points.size(); ++i) {
120  gsl::at(*result, i) = gsl::at(target_points, i);
121  }
122  return result;
123 }
124 
125 template <typename T, size_t Dim, typename Map>
126 auto apply_inverse_map(
127  const Map& the_map, const std::array<T, Dim>& target_points, const double t,
128  const std::unordered_map<
130  functions_of_time,
131  const std::true_type
132  /*is_time_dependent*/) {
133  ASSERT(not functions_of_time.empty(),
134  "A function of time must be present if the maps are time-dependent.");
135  ASSERT(
136  [t]() noexcept {
138  const bool isnan = std::isnan(t);
140  return not isnan;
141  }(),
142  "The time must not be NaN for time-dependent maps.");
143  return the_map.inverse(target_points, t, functions_of_time);
144 }
145 // @}
146 
147 // @{
148 /// Compute the frame velocity
149 template <typename T, size_t Dim, typename Map>
150 auto apply_frame_velocity(
151  const Map& /*the_map*/, const std::array<T, Dim>& source_points,
152  const double /*t*/,
153  const std::unordered_map<
155  /*functions_of_time*/,
156  const std::false_type /*is_time_independent*/) {
157  return make_array<Map::dim, tt::remove_cvref_wrap_t<T>>(
158  make_with_value<tt::remove_cvref_wrap_t<T>>(
159  dereference_wrapper(source_points[0]), 0.0));
160 }
161 
162 template <typename T, size_t Dim, typename Map>
163 auto apply_frame_velocity(
164  const Map& the_map, const std::array<T, Dim>& source_points, const double t,
165  const std::unordered_map<
167  functions_of_time,
168  const std::true_type
169  /*is_time_dependent*/) {
170  ASSERT(not functions_of_time.empty(),
171  "A function of time must be present if the maps are time-dependent.");
172  ASSERT(
173  [t]() noexcept {
175  const bool isnan = std::isnan(t);
177  return not isnan;
178  }(),
179  "The time must not be NaN for time-dependent maps.");
180  return the_map.frame_velocity(source_points, t, functions_of_time);
181 }
182 // @}
183 
184 // @{
185 /// Compute the Jacobian
186 template <typename T, size_t Dim, typename Map>
187 auto apply_jacobian(
188  const Map& the_map, const std::array<T, Dim>& source_points,
189  const double /*t*/,
190  const std::unordered_map<
192  /*functions_of_time*/,
193  const std::false_type /*is_time_independent*/) {
194  if (LIKELY(not the_map.is_identity())) {
195  return the_map.jacobian(source_points);
196  }
197  return identity<Dim>(dereference_wrapper(source_points[0]));
198 }
199 
200 template <typename T, size_t Dim, typename Map>
201 auto apply_jacobian(
202  const Map& the_map, const std::array<T, Dim>& source_points, const double t,
203  const std::unordered_map<
205  functions_of_time,
206  const std::true_type
207  /*is_time_dependent*/) {
208  ASSERT(not functions_of_time.empty(),
209  "A function of time must be present if the maps are time-dependent.");
210  ASSERT(
211  [t]() noexcept {
213  const bool isnan = std::isnan(t);
215  return not isnan;
216  }(),
217  "The time must not be NaN for time-dependent maps.");
218  if (LIKELY(not the_map.is_identity())) {
219  return the_map.jacobian(source_points, t, functions_of_time);
220  }
221  return identity<Dim>(dereference_wrapper(source_points[0]));
222 }
223 // @}
224 
225 // @{
226 /// Compute the Jacobian
227 template <typename T, size_t Dim, typename Map>
228 auto apply_inverse_jacobian(
229  const Map& the_map, const std::array<T, Dim>& source_points,
230  const double /*t*/,
231  const std::unordered_map<
233  /*functions_of_time*/,
234  const std::false_type /*is_time_independent*/) {
235  if (LIKELY(not the_map.is_identity())) {
236  return the_map.inv_jacobian(source_points);
237  }
238  return identity<Dim>(dereference_wrapper(source_points[0]));
239 }
240 
241 template <typename T, size_t Dim, typename Map>
242 auto apply_inverse_jacobian(
243  const Map& the_map, const std::array<T, Dim>& source_points, const double t,
244  const std::unordered_map<
246  functions_of_time,
247  const std::true_type
248  /*is_time_dependent*/) {
249  ASSERT(not functions_of_time.empty(),
250  "A function of time must be present if the maps are time-dependent.");
251  ASSERT(
252  [t]() noexcept {
254  const bool isnan = std::isnan(t);
256  return not isnan;
257  }(),
258  "The time must not be NaN for time-dependent maps.");
259  if (LIKELY(not the_map.is_identity())) {
260  return the_map.inv_jacobian(source_points, t, functions_of_time);
261  }
262  return identity<Dim>(dereference_wrapper(source_points[0]));
263 }
264 // @}
265 } // namespace CoordinateMap_detail
266 } // namespace domain
FloatingPointExceptions.hpp
gsl::at
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:125
std::false_type
std::string
disable_floating_point_exceptions
void disable_floating_point_exceptions()
Definition: FloatingPointExceptions.cpp:37
enable_floating_point_exceptions
void enable_floating_point_exceptions()
Definition: FloatingPointExceptions.cpp:27
std::isnan
T isnan(T... args)
array
memory
std::decay_t
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
DereferenceWrapper.hpp
Gsl.hpp
Tensor.hpp
LIKELY
#define LIKELY(x)
Definition: Gsl.hpp:67
dereference_wrapper
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
std::unique_ptr
unordered_map
type_traits
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
string