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