ProductMaps.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines the class templates ProductOf2Maps and ProductOf3Maps.
6 
7 #pragma once
8 
9 #include <array>
10 #include <boost/none.hpp>
11 #include <boost/optional.hpp>
12 #include <cstddef>
13 #include <functional>
14 #include <utility>
15 
19 #include "Utilities/TMPL.hpp"
20 #include "Utilities/TypeTraits.hpp"
21 
22 /// \cond
23 namespace PUP {
24 class er;
25 } // namespace PUP
26 /// \endcond
27 
28 namespace domain {
29 namespace CoordinateMaps {
30 
31 namespace product_detail {
32 template <typename T, size_t Size, typename Map1, typename Map2,
33  typename Function, size_t... Is, size_t... Js>
35  const std::array<T, Size>& coords, const Map1& map1, const Map2& map2,
36  const Function func, std::integer_sequence<size_t, Is...> /*meta*/,
37  std::integer_sequence<size_t, Js...> /*meta*/) noexcept {
38  using UnwrappedT = tt::remove_cvref_wrap_t<T>;
39  return {
40  {func(
42  {coords[Is]...}},
43  map1)[Is]...,
44  func(
46  {coords[Map1::dim + Js]...}},
47  map2)[Js]...}};
48 }
49 
50 template <size_t Size, typename Map1, typename Map2, typename Function,
51  size_t... Is, size_t... Js>
52 boost::optional<std::array<double, Size>> apply_inverse(
53  const std::array<double, Size>& coords, const Map1& map1, const Map2& map2,
54  const Function func, std::integer_sequence<size_t, Is...> /*meta*/,
55  std::integer_sequence<size_t, Js...> /*meta*/) noexcept {
56  auto map1_func =
57  func(std::array<double, sizeof...(Is)>{{coords[Is]...}}, map1);
58  auto map2_func = func(
59  std::array<double, sizeof...(Js)>{{coords[Map1::dim + Js]...}}, map2);
60  if (map1_func and map2_func) {
61  return {{{map1_func.get()[Is]..., map2_func.get()[Js]...}}};
62  } else {
63  return boost::none;
64  }
65 }
66 
67 template <typename T, size_t Size, typename Map1, typename Map2,
68  typename Function, size_t... Is, size_t... Js>
69 tnsr::Ij<tt::remove_cvref_wrap_t<T>, Size, Frame::NoFrame> apply_jac(
70  const std::array<T, Size>& source_coords, const Map1& map1,
71  const Map2& map2, const Function func,
73  std::integer_sequence<size_t, Js...> /*meta*/) noexcept {
74  using UnwrappedT = tt::remove_cvref_wrap_t<T>;
75  auto map1_jac = func(
77  {source_coords[Is]...}},
78  map1);
79  auto map2_jac = func(
81  {source_coords[Map1::dim + Js]...}},
82  map2);
83  tnsr::Ij<UnwrappedT, Size, Frame::NoFrame> jac{
84  make_with_value<UnwrappedT>(dereference_wrapper(source_coords[0]), 0.0)};
85  for (size_t i = 0; i < Map1::dim; ++i) {
86  for (size_t j = 0; j < Map1::dim; ++j) {
87  jac.get(i, j) = std::move(map1_jac.get(i, j));
88  }
89  }
90  for (size_t i = 0; i < Map2::dim; ++i) {
91  for (size_t j = 0; j < Map2::dim; ++j) {
92  jac.get(Map1::dim + i, Map1::dim + j) = std::move(map2_jac.get(i, j));
93  }
94  }
95  return jac;
96 }
97 } // namespace product_detail
98 
99 /// \ingroup CoordinateMapsGroup
100 /// \brief Product of two codimension=0 CoordinateMaps.
101 ///
102 /// \tparam Map1 the map for the first coordinate(s)
103 /// \tparam Map2 the map for the second coordinate(s)
104 template <typename Map1, typename Map2>
106  public:
107  static constexpr size_t dim = Map1::dim + Map2::dim;
108  using map_list = tmpl::list<Map1, Map2>;
109  static_assert(dim == 2 or dim == 3,
110  "Only 2D and 3D maps are supported by ProductOf2Maps");
111 
112  // Needed for Charm++ serialization
113  ProductOf2Maps() = default;
114 
115  ProductOf2Maps(Map1 map1, Map2 map2) noexcept;
116 
117  template <typename T>
118  std::array<tt::remove_cvref_wrap_t<T>, dim> operator()(
119  const std::array<T, dim>& source_coords) const noexcept;
120 
121  boost::optional<std::array<double, dim>> inverse(
122  const std::array<double, dim>& target_coords) const noexcept;
123 
124  template <typename T>
125  tnsr::Ij<tt::remove_cvref_wrap_t<T>, dim, Frame::NoFrame> inv_jacobian(
126  const std::array<T, dim>& source_coords) const noexcept;
127 
128  template <typename T>
129  tnsr::Ij<tt::remove_cvref_wrap_t<T>, dim, Frame::NoFrame> jacobian(
130  const std::array<T, dim>& source_coords) const noexcept;
131 
132  // clang-tidy: google-runtime-references
133  void pup(PUP::er& p); // NOLINT
134 
135  bool is_identity() const noexcept { return is_identity_; }
136 
137  private:
138  friend bool operator==(const ProductOf2Maps& lhs,
139  const ProductOf2Maps& rhs) noexcept {
140  return lhs.map1_ == rhs.map1_ and lhs.map2_ == rhs.map2_ and
141  lhs.is_identity_ == rhs.is_identity_;
142  }
143 
144  Map1 map1_;
145  Map2 map2_;
146  bool is_identity_ = false;
147 };
148 
149 template <typename Map1, typename Map2>
150 ProductOf2Maps<Map1, Map2>::ProductOf2Maps(Map1 map1, Map2 map2) noexcept
151  : map1_(std::move(map1)),
152  map2_(std::move(map2)),
153  is_identity_(map1_.is_identity() and map2_.is_identity()) {}
154 
155 template <typename Map1, typename Map2>
156 template <typename T>
157 std::array<tt::remove_cvref_wrap_t<T>, ProductOf2Maps<Map1, Map2>::dim>
159  const std::array<T, dim>& source_coords) const noexcept {
160  return product_detail::apply_call(
161  source_coords, map1_,
162  map2_, [](const auto& point,
163  const auto& map) noexcept { return map(point); },
166 }
167 
168 template <typename Map1, typename Map2>
169 boost::optional<std::array<double, ProductOf2Maps<Map1, Map2>::dim>>
171  const std::array<double, dim>& target_coords) const noexcept {
172  return product_detail::apply_inverse(
173  target_coords, map1_,
174  map2_, [](const auto& point,
175  const auto& map) noexcept { return map.inverse(point); },
178 }
179 
180 template <typename Map1, typename Map2>
181 template <typename T>
182 tnsr::Ij<tt::remove_cvref_wrap_t<T>, ProductOf2Maps<Map1, Map2>::dim,
185  const std::array<T, dim>& source_coords) const noexcept {
186  return product_detail::apply_jac(
187  source_coords, map1_,
188  map2_, [](const auto& point,
189  const auto& map) noexcept { return map.inv_jacobian(point); },
192 }
193 
194 template <typename Map1, typename Map2>
195 template <typename T>
196 tnsr::Ij<tt::remove_cvref_wrap_t<T>, ProductOf2Maps<Map1, Map2>::dim,
199  const std::array<T, dim>& source_coords) const noexcept {
200  return product_detail::apply_jac(
201  source_coords, map1_,
202  map2_, [](const auto& point,
203  const auto& map) noexcept { return map.jacobian(point); },
206 }
207 
208 template <typename Map1, typename Map2>
209 void ProductOf2Maps<Map1, Map2>::pup(PUP::er& p) {
210  p | map1_;
211  p | map2_;
212  p | is_identity_;
213 }
214 
215 template <typename Map1, typename Map2>
216 bool operator!=(const ProductOf2Maps<Map1, Map2>& lhs,
217  const ProductOf2Maps<Map1, Map2>& rhs) noexcept {
218  return not(lhs == rhs);
219 }
220 
221 /// \ingroup CoordinateMaps
222 /// \brief Product of three one-dimensional CoordinateMaps.
223 template <typename Map1, typename Map2, typename Map3>
225  public:
226  static constexpr size_t dim = Map1::dim + Map2::dim + Map3::dim;
227  using map_list = tmpl::list<Map1, Map2, Map3>;
228  static_assert(dim == 3, "Only 3D maps are implemented for ProductOf3Maps");
229 
230  // Needed for Charm++ serialization
231  ProductOf3Maps() = default;
232 
233  ProductOf3Maps(Map1 map1, Map2 map2, Map3 map3) noexcept;
234 
235  template <typename T>
236  std::array<tt::remove_cvref_wrap_t<T>, dim> operator()(
237  const std::array<T, dim>& source_coords) const noexcept;
238 
239  boost::optional<std::array<double, dim>> inverse(
240  const std::array<double, dim>& target_coords) const noexcept;
241 
242  template <typename T>
243  tnsr::Ij<tt::remove_cvref_wrap_t<T>, dim, Frame::NoFrame> inv_jacobian(
244  const std::array<T, dim>& source_coords) const noexcept;
245 
246  template <typename T>
247  tnsr::Ij<tt::remove_cvref_wrap_t<T>, dim, Frame::NoFrame> jacobian(
248  const std::array<T, dim>& source_coords) const noexcept;
249 
250  // clang-tidy: google-runtime-references
251  void pup(PUP::er& p) noexcept; // NOLINT
252 
253  bool is_identity() const noexcept { return is_identity_; }
254 
255  private:
256  friend bool operator==(const ProductOf3Maps& lhs,
257  const ProductOf3Maps& rhs) noexcept {
258  return lhs.map1_ == rhs.map1_ and lhs.map2_ == rhs.map2_ and
259  lhs.map3_ == rhs.map3_ and lhs.is_identity_ == rhs.is_identity_;
260  }
261 
262  Map1 map1_;
263  Map2 map2_;
264  Map3 map3_;
265  bool is_identity_ = false;
266 };
267 
268 template <typename Map1, typename Map2, typename Map3>
270  Map3 map3) noexcept
271  : map1_(std::move(map1)),
272  map2_(std::move(map2)),
273  map3_(std::move(map3)),
274  is_identity_(map1_.is_identity() and map2_.is_identity() and
275  map3_.is_identity()) {}
276 
277 template <typename Map1, typename Map2, typename Map3>
278 template <typename T>
279 std::array<tt::remove_cvref_wrap_t<T>, ProductOf3Maps<Map1, Map2, Map3>::dim>
281  const std::array<T, dim>& source_coords) const noexcept {
282  using UnwrappedT = tt::remove_cvref_wrap_t<T>;
284  {source_coords[0]}})[0],
286  {source_coords[1]}})[0],
288  {source_coords[2]}})[0]}};
289 }
290 
291 template <typename Map1, typename Map2, typename Map3>
292 boost::optional<std::array<double, ProductOf3Maps<Map1, Map2, Map3>::dim>>
294  const std::array<double, dim>& target_coords) const noexcept {
295  auto c1 = map1_.inverse(std::array<double, 1>{{target_coords[0]}});
296  auto c2 = map2_.inverse(std::array<double, 1>{{target_coords[1]}});
297  auto c3 = map3_.inverse(std::array<double, 1>{{target_coords[2]}});
298  if (c1 and c2 and c3) {
299  return {{{c1.get()[0], c2.get()[0], c3.get()[0]}}};
300  } else {
301  return boost::none;
302  }
303 }
304 
305 template <typename Map1, typename Map2, typename Map3>
306 template <typename T>
307 tnsr::Ij<tt::remove_cvref_wrap_t<T>, ProductOf3Maps<Map1, Map2, Map3>::dim,
310  const std::array<T, dim>& source_coords) const noexcept {
311  using UnwrappedT = tt::remove_cvref_wrap_t<T>;
312  tnsr::Ij<UnwrappedT, dim, Frame::NoFrame> inv_jacobian_matrix{
313  make_with_value<UnwrappedT>(dereference_wrapper(source_coords[0]), 0.0)};
314  get<0, 0>(inv_jacobian_matrix) = get<0, 0>(map1_.inv_jacobian(
316  {source_coords[0]}}));
317  get<1, 1>(inv_jacobian_matrix) = get<0, 0>(map2_.inv_jacobian(
319  {source_coords[1]}}));
320  get<2, 2>(inv_jacobian_matrix) = get<0, 0>(map3_.inv_jacobian(
322  {source_coords[2]}}));
323  return inv_jacobian_matrix;
324 }
325 template <typename Map1, typename Map2, typename Map3>
326 template <typename T>
327 tnsr::Ij<tt::remove_cvref_wrap_t<T>, ProductOf3Maps<Map1, Map2, Map3>::dim,
330  const std::array<T, dim>& source_coords) const noexcept {
331  using UnwrappedT = tt::remove_cvref_wrap_t<T>;
332  tnsr::Ij<UnwrappedT, dim, Frame::NoFrame> jacobian_matrix{
333  make_with_value<UnwrappedT>(dereference_wrapper(source_coords[0]), 0.0)};
334  get<0, 0>(jacobian_matrix) = get<0, 0>(
336  {source_coords[0]}}));
337  get<1, 1>(jacobian_matrix) = get<0, 0>(
339  {source_coords[1]}}));
340  get<2, 2>(jacobian_matrix) = get<0, 0>(
342  {source_coords[2]}}));
343  return jacobian_matrix;
344 }
345 template <typename Map1, typename Map2, typename Map3>
346 void ProductOf3Maps<Map1, Map2, Map3>::pup(PUP::er& p) noexcept {
347  p | map1_;
348  p | map2_;
349  p | map3_;
350  p | is_identity_;
351 }
352 
353 template <typename Map1, typename Map2, typename Map3>
354 bool operator!=(const ProductOf3Maps<Map1, Map2, Map3>& lhs,
355  const ProductOf3Maps<Map1, Map2, Map3>& rhs) noexcept {
356  return not(lhs == rhs);
357 }
358 } // namespace CoordinateMaps
359 } // namespace domain
Defines function dereference_wrapper.
Definition: Strahlkorper.hpp:14
Definition: BlockId.hpp:16
Product of three one-dimensional CoordinateMaps.
Definition: ProductMaps.hpp:224
Defines classes for Tensor.
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
Wraps the template metaprogramming library used (brigand)
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
Product of two codimension=0 CoordinateMaps.
Definition: ProductMaps.hpp:105
Defines type traits, some of which are future STL type_traits header.
Defines make_with_value.