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 <limits>
9 : #include <optional>
10 :
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Domain/CoordinateMaps/Distribution.hpp"
13 : #include "Domain/Structure/OrientationMap.hpp"
14 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
15 :
16 : /// \cond
17 : namespace PUP {
18 : class er;
19 : } // namespace PUP
20 : /// \endcond
21 :
22 : namespace domain {
23 : namespace CoordinateMaps {
24 :
25 : /*!
26 : * \ingroup CoordinateMapsGroup
27 : *
28 : * \brief A reorientable map from the cube to a frustum.
29 : * \image html frustum_table1.png "A frustum with rectangular bases."
30 : *
31 : * ### Description and Specifiable Parameters
32 : * A map from the logical cube to the volume determined by interpolating between
33 : * two parallel rectangles each perpendicular to the \f$z\f$-axis. A Frustum map
34 : * \f$\vec{x}(\xi,\eta,\zeta)\f$ is determined by specifying two
35 : * heights \f$z_1 = \f$ `lower_bound` and \f$z_2 = \f$ `upper_bound` for the
36 : * positions of the rectangular bases of the frustum along the \f$z-\f$axis,
37 : * and the sizes of the rectangles are determined by the eight values passed to
38 : * `face_vertices`:
39 : *
40 : * \f{align*}{
41 : * &\textrm{lower x upper base} : x^{+\zeta}_{-\xi,-\eta} = x(-1,-1,1)\\
42 : * &\textrm{lower x lower base} : x^{-\zeta}_{-\xi,-\eta} = x(-1,-1,-1)\\
43 : * &\textrm{upper x upper base} : x^{+\zeta}_{+\xi,+\eta} = x(1,1,1)\\
44 : * &\textrm{upper x lower base} : x^{-\zeta}_{+\xi,+\eta} = x(1,1,-1)\\
45 : * &\textrm{lower y upper base} : y^{+\zeta}_{-\xi,-\eta} = y(-1,-1,1)\\
46 : * &\textrm{lower y lower base} : y^{-\zeta}_{-\xi,-\eta} = y(-1,-1,-1)\\
47 : * &\textrm{upper y upper base} : y^{+\zeta}_{+\xi,+\eta} = y(1,1,1)\\
48 : * &\textrm{upper y lower base} : y^{-\zeta}_{+\xi,+\eta} = y(1,1,-1)\f}
49 : *
50 : * \note As an example, consider a frustum along the z-axis, with the lower base
51 : * starting at (x,y) = (-2.0,3.0) and extending to
52 : * (2.0,5.0), and with the upper base extending from (0.0,1.0) to (1.0,3.0).
53 : * The corresponding value for `face_vertices` is `{{{{-2.0,3.0}}, {{2.0,5.0}},
54 : * {{0.0,1.0}}, {{1.0,3.0}}}}`.
55 : *
56 : * In the case where the two rectangles are geometrically similar, the volume
57 : * formed is a geometric frustum. However, this coordinate map generalizes to
58 : * rectangles which need not be similar.
59 : * The user may reorient the frustum by passing an OrientationMap to the
60 : * constructor. If `with_equiangular_map` is true, then this coordinate map
61 : * applies a tangent function mapping to the logical \f$\xi\f$ and \f$\eta\f$
62 : * coordinates. We then refer to the generalized
63 : * logical coordinates as \f$\Xi\f$ and \f$\mathrm{H}\f$. If
64 : * `projective_scale_factor` is set to a quantity other than unity, then this
65 : * coordinate map applies a rational function mapping to the logical \f$\zeta\f$
66 : * coordinate. This generalized logical coordinate is referred to as
67 : * \f$\mathrm{Z}\f$. If `auto_projective_scale_factor` is `true`, the
68 : * user-specified `projective_scale_factor` is ignored and an appropriate value
69 : * for `projective_scale_factor` is computed based on the values passed to
70 : * `face_vertices`. See the
71 : * [page on redistributing gridpoints](@ref redistributing_gridpoints)
72 : * to see more detailed information on equiangular variables and projective
73 : * scaling.
74 : *
75 : * ### Coordinate Map and Jacobian
76 : * In terms of the `face_vertices` variables, we define the
77 : * following auxiliary variables:
78 : *
79 : * \f{align*}{\Sigma^x &= \frac{1}{4}(x^{+\zeta}_{-\xi,-\eta} +
80 : * x^{+\zeta}_{+\xi,+\eta} + x^{-\zeta}_{-\xi,-\eta} + x^{-\zeta}_{+\xi,+\eta})
81 : * \\
82 : * \Delta^x_{\zeta} &= \frac{1}{4}(x^{+\zeta}_{-\xi,-\eta} +
83 : * x^{+\zeta}_{+\xi,+\eta} - x^{-\zeta}_{-\xi,-\eta} - x^{-\zeta}_{+\xi,+\eta})
84 : * \\
85 : * \Delta^x_{\xi} &= \frac{1}{4}(x^{+\zeta}_{+\xi,+\eta} -
86 : * x^{+\zeta}_{-\xi,-\eta} + x^{-\zeta}_{+\xi,+\eta} - x^{-\zeta}_{-\xi,-\eta})
87 : * \\
88 : * \Delta^x_{\xi\zeta} &= \frac{1}{4}(x^{+\zeta}_{+\xi,+\eta} -
89 : * x^{+\zeta}_{-\xi,-\eta} - x^{-\zeta}_{+\xi,+\eta} + x^{-\zeta}_{-\xi,-\eta})
90 : * \\
91 : * \Sigma^y &= \frac{1}{4}(y^{+\zeta}_{-\xi,-\eta} +
92 : * y^{+\zeta}_{+\xi,+\eta} + y^{-\zeta}_{-\xi,-\eta} + y^{-\zeta}_{+\xi,+\eta})
93 : * \\
94 : * \Delta^y_{\zeta} &= \frac{1}{4}(y^{+\zeta}_{-\xi,-\eta} +
95 : * y^{+\zeta}_{+\xi,+\eta} - y^{-\zeta}_{-\xi,-\eta} - y^{-\zeta}_{+\xi,+\eta})
96 : * \\
97 : * \Delta^y_{\eta} &= \frac{1}{4}(y^{+\zeta}_{+\xi,+\eta} -
98 : * y^{+\zeta}_{-\xi,-\eta} + y^{-\zeta}_{+\xi,+\eta} - y^{-\zeta}_{-\xi,-\eta})
99 : * \\
100 : * \Delta^y_{\eta\zeta} &= \frac{1}{4}(y^{+\zeta}_{+\xi,+\eta} -
101 : * y^{+\zeta}_{-\xi,-\eta} - y^{-\zeta}_{+\xi,+\eta} + y^{-\zeta}_{-\xi,-\eta})
102 : * \\
103 : * \Sigma^z &= \frac{z_1 + z_2}{2}\\
104 : * \Delta^z &= \frac{z_2 - z_1}{2}
105 : * \f}
106 : *
107 : * The full map is then given by:
108 : * \f[\vec{x}(\xi,\eta,\zeta) = \begin{bmatrix}
109 : * \Sigma^x + \Delta^x_{\xi}\Xi + (\Delta^x_{\zeta} +
110 : * \Delta^x_{\xi\zeta}\Xi)\mathrm{Z}\\
111 : * \Sigma^y + \Delta^y_{\eta}\mathrm{H} + (\Delta^y_{\zeta} +
112 : * \Delta^y_{\eta\zeta}\mathrm{H})\mathrm{Z}\\
113 : * \Sigma^z + \Delta^z\mathrm{Z}\\
114 : * \end{bmatrix}\f]
115 : *
116 : * With Jacobian: \f[J =
117 : * \begin{bmatrix}
118 : * (\Delta^x_{\xi} + \Delta^x_{\xi\zeta}\mathrm{Z})\Xi' &
119 : * 0 & (\Delta^x_{\zeta}+ \Delta^x_{\xi\zeta}\Xi)\mathrm{Z}' \\
120 : * 0 & (\Delta^y_{\eta} + \Delta^y_{\eta\zeta}\mathrm{Z})\mathrm{H}' &
121 : * (\Delta^y_{\zeta} + \Delta^y_{\eta\zeta}\mathrm{H})\mathrm{Z}'\\
122 : * 0 & 0 & \Delta^z\mathrm{Z}' \\
123 : * \end{bmatrix}
124 : * \f]
125 : *
126 : * ### Suggested values for the projective scale factor
127 : * \note This section assumes familiarity with projective scaling as discussed
128 : * on the [page on redistributing gridpoints](@ref redistributing_gridpoints).
129 : *
130 : * When constructing a Frustum map, it is not immediately obvious what value of
131 : * \f$w_{\delta}\f$ to use in the projective map; here we present a choice of
132 : * \f$w_{\delta}\f$ that will produce minimal grid distortion. We cover two
133 : * cases: the first is the special case in which the two rectangular Frustum
134 : * bases are related by a simple scale factor; the second is the general case.
135 : *
136 : * \image html ProjectionOfCube.png "Trapezoids from squares. (Davide Cervone)"
137 : *
138 : * As seen in Cervone's [Cubes and Hypercubes Rotating]
139 : * (http://www.math.union.edu/~dpvc/math/4D/rotation/welcome.html), there is a
140 : * special case in which the inverse projection of a trapezoid is not another
141 : * trapezoid, but a rectangle where the bases are congruent.
142 : * Most often one will want to use the special value of \f$w_{\delta}\f$ where
143 : * this occurs. This value of \f$w_{\delta}\f$ corresponds to the projective
144 : * transformation mapping a rectangular prism in an ambient 4D space with
145 : * congruent faces at \f$w=1\f$ and \f$w=w_{\delta}\f$ to the frustum in the
146 : * plane \f$w=1\f$:
147 : *
148 : * \f[w_{\delta} = \frac{L_1}{L_2}\f]
149 : *
150 : * where \f$\frac{L_1}{L_2}\f$ is the ratio between any pair of
151 : * corresponding side lengths of the \f$z_1\f$ and \f$z_2\f$-bases of the
152 : * frustum, respectively.
153 : *
154 : * For the general case one will want to use the value:
155 : * \f[w_{\delta} =
156 : * \frac{\sqrt{(x^{-\zeta}_{+\xi,+\eta} - x^{-\zeta}_{-\xi,+\eta})
157 : * (y^{-\zeta}_{+\xi,+\eta} - y^{-\zeta}_{+\xi,-\eta})}}
158 : * {\sqrt{(x^{+\zeta}_{+\xi,+\eta} - x^{+\zeta}_{-\xi,+\eta})
159 : * (y^{+\zeta}_{+\xi,+\eta} - y^{+\zeta}_{+\xi,-\eta})}}\f]
160 : *
161 : * This is the value for \f$w_{\delta}\f$ used by this CoordinateMap when
162 : * `auto_projective_scale_factor` is `true`.
163 : *
164 : * ### Bulged Frustum Coordinate Map and Jacobian
165 : * Each of the frustum faces in the frustum map given above are flat, but the
166 : * upper +z face of the frustum can be bulged out by setting a non-zero value
167 : * for the `sphericity`, where a value of `0.0` corresponds to the usual flat-
168 : * face, and a value of `1.0` corresponds to a value of fully spherical. Using
169 : * OrientationMaps allows the user to create a set of frustums that fully cover
170 : * a spherical surface. The radius of the sphere is determined by the corner of
171 : * the frustum that is furthest from the origin.
172 : *
173 : * The full map is given by:
174 : * \f[\vec{x}(\xi,\eta,\zeta) = \begin{bmatrix}
175 : * \Sigma^x + \Delta^x_{\xi}\Xi + (\Delta^x_{\zeta} +
176 : * \Delta^x_{\xi\zeta}\Xi)\mathrm{Z}\\
177 : * \Sigma^y + \Delta^y_{\eta}\mathrm{H} + (\Delta^y_{\zeta} +
178 : * \Delta^y_{\eta\zeta}\mathrm{H})\mathrm{Z}\\
179 : * \Sigma^z + \Delta^z\mathrm{Z}\\
180 : * \end{bmatrix}
181 : * + s \frac{1 +
182 : * \mathrm{Z}}{2}\left(\frac{R}{|\vec{\sigma}_{\mathrm{+z}}|}-1\right)
183 : * \vec{\sigma}_{\mathrm{+z}}
184 : * \f]
185 : *
186 : * where \f$R\f$ is the radius, \f$s\f$ is the `sphericity`, and
187 : * \f$\vec{\sigma}_{\mathrm{+z}}\f$ is given by:
188 : * \f[
189 : * \vec{\sigma}_{\mathrm{+z}} =
190 : * \begin{bmatrix}
191 : * \Sigma^x + \Delta^x_{\xi}\Xi + \Delta^x_{\zeta} +
192 : * \Delta^x_{\xi\zeta}\Xi\\
193 : * \Sigma^y + \Delta^y_{\eta}\mathrm{H} + \Delta^y_{\zeta} +
194 : * \Delta^y_{\eta\zeta}\mathrm{H}\\
195 : * \Sigma^z + \Delta^z\\
196 : * \end{bmatrix}.
197 : * \f]
198 : *
199 : * The Jacobian is:
200 : * \f[J =
201 : * \begin{bmatrix}
202 : * (\Delta^x_{\xi} + \Delta^x_{\xi\zeta}\mathrm{Z})\Xi' &
203 : * 0 & (\Delta^x_{\zeta}+ \Delta^x_{\xi\zeta}\Xi)\mathrm{Z}' \\
204 : * 0 & (\Delta^y_{\eta} + \Delta^y_{\eta\zeta}\mathrm{Z})\mathrm{H}' &
205 : * (\Delta^y_{\zeta} + \Delta^y_{\eta\zeta}\mathrm{H})\mathrm{Z}'\\
206 : * 0 & 0 & \Delta^z\mathrm{Z}' \\
207 : * \end{bmatrix}
208 : * +
209 : * \frac{s}{2} \left\{ \left(\frac{R}{|\vec{\sigma}_{\mathrm{+z}}|}-1\right)
210 : * \vec{\sigma}_{\mathrm{+z}}\hat{\zeta}^{\intercal}\mathrm{Z}'+
211 : * (1+\mathrm{Z})\left(\frac{R}{|\vec{\sigma}_{\mathrm{+z}}|}\left(
212 : * \mathbb{1}-\frac{\vec{\sigma}_{\mathrm{+z}}
213 : * \vec{\sigma}_{\mathrm{+z}}^{\intercal}}
214 : * {|\vec{\sigma}_{\mathrm{+z}}|^2}\right)-\mathbb{1}\right)
215 : * J_{\sigma}\right\},
216 : * \f]
217 : *
218 : * where \f$\hat{\zeta}\f$ is the row vector \f$[0, 0, 1]\f$, and
219 : * \f$J_{\sigma}\f$ is the Jacobian of the upper +z surface, given by:
220 : * \f[
221 : * J_{\sigma} =
222 : * \begin{bmatrix}
223 : * (\Delta^x_{\xi} + \Delta^x_{\xi\zeta})\Xi' &
224 : * 0 & 0 \\
225 : * 0 & (\Delta^y_{\eta} + \Delta^y_{\eta\zeta})\mathrm{H}' &
226 : * 0 \\
227 : * 0 & 0 & 0 \\
228 : * \end{bmatrix}.
229 : * \f]
230 : *
231 : * ### Using the Frustum Transition between Equiangular Maps
232 : * Using the parameter \f$\phi\f$, the user can specify whether to use the
233 : * full equiangular map on the upper +z face, or whether to use only the upper
234 : * or lower half of the equiangular map on the upper +z face. This capability
235 : * is used in the BinaryCompactObject Domain, where two equiangular maps around
236 : * each spherical object must be transitioned into a single equiangular map at
237 : * the outer boundary. The transition region that interpolates between these
238 : * two regions is the layer of Blocks with the Frustum maps.
239 : *
240 : * The full map is obtained by setting
241 : * \f$\Xi = \frac{1 + \mathrm{Z}}{2}\Xi_{\mathrm{upper}}(\xi,\zeta) +
242 : * \frac{1 - \mathrm{Z}}{2}\Xi_{0}(\xi)\f$
243 : *
244 : * for the volume map \f$\vec{x}(\xi,\eta,\zeta)\f$, and
245 : * \f$\Xi = \Xi_{\mathrm{upper}}\f$ for the surface map
246 : * \f$\vec{\sigma}(\xi,\eta)\f$, where \f$\Xi_{\mathrm{upper}}\f$ is the
247 : * \f$\phi\f$-dependent generalized equiangular map, which corresponds to the
248 : * lower half of the equiangular map when \f$\phi = -1\f$,
249 : * and to the upper half of the equiangular map when \f$\phi = 1 \f$.
250 : * It is given by:
251 : *
252 : * \f[\Xi_{\mathrm{upper}} = (1 +
253 : * \phi^2)\tan(\frac{\pi}{4}\frac{\xi+\phi}{1+\phi^2})-\phi\f]
254 : *
255 : * For \f$\phi=0\f$ this generalized map simplifies to the original equiangular
256 : * map \f$Xi_0 = \tan(\frac{\pi}{4}\xi)\f$.
257 : *
258 : * This function is compatible with the ability to bulge out the Frustum; the
259 : * map and Jacobian remain unchanged, modulo this substitution.
260 : */
261 1 : class Frustum {
262 : public:
263 0 : static constexpr size_t dim = 3;
264 : /*!
265 : * Constructs a frustum.
266 : * \param face_vertices An array of four 2D vertices that define the (x, y)
267 : * coordinates of the upper and lower base of the frustum.
268 : * \param lower_bound z distance from the origin to the lower base of the
269 : * frustum.
270 : * \param upper_bound z distance from the origin to the upper base of the
271 : * frustum.
272 : * \param orientation_of_frustum The orientation of the frustum in 3D space.
273 : * \param with_equiangular_map Determines whether to apply a tangent function
274 : * mapping to the logical coordinates (true) or not (false).
275 : * \param zeta_distribution Whether to apply a linear, logarithmic, or
276 : * projective mapping to the logical zeta coordinate.
277 : * \param distribution_value In the case of a projective map, the projective
278 : * scale factor, otherwise unused.
279 : * \param sphericity Value between 0 and 1 which determines whether the
280 : * surface of the upper base of the Frustum is flat (value of 0), spherical
281 : * (value of 1), or somewhere in between.
282 : * \param transition_phi Determines whether the equiangular map used conforms
283 : * to a lower half wedge (value of -1), an upper half wedge (value of +1), or
284 : * a full wedge (value of 0).
285 : * \param opening_angle The opening angle of the wedge the frustum will
286 : * conform to (have conforming boundaries with).
287 : */
288 1 : Frustum(const std::array<std::array<double, 2>, 4>& face_vertices,
289 : double lower_bound, double upper_bound,
290 : OrientationMap<3> orientation_of_frustum,
291 : bool with_equiangular_map = false,
292 : Distribution zeta_distribution = Distribution::Linear,
293 : std::optional<double> distribution_value = std::nullopt,
294 : double sphericity = 0.0, double transition_phi = 0.0,
295 : double opening_angle = M_PI_2);
296 0 : Frustum() = default;
297 0 : ~Frustum() = default;
298 0 : Frustum(Frustum&&) = default;
299 0 : Frustum(const Frustum&) = default;
300 0 : Frustum& operator=(const Frustum&) = default;
301 0 : Frustum& operator=(Frustum&&) = default;
302 :
303 : template <typename T>
304 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
305 : const std::array<T, 3>& source_coords) const;
306 :
307 : /// Returns std::nullopt if \f$z\f$ is at or beyond the \f$z\f$-coordinate of
308 : /// the apex of the pyramid, tetrahedron, or triangular prism that is
309 : /// formed by extending the `Frustum` (for a
310 : /// \f$z\f$-oriented `Frustum`).
311 : /// The inverse function is only callable with doubles because the inverse
312 : /// might fail if called for a point out of range, and it is unclear
313 : /// what should happen if the inverse were to succeed for some points in a
314 : /// DataVector but fail for other points.
315 1 : std::optional<std::array<double, 3>> inverse(
316 : const std::array<double, 3>& target_coords) const;
317 :
318 : template <typename T>
319 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
320 : const std::array<T, 3>& source_coords) const;
321 :
322 : template <typename T>
323 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
324 : const std::array<T, 3>& source_coords) const;
325 :
326 : // NOLINTNEXTLINE(google-runtime-references)
327 0 : void pup(PUP::er& p);
328 :
329 0 : bool is_identity() const { return is_identity_; }
330 :
331 : private:
332 0 : friend bool operator==(const Frustum& lhs, const Frustum& rhs);
333 :
334 0 : OrientationMap<3> orientation_of_frustum_{};
335 0 : bool with_equiangular_map_{false};
336 0 : bool is_identity_{false};
337 0 : Distribution zeta_distribution_ = Distribution::Linear;
338 0 : double zeta_distribution_value_{std::numeric_limits<double>::signaling_NaN()};
339 0 : double sigma_x_{std::numeric_limits<double>::signaling_NaN()};
340 0 : double delta_x_zeta_{std::numeric_limits<double>::signaling_NaN()};
341 0 : double delta_x_xi_{std::numeric_limits<double>::signaling_NaN()};
342 0 : double delta_x_xi_zeta_{std::numeric_limits<double>::signaling_NaN()};
343 0 : double sigma_y_{std::numeric_limits<double>::signaling_NaN()};
344 0 : double delta_y_zeta_{std::numeric_limits<double>::signaling_NaN()};
345 0 : double delta_y_eta_{std::numeric_limits<double>::signaling_NaN()};
346 0 : double delta_y_eta_zeta_{std::numeric_limits<double>::signaling_NaN()};
347 0 : double sigma_z_{std::numeric_limits<double>::signaling_NaN()};
348 0 : double delta_z_zeta_{std::numeric_limits<double>::signaling_NaN()};
349 0 : double w_plus_{std::numeric_limits<double>::signaling_NaN()};
350 0 : double w_minus_{std::numeric_limits<double>::signaling_NaN()};
351 0 : double sphericity_{std::numeric_limits<double>::signaling_NaN()};
352 0 : double radius_{std::numeric_limits<double>::signaling_NaN()};
353 0 : double phi_{std::numeric_limits<double>::signaling_NaN()};
354 0 : double half_opening_angle_{std::numeric_limits<double>::signaling_NaN()};
355 0 : double one_over_tan_half_opening_angle_{
356 : std::numeric_limits<double>::signaling_NaN()};
357 0 : double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
358 : };
359 :
360 0 : bool operator!=(const Frustum& lhs, const Frustum& rhs);
361 : } // namespace CoordinateMaps
362 : } // namespace domain
|
