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