ElementMap.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 <limits>
9 #include <memory>
10 
12 #include "Domain/ElementId.hpp"
13 #include "Utilities/Gsl.hpp"
14 #include "Utilities/MakeArray.hpp"
15 
16 namespace PUP {
17 class er;
18 } // namespace PUP
19 /// \cond
20 template <typename SourceFrame, typename TargetFrame, size_t Dim>
21 class CoordinateMapBase;
22 /// \endcond
23 
24 /*!
25  * \ingroup ComputationalDomainGroup
26  * \brief The CoordinateMap for the Element from the Logical frame to the
27  * `TargetFrame`
28  *
29  * An ElementMap takes a CoordinateMap for a Block and an ElementId as input,
30  * and then "prepends" the correct affine map to the CoordinateMap so that the
31  * map corresponds to the coordinate map for the Element rather than the Block.
32  * This allows DomainCreators to only specify the maps for the Blocks without
33  * worrying about how the domain may be decomposed beyond that.
34  */
35 template <size_t Dim, typename TargetFrame>
36 class ElementMap {
37  public:
38  static constexpr size_t dim = Dim;
40  using target_frame = TargetFrame;
41 
42  /// \cond HIDDEN_SYMBOLS
43  ElementMap() = default;
44  /// \endcond
45 
46  ElementMap(
47  ElementId<Dim> element_id,
49  block_map) noexcept;
50 
52  noexcept {
53  return *block_map_;
54  }
55 
56  const ElementId<Dim>& element_id() const noexcept { return element_id_; }
57 
58  template <typename T>
59  tnsr::I<T, Dim, TargetFrame> operator()(
60  tnsr::I<T, Dim, Frame::Logical> source_point) const noexcept {
61  apply_affine_transformation_to_point(source_point);
62  return block_map_->operator()(std::move(source_point));
63  }
64 
65  template <typename T>
66  tnsr::I<T, Dim, Frame::Logical> inverse(
67  tnsr::I<T, Dim, TargetFrame> target_point) const noexcept {
68  auto source_point{
69  block_map_->inverse(std::move(target_point)).get()};
70  // Apply the affine map to the points
71  for (size_t d = 0; d < Dim; ++d) {
72  source_point.get(d) =
73  source_point.get(d) * gsl::at(map_inverse_slope_, d) +
74  gsl::at(map_inverse_offset_, d);
75  }
76  return source_point;
77  }
78 
79  template <typename T>
80  InverseJacobian<T, Dim, Frame::Logical, TargetFrame> inv_jacobian(
81  tnsr::I<T, Dim, Frame::Logical> source_point) const noexcept {
82  apply_affine_transformation_to_point(source_point);
83  auto inv_jac = block_map_->inv_jacobian(std::move(source_point));
84  for (size_t d = 0; d < Dim; ++d) {
85  for (size_t i = 0; i < Dim; ++i) {
86  inv_jac.get(d, i) *= gsl::at(inverse_jacobian_, d);
87  }
88  }
89  return inv_jac;
90  }
91 
92  template <typename T>
93  Jacobian<T, Dim, Frame::Logical, TargetFrame> jacobian(
94  tnsr::I<T, Dim, Frame::Logical> source_point) const noexcept {
95  apply_affine_transformation_to_point(source_point);
96  auto jac = block_map_->jacobian(std::move(source_point));
97  for (size_t d = 0; d < Dim; ++d) {
98  for (size_t i = 0; i < Dim; ++i) {
99  jac.get(i, d) *= gsl::at(jacobian_, d);
100  }
101  }
102  return jac;
103  }
104 
105  // clang-tidy: do not use references
106  void pup(PUP::er& p) noexcept; // NOLINT
107 
108  private:
109  template <typename T>
110  void apply_affine_transformation_to_point(
111  tnsr::I<T, Dim, Frame::Logical>& source_point) const noexcept {
112  for (size_t d = 0; d < Dim; ++d) {
113  source_point.get(d) = source_point.get(d) * gsl::at(map_slope_, d) +
114  gsl::at(map_offset_, d);
115  }
116  }
117 
119  block_map_{nullptr};
120  ElementId<Dim> element_id_{};
121  // map_slope_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) -
122  // segment_ids[i].endpoint(Side::Lower))
123  std::array<double, Dim> map_slope_{
124  make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
125  // map_offset_[i] = 0.5 * (segment_ids[i].endpoint(Side::Upper) +
126  // segment_ids[i].endpoint(Side::Lower))
127  std::array<double, Dim> map_offset_{
128  make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
129  // map_inverse_slope_[i] = 1.0 / map_slope_[i]
130  std::array<double, Dim> map_inverse_slope_{
131  make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
132  // map_inverse_offset_[i] = -map_offset_[i] / map_slope_[i]
133  std::array<double, Dim> map_inverse_offset_{
134  make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
135  // Note: The Jacobian is diagonal
136  std::array<double, Dim> jacobian_{map_slope_};
137  std::array<double, Dim> inverse_jacobian_{map_inverse_slope_};
138 };
Definition: Strahlkorper.hpp:14
T signaling_NaN(T... args)
Defines function make_array.
Abstract base class for CoordinateMap.
Definition: CoordinateMap.hpp:49
Defines class ElementId.
The CoordinateMap for the Element from the Logical frame to the TargetFrame
Definition: ElementMap.hpp:36
Definition: IndexType.hpp:42
Defines classes for Tensor.
Defines functions and classes from the GSL.
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