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 <cstddef>
8 : #include <iosfwd>
9 :
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Domain/Structure/Direction.hpp"
12 : #include "Domain/Structure/SegmentId.hpp"
13 : #include "Domain/Structure/Side.hpp"
14 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
15 : #include "Utilities/Gsl.hpp"
16 :
17 : namespace PUP {
18 : class er;
19 : } // namespace PUP
20 :
21 : /*!
22 : * \ingroup ComputationalDomainGroup
23 : * \brief A mapping of the logical coordinate axes of a host to the logical
24 : * coordinate axes of a neighbor of the host.
25 : *
26 : * Given a `size_t dimension`, a `Direction`, a `SegmentId`, or a `Mesh` of the
27 : * host, an `OrientationMap` will give the corresponding value in the neighbor.
28 : *
29 : * \tparam VolumeDim the dimension of the blocks.
30 : *
31 : * See the [tutorial](@ref tutorial_orientations) for information on how
32 : * OrientationMaps are used and constructed.
33 : *
34 : * \note If there is no discrete rotation between logical coordinates (e.g. the
35 : * angular coordinates of a spherical shell abutting a wedge of a cubed sphere)
36 : * specify Direction<VolumeDim>::self() as the mapped direction
37 : *
38 : */
39 : template <size_t VolumeDim>
40 1 : class OrientationMap {
41 : public:
42 0 : static constexpr uint16_t aligned_mask = 0b1000000000000000;
43 0 : static constexpr uint16_t version_mask = 0b0111000000000000;
44 :
45 : /// \brief Creates an OrientationMap in an uninitialized state.
46 : ///
47 : /// This can be helpful for debugging code. If you would like the identity
48 : /// map, please use `create_aligned()`.
49 1 : OrientationMap();
50 : /// Mapped directions relative to the positive (`Side::Upper`) direction in
51 : /// each logical direction.
52 1 : explicit OrientationMap(
53 : std::array<Direction<VolumeDim>, VolumeDim> mapped_directions);
54 0 : OrientationMap(
55 : const std::array<Direction<VolumeDim>, VolumeDim>& directions_in_host,
56 : const std::array<Direction<VolumeDim>, VolumeDim>&
57 : directions_in_neighbor);
58 0 : ~OrientationMap() = default;
59 0 : OrientationMap(const OrientationMap&) = default;
60 0 : OrientationMap& operator=(const OrientationMap&) = default;
61 0 : OrientationMap(OrientationMap&& /*rhs*/) = default;
62 0 : OrientationMap& operator=(OrientationMap&& /*rhs*/) = default;
63 :
64 : /// Creates an OrientationMap that is the identity map on directions.
65 : /// `is_aligned()` is `true` in this case.
66 1 : static OrientationMap<VolumeDim> create_aligned();
67 :
68 : /// True when mapped(Direction) == Direction
69 1 : bool is_aligned() const {
70 : ASSERT(bit_field_ != static_cast<uint16_t>(0b1 << 15),
71 : "Cannot use a default-constructed OrientationMap");
72 : return (bit_field_ bitand aligned_mask) == aligned_mask;
73 : }
74 :
75 : /// The corresponding dimension in the neighbor.
76 1 : size_t operator()(const size_t dim) const {
77 : ASSERT(bit_field_ != static_cast<uint16_t>(0b1 << 15),
78 : "Cannot use a default-constructed OrientationMap");
79 : const auto neighbor_direction = get_direction(dim);
80 : ASSERT(neighbor_direction.side() != Side::Self,
81 : "There is no corresponding dimension");
82 : return neighbor_direction.dimension();
83 : }
84 :
85 : /// The corresponding direction in the neighbor.
86 1 : Direction<VolumeDim> operator()(const Direction<VolumeDim>& direction) const {
87 : ASSERT(bit_field_ != static_cast<uint16_t>(0b1 << 15),
88 : "Cannot use a default-constructed OrientationMap");
89 : return direction.side() == Side::Upper
90 : ? get_direction(direction.dimension())
91 : : get_direction(direction.dimension()).opposite();
92 : }
93 :
94 : /// The corresponding SegmentIds in the neighbor.
95 1 : std::array<SegmentId, VolumeDim> operator()(
96 : const std::array<SegmentId, VolumeDim>& segmentIds) const;
97 :
98 : /// The corresponding Mesh in the neighbor
99 1 : Mesh<VolumeDim> operator()(const Mesh<VolumeDim>& mesh) const;
100 :
101 : /// An array whose elements are permuted such that
102 : /// `result[this->operator()(d)] = array_to_permute[d]`.
103 : ///
104 : /// \note the permutation depends only on how the dimension is mapped
105 : /// and ignores the side of the mapped direction.
106 : template <typename T>
107 1 : std::array<T, VolumeDim> permute_to_neighbor(
108 : const std::array<T, VolumeDim>& array_to_permute) const;
109 :
110 : /// An array whose elements are permuted such that
111 : /// `result[d] = array_in_neighbor[this->operator()(d)]`
112 : ///
113 : /// \note the permutation depends only on how the dimension is mapped
114 : /// and ignores the side of the mapped direction.
115 : template <typename T>
116 1 : std::array<T, VolumeDim> permute_from_neighbor(
117 : const std::array<T, VolumeDim>& array_in_neighbor) const;
118 :
119 : /// The corresponding Orientation of the host in the frame of the neighbor.
120 1 : OrientationMap<VolumeDim> inverse_map() const;
121 :
122 : /// Serialization for Charm++
123 : // NOLINTNEXTLINE(google-runtime-references)
124 1 : void pup(PUP::er& p);
125 :
126 : private:
127 0 : friend bool operator==(const OrientationMap& lhs, const OrientationMap& rhs) {
128 : return lhs.bit_field_ == rhs.bit_field_;
129 : }
130 :
131 0 : Direction<VolumeDim> get_direction(size_t dim) const;
132 0 : void set_direction(size_t dim, const Direction<VolumeDim>& direction);
133 0 : void set_aligned(bool is_aligned);
134 0 : std::set<size_t> set_of_dimensions() const;
135 :
136 0 : uint16_t bit_field_{0b1 << 15};
137 : };
138 :
139 : /// Output operator for OrientationMap.
140 : template <size_t VolumeDim>
141 1 : std::ostream& operator<<(std::ostream& os,
142 : const OrientationMap<VolumeDim>& orientation);
143 :
144 : template <size_t VolumeDim>
145 0 : bool operator!=(const OrientationMap<VolumeDim>& lhs,
146 : const OrientationMap<VolumeDim>& rhs) {
147 : return not(lhs == rhs);
148 : }
149 :
150 : template <size_t VolumeDim>
151 : template <typename T>
152 : std::array<T, VolumeDim> OrientationMap<VolumeDim>::permute_to_neighbor(
153 : const std::array<T, VolumeDim>& array_to_permute) const {
154 : std::array<T, VolumeDim> array_in_neighbor = array_to_permute;
155 : if (is_aligned() or VolumeDim <= 1) {
156 : return array_in_neighbor;
157 : }
158 : for (size_t i = 0; i < VolumeDim; i++) {
159 : gsl::at(array_in_neighbor, this->operator()(i)) =
160 : gsl::at(array_to_permute, i);
161 : }
162 : return array_in_neighbor;
163 : }
164 :
165 : template <size_t VolumeDim>
166 : template <typename T>
167 : std::array<T, VolumeDim> OrientationMap<VolumeDim>::permute_from_neighbor(
168 : const std::array<T, VolumeDim>& array_in_neighbor) const {
169 : std::array<T, VolumeDim> result = array_in_neighbor;
170 : if (not is_aligned() and VolumeDim > 1) {
171 : for (size_t i = 0; i < VolumeDim; i++) {
172 : gsl::at(result, i) = gsl::at(array_in_neighbor, this->operator()(i));
173 : }
174 : }
175 : return result;
176 : }
177 :
178 : /// \ingroup ComputationalDomainGroup
179 : /// `OrientationMap`s define an active rotation of the logical axes that bring
180 : /// the axes of a host block into alignment with the logical axes of the
181 : /// neighbor block. `discrete_rotation` applies this active rotation on the
182 : /// coordinates as opposed to the axes.
183 : /// For a two-dimensional example, consider a host block and a neighbor block,
184 : /// where the OrientationMap between them is \f$\{-\eta,+\xi\}\f$. A quarter-
185 : /// turn counterclockwise of the host block's logical axes would bring them into
186 : /// alignment with those of the neighbor. That is, after this active rotation,
187 : /// the blocks would be Aligned. Now consider a point A with coordinates
188 : /// (+1.0,-0.5). An active quarter-turn rotation counter-clockwise about the
189 : /// origin, keeping the axes fixed, brings point A into the coordinates
190 : /// (+0.5,+1.0). This is how `discrete_rotation` interprets the
191 : /// `OrientationMap` passed to it.
192 : template <size_t VolumeDim, typename T>
193 1 : std::array<T, VolumeDim> discrete_rotation(
194 : const OrientationMap<VolumeDim>& rotation,
195 : std::array<T, VolumeDim> source_coords);
196 :
197 : /*!
198 : * \ingroup ComputationalDomainGroup
199 : * \brief Computes the Jacobian of the transformation that is computed by
200 : * `discrete_rotation()`
201 : *
202 : * \note This always returns a `double` because the Jacobian is spatially
203 : * constant.
204 : */
205 : template <size_t VolumeDim>
206 1 : tnsr::Ij<double, VolumeDim, Frame::NoFrame> discrete_rotation_jacobian(
207 : const OrientationMap<VolumeDim>& orientation);
208 :
209 : /*!
210 : * \ingroup ComputationalDomainGroup
211 : * \brief Computes the inverse Jacobian of the transformation that is computed
212 : * by `discrete_rotation()`
213 : */
214 : template <size_t VolumeDim>
215 1 : tnsr::Ij<double, VolumeDim, Frame::NoFrame> discrete_rotation_inverse_jacobian(
216 : const OrientationMap<VolumeDim>& orientation);
|
