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 <memory>
8 : #include <string>
9 : #include <type_traits>
10 : #include <unordered_map>
11 :
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/Tensor/Identity.hpp"
14 : #include "DataStructures/Tensor/Tensor.hpp"
15 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
16 : #include "Utilities/DereferenceWrapper.hpp"
17 : #include "Utilities/ErrorHandling/FloatingPointExceptions.hpp"
18 : #include "Utilities/Gsl.hpp"
19 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
20 :
21 : namespace domain {
22 : namespace CoordinateMap_detail {
23 : /// Returns a value that aliases `x` where that is cheaper than copying: a
24 : /// non-owning view for `DataVector` and a copy for arithmetic types such as
25 : /// `double`. Coordinate maps that pass (a subset of) their source coordinates
26 : /// on to another map use this function to avoid copying the underlying data.
27 : /// Coordinate maps take their source coordinates by const reference, so the
28 : /// data is never modified through the view. The view must not outlive `x`.
29 : template <typename T>
30 : T view_or_copy(const T& x) {
31 : if constexpr (std::is_same_v<T, DataVector>) {
32 : return {const_cast<double*>(x.data()), x.size()}; // NOLINT
33 : } else {
34 : return x;
35 : }
36 : }
37 : /// @{
38 : /// Call the map passing in the time and FunctionsOfTime if the map is
39 : /// time-dependent
40 : template <typename T, size_t Dim, typename Map>
41 : void apply_map(
42 : const gsl::not_null<std::array<T, Dim>*> t_map_point, const Map& the_map,
43 : const double /*t*/,
44 : const std::unordered_map<
45 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
46 : /*functions_of_time*/,
47 : const std::false_type /*is_time_independent*/) {
48 : if (LIKELY(not the_map.is_identity())) {
49 : *t_map_point = the_map(*t_map_point);
50 : }
51 : }
52 :
53 : template <typename T, size_t Dim, typename Map>
54 : void apply_map(
55 : const gsl::not_null<std::array<T, Dim>*> t_map_point, const Map& the_map,
56 : const double t,
57 : const std::unordered_map<
58 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
59 : functions_of_time,
60 : const std::true_type
61 : /*is_time_dependent*/) {
62 : ASSERT(not functions_of_time.empty(),
63 : "A function of time must be present if the maps are time-dependent.");
64 : ASSERT(
65 : [t]() {
66 : const ScopedFpeState disable_fpes(false);
67 : return not std::isnan(t);
68 : }(),
69 : "The time must not be NaN for time-dependent maps.");
70 : *t_map_point = the_map(*t_map_point, t, functions_of_time);
71 : }
72 :
73 : template <typename T, size_t Dim, typename Map>
74 : auto apply_map(
75 : const Map& the_map, const std::array<T, Dim>& source_points,
76 : const double /*t*/,
77 : const std::unordered_map<
78 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
79 : /*functions_of_time*/,
80 : const std::false_type /*is_time_independent*/) {
81 : if (LIKELY(not the_map.is_identity())) {
82 : return the_map(source_points);
83 : }
84 : std::decay_t<decltype(the_map(source_points))> result{};
85 : for (size_t i = 0; i < result.size(); ++i) {
86 : gsl::at(result, i) = gsl::at(source_points, i);
87 : }
88 : return result;
89 : }
90 :
91 : template <typename T, size_t Dim, typename Map>
92 : auto apply_map(
93 : const Map& the_map, const std::array<T, Dim>& source_points, const double t,
94 : const std::unordered_map<
95 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
96 : functions_of_time,
97 : const std::true_type
98 : /*is_time_dependent*/) {
99 : // Note: We don't forward to the return-by-not-null version to avoid
100 : // additional allocations of the target points array. That is, we would
101 : // allocate the target points array once here, and then again inside the call
102 : // to the coordinate map.
103 : ASSERT(not functions_of_time.empty(),
104 : "A function of time must be present if the maps are time-dependent.");
105 : ASSERT(
106 : [t]() {
107 : const ScopedFpeState disable_fpes(false);
108 : return not std::isnan(t);
109 : }(),
110 : "The time must not be NaN for time-dependent maps.");
111 : return the_map(source_points, t, functions_of_time);
112 : }
113 : /// @}
114 :
115 : /// @{
116 : template <typename T, size_t Dim, typename Map>
117 : auto apply_inverse_map(
118 : const Map& the_map, const std::array<T, Dim>& target_points,
119 : const double /*t*/,
120 : const std::unordered_map<
121 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
122 : /*functions_of_time*/,
123 : const std::false_type /*is_time_independent*/) {
124 : if (LIKELY(not the_map.is_identity())) {
125 : return the_map.inverse(target_points);
126 : }
127 : using UnwrappedT = tt::remove_cvref_wrap_t<T>;
128 : std::decay_t<decltype(the_map.inverse(target_points))> result{
129 : std::array<UnwrappedT, Dim>{}};
130 : for (size_t i = 0; i < target_points.size(); ++i) {
131 : gsl::at(*result, i) = gsl::at(target_points, i);
132 : }
133 : return result;
134 : }
135 :
136 : template <typename T, size_t Dim, typename Map>
137 : auto apply_inverse_map(
138 : const Map& the_map, const std::array<T, Dim>& target_points, const double t,
139 : const std::unordered_map<
140 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
141 : functions_of_time,
142 : const std::true_type
143 : /*is_time_dependent*/) {
144 : ASSERT(not functions_of_time.empty(),
145 : "A function of time must be present if the maps are time-dependent.");
146 : ASSERT(
147 : [t]() {
148 : const ScopedFpeState disable_fpes(false);
149 : return not std::isnan(t);
150 : }(),
151 : "The time must not be NaN for time-dependent maps.");
152 : return the_map.inverse(target_points, t, functions_of_time);
153 : }
154 : /// @}
155 :
156 : /// @{
157 : /// Compute the frame velocity
158 : template <typename T, size_t Dim, typename Map>
159 : auto apply_frame_velocity(
160 : const Map& /*the_map*/, const std::array<T, Dim>& source_points,
161 : const double /*t*/,
162 : const std::unordered_map<
163 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
164 : /*functions_of_time*/,
165 : const std::false_type /*is_time_independent*/) {
166 : return make_array<Map::dim, tt::remove_cvref_wrap_t<T>>(
167 : make_with_value<tt::remove_cvref_wrap_t<T>>(
168 : dereference_wrapper(source_points[0]), 0.0));
169 : }
170 :
171 : template <typename T, size_t Dim, typename Map>
172 : auto apply_frame_velocity(
173 : const Map& the_map, const std::array<T, Dim>& source_points, const double t,
174 : const std::unordered_map<
175 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
176 : functions_of_time,
177 : const std::true_type
178 : /*is_time_dependent*/) {
179 : ASSERT(not functions_of_time.empty(),
180 : "A function of time must be present if the maps are time-dependent.");
181 : ASSERT(
182 : [t]() {
183 : const ScopedFpeState disable_fpes(false);
184 : return not std::isnan(t);
185 : }(),
186 : "The time must not be NaN for time-dependent maps.");
187 : return the_map.frame_velocity(source_points, t, functions_of_time);
188 : }
189 : /// @}
190 :
191 : /// @{
192 : /// Compute the Jacobian
193 : template <typename T, size_t Dim, typename Map>
194 : auto apply_jacobian(
195 : const Map& the_map, const std::array<T, Dim>& source_points,
196 : const double /*t*/,
197 : const std::unordered_map<
198 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
199 : /*functions_of_time*/,
200 : const std::false_type /*is_time_independent*/) {
201 : if (LIKELY(not the_map.is_identity())) {
202 : return the_map.jacobian(source_points);
203 : }
204 : return identity<Dim>(dereference_wrapper(source_points[0]));
205 : }
206 :
207 : template <typename T, size_t Dim, typename Map>
208 : auto apply_jacobian(
209 : const Map& the_map, const std::array<T, Dim>& source_points, const double t,
210 : const std::unordered_map<
211 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
212 : functions_of_time,
213 : const std::true_type
214 : /*is_time_dependent*/) {
215 : ASSERT(not functions_of_time.empty(),
216 : "A function of time must be present if the maps are time-dependent.");
217 : ASSERT(
218 : [t]() {
219 : const ScopedFpeState disable_fpes(false);
220 : return not std::isnan(t);
221 : }(),
222 : "The time must not be NaN for time-dependent maps.");
223 : if (LIKELY(not the_map.is_identity())) {
224 : return the_map.jacobian(source_points, t, functions_of_time);
225 : }
226 : return identity<Dim>(dereference_wrapper(source_points[0]));
227 : }
228 : /// @}
229 :
230 : /// @{
231 : /// Compute the Jacobian
232 : template <typename T, size_t Dim, typename Map>
233 : auto apply_inverse_jacobian(
234 : const Map& the_map, const std::array<T, Dim>& source_points,
235 : const double /*t*/,
236 : const std::unordered_map<
237 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
238 : /*functions_of_time*/,
239 : const std::false_type /*is_time_independent*/) {
240 : if (LIKELY(not the_map.is_identity())) {
241 : return the_map.inv_jacobian(source_points);
242 : }
243 : return identity<Dim>(dereference_wrapper(source_points[0]));
244 : }
245 :
246 : template <typename T, size_t Dim, typename Map>
247 : auto apply_inverse_jacobian(
248 : const Map& the_map, const std::array<T, Dim>& source_points, const double t,
249 : const std::unordered_map<
250 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
251 : functions_of_time,
252 : const std::true_type
253 : /*is_time_dependent*/) {
254 : ASSERT(not functions_of_time.empty(),
255 : "A function of time must be present if the maps are time-dependent.");
256 : ASSERT(
257 : [t]() {
258 : const ScopedFpeState disable_fpes(false);
259 : return not std::isnan(t);
260 : }(),
261 : "The time must not be NaN for time-dependent maps.");
262 : if (LIKELY(not the_map.is_identity())) {
263 : return the_map.inv_jacobian(source_points, t, functions_of_time);
264 : }
265 : return identity<Dim>(dereference_wrapper(source_points[0]));
266 : }
267 : /// @}
268 : } // namespace CoordinateMap_detail
269 : } // namespace domain
|
