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::CoordinateMaps {
23 :
24 : namespace detail {
25 : // This mapping can be deleted once the 2D and 3D wedges are oriented the same
26 : // (see issue https://github.com/sxs-collaboration/spectre/issues/2988)
27 : template <size_t Dim>
28 : struct WedgeCoordOrientation;
29 : template <>
30 : struct WedgeCoordOrientation<2> {
31 : static constexpr size_t radial_coord = 0;
32 : static constexpr size_t polar_coord = 1;
33 : static constexpr size_t azimuth_coord = 2; // unused
34 : };
35 : template <>
36 : struct WedgeCoordOrientation<3> {
37 : static constexpr size_t radial_coord = 2;
38 : static constexpr size_t polar_coord = 0;
39 : static constexpr size_t azimuth_coord = 1;
40 : };
41 : } // namespace detail
42 :
43 : /*!
44 : * \ingroup CoordinateMapsGroup
45 : *
46 : * \brief Map from a square or cube to a wedge.
47 : * \image html Shell.png "A shell can be constructed out of six wedges."
48 : *
49 : * \details The mapping that goes from a reference cube (in 3D) or square (in
50 : * 2D) to a wedge centered on a coordinate axis covering a volume between an
51 : * inner surface and outer surface. Each surface can be given a curvature
52 : * between flat (a sphericity of 0) or spherical (a sphericity of 1).
53 : *
54 : * In 2D, the first logical coordinate corresponds to the radial coordinate,
55 : * and the second logical coordinate corresponds to the angular coordinate. In
56 : * 3D, the first two logical coordinates correspond to the two angular
57 : * coordinates, and the third to the radial coordinate. This difference
58 : * originates from separate implementations for the 2D and 3D map that were
59 : * merged. The 3D implementation can be changed to use the first logical
60 : * coordinate as the radial direction, but this requires propagating the change
61 : * through the rest of the domain code (see issue
62 : * https://github.com/sxs-collaboration/spectre/issues/2988).
63 : *
64 : * The following documentation is for the **centered** 3D map, as we will defer
65 : * the dicussion of `Wedge`s with a `focal_offset_` to a later section. The 2D
66 : * map is obtained by setting either of the two angular coordinates to zero
67 : * (and using \f$\xi\f$ as the radial coordinate). Note that there is also a
68 : * normalization factor of $\sqrt{3}$ that appears in multiple expressions in
69 : * the 3D case that becomes $\sqrt{2}$ in the 2D case.
70 : *
71 : * The Wedge map is constructed by linearly interpolating between a bulged
72 : * face of radius `radius_inner_` to a bulged face of radius `radius_outer_`,
73 : * where the radius of each bulged face is defined to be the radius of the
74 : * sphere circumscribing the bulge.
75 : *
76 : * We make a choice here as to whether we wish to use the logical coordinates
77 : * parameterizing these surface as they are, in which case we have the
78 : * equidistant choice of coordinates, or whether to apply a tangent map to them
79 : * which leads us to the equiangular choice of coordinates. `Wedge`s have
80 : * variable `opening_angles_` which, for centered `Wedge`s, are the angular
81 : * sizes of the wedge in the $\xi$ and $\eta$ directions (for the 3D case) in
82 : * the target frame. By default, `Wedge`s have opening angles of $\pi/2$, so we
83 : * will discuss that case here and defer both the discussion of generalized
84 : * opening angles and the interaction between opening angles and non-zero focal
85 : * offsets for later sections.
86 : *
87 : * For a Wedge with $\xi$ and $\eta$ opening angles of $\pi/2$, the
88 : * equiangular coordinates in terms of the logical coordinates are:
89 : *
90 : * \begin{align}
91 : * \textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4)
92 : * \label{eq:equiangular_xi_pi_over_2}
93 : * \end{align}
94 : *
95 : * \begin{align}
96 : * \textrm{equiangular eta} :
97 : * \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)
98 : * \label{eq:equiangular_eta_pi_over_2}
99 : * \end{align}
100 : *
101 : * With derivatives:
102 : *
103 : * \begin{align}
104 : * \Xi'(\xi) &= \frac{\pi}{4}(1+\Xi^2) \\
105 : * \mathrm{H}'(\eta) &= \frac{\pi}{4}(1+\mathrm{H}^2)
106 : * \end{align}
107 : *
108 : * The equidistant coordinates are:
109 : *
110 : * \begin{align}
111 : * \textrm{equidistant xi} : \Xi = \xi \\
112 : * \textrm{equidistant eta} : \mathrm{H} = \eta
113 : * \end{align}
114 : *
115 : * with derivatives:
116 : *
117 : * \begin{align}
118 : * \Xi'(\xi) &= 1 \\
119 : * \mathrm{H}'(\eta) &= 1
120 : * \end{align}
121 : *
122 : * We also define the variable \f$\rho\f$, given by:
123 : *
124 : * \begin{align}
125 : * \textrm{rho} : \rho = \sqrt{1+\Xi^2+\mathrm{H}^2}
126 : * \end{align}
127 : *
128 : * ### The Spherical Face Map
129 : * The surface map for the spherical face of radius \f$R\f$ lying in the
130 : * \f$+z\f$ direction in either choice of coordinates is then given by:
131 : *
132 : * \begin{align}
133 : * \vec{\sigma}_{spherical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})
134 : * \end{align}
135 : *
136 : * Where
137 : *
138 : * \begin{align}
139 : * \vec{x}(\xi,\eta) =
140 : * \begin{bmatrix}
141 : * x(\xi,\eta) \\
142 : * y(\xi,\eta) \\
143 : * z(\xi,\eta) \\
144 : * \end{bmatrix} =
145 : * \frac{R}{\rho}
146 : * \begin{bmatrix}
147 : * \Xi \\
148 : * \mathrm{H} \\
149 : * 1 \\
150 : * \end{bmatrix}
151 : * \end{align}
152 : *
153 : * ### The Bulged Face Map
154 : * The bulged surface is itself constructed by linearly interpolating between
155 : * a cubical face and a spherical face. The surface map for the cubical face
156 : * of side length \f$2L\f$ lying in the \f$+z\f$ direction is given by:
157 : *
158 : * \begin{align}
159 : * \vec{\sigma}_{cubical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})
160 : * \end{align}
161 : *
162 : * Where
163 : *
164 : * \begin{align}
165 : * \vec{x}(\xi,\eta) =
166 : * \begin{bmatrix}
167 : * x(\xi,\eta) \\
168 : * y(\xi,\eta) \\
169 : * L \\
170 : * \end{bmatrix} =
171 : * L\begin{bmatrix}
172 : * \Xi \\
173 : * \mathrm{H} \\
174 : * 1 \\
175 : * \end{bmatrix}
176 : * \end{align}
177 : *
178 : * To construct the bulged map we interpolate between this cubical face map
179 : * and a spherical face map of radius \f$R\f$, with the interpolation
180 : * parameter being \f$s\f$, called the *sphericity* and which ranges from
181 : * 0 to 1, with 0 corresponding to a flat surface and 1 corresponding to a
182 : * spherical surface. The surface map for the bulged face lying in the \f$+z\f$
183 : * direction is then given by:
184 : *
185 : * \begin{align}
186 : * \vec{\sigma}_{bulged}(\xi,\eta) =
187 : * \left\{(1-s)L +
188 : * \frac{sR}{\rho}\right\}
189 : * \begin{bmatrix}
190 : * \Xi \\
191 : * \mathrm{H} \\
192 : * 1 \\
193 : * \end{bmatrix}
194 : * \end{align}
195 : *
196 : * We constrain $L$ by demanding that the spherical face circumscribe the cube.
197 : * With this condition, we have \f$L = R/\sqrt3\f$.
198 : * \note This differs from the choice in SpEC where it is demanded that the
199 : * surfaces touch at the cube face centers, which leads to \f$L = R\f$.
200 : *
201 : * ### The Full Volume Map
202 : * The final map for the wedge which lies along the \f$+z\f$ axis is obtained
203 : * by interpolating between the two surfaces with the interpolation parameter
204 : * being the logical coordinate \f$\zeta\f$. For a wedge whose gridpoints are
205 : * **linearly** distributed in the radial direction (`radial_distribution_` is
206 : * \ref domain::CoordinateMaps::Distribution
207 : * "domain::CoordinateMaps::Distribution::Linear"), this interpolation results
208 : * in the following map:
209 : *
210 : * \begin{align}
211 : * \vec{x}(\xi,\eta,\zeta) =
212 : * \frac{1}{2}\left\{
213 : * (1-\zeta)\Big[
214 : * (1-s_{inner})\frac{R_{inner}}{\sqrt 3} +
215 : * s_{inner}\frac{R_{inner}}{\rho}
216 : * \Big] +
217 : * (1+\zeta)\Big[
218 : * (1-s_{outer})\frac{R_{outer}}{\sqrt 3} +
219 : * s_{outer}\frac{R_{outer}}{\rho}
220 : * \Big]
221 : * \right\}
222 : * \begin{bmatrix}
223 : * \Xi \\
224 : * \mathrm{H} \\
225 : * 1 \\
226 : * \end{bmatrix}
227 : * \end{align}
228 : *
229 : * We will define the variables \f$F(\zeta)\f$ and \f$S(\zeta)\f$, the frustum
230 : * and sphere factors (in the linear case):
231 : *
232 : * \begin{align}
233 : * F(\zeta) &= F_0 + F_1\zeta \label{eq:frustum_factor} \\
234 : * S(\zeta) &= S_0 + S_1\zeta \label{eq:sphere_factor}
235 : * \end{align}
236 : *
237 : * Where
238 : *
239 : * \begin{align}
240 : * F_0 &=
241 : * \frac{1}{2} \big\{
242 : * (1-s_{outer})R_{outer} + (1-s_{inner})R_{inner}
243 : * \big\} \label{eq:frustum_zero_linear} \\
244 : * F_1 &= \partial_{\zeta}F
245 : * = \frac{1}{2} \big\{
246 : * (1-s_{outer})R_{outer} - (1-s_{inner})R_{inner}
247 : * \big\} \label{eq:frustum_rate_linear} \\
248 : * S_0 &=
249 : * \frac{1}{2} \big\{
250 : * s_{outer}R_{outer} + s_{inner}R_{inner}
251 : * \big\} \label{eq:sphere_zero_linear} \\
252 : * S_1 &= \partial_{\zeta}S
253 : * = \frac{1}{2} \big\{ s_{outer}R_{outer} - s_{inner}R_{inner}\big\}
254 : * \label{eq:sphere_rate_linear}
255 : * \end{align}
256 : *
257 : * The map can then be rewritten as:
258 : *
259 : * \begin{align}
260 : * \vec{x}(\xi,\eta,\zeta) =
261 : * \left\{
262 : * \frac{F(\zeta)}{\sqrt 3} + \frac{S(\zeta)}{\rho}
263 : * \right\}
264 : * \begin{bmatrix}
265 : * \Xi \\
266 : * \mathrm{H} \\
267 : * 1 \\
268 : * \end{bmatrix}
269 : * \end{align}
270 : *
271 : * The inverse map is given by:
272 : *
273 : * \begin{align}
274 : * \xi &= \frac{x}{z} \\
275 : * \eta &= \frac{y}{z} \\
276 : * \zeta &= \frac{z - \left(\frac{F_0}{\sqrt{3}} + \frac{S_0}{\rho}\right)}
277 : * {\left(\frac{F_1}{\sqrt{3}} + \frac{S_1}{\rho}\right)}
278 : * \end{align}
279 : *
280 : * We provide some common derivatives:
281 : *
282 : * \f{align}
283 : * \partial_{\xi}z &= \frac{-S(\zeta)\Xi\Xi'}{\rho^3} \\
284 : * \partial_{\eta}z &= \frac{-S(\zeta)\mathrm{H}\mathrm{H}'}{\rho^3} \\
285 : * \partial_{\zeta}z &= \frac{F'}{\sqrt 3} + \frac{S'(\zeta)}{\rho}
286 : * \f}
287 : *
288 : * The Jacobian then is:
289 : *
290 : * \begin{align}
291 : * J =
292 : * \begin{bmatrix}
293 : * \Xi'z + \Xi\partial_{\xi}z &
294 : * \Xi\partial_{\eta}z &
295 : * \Xi\partial_{\zeta}z \\
296 : * \mathrm{H}\partial_{\xi}z &
297 : * \mathrm{H}'z + \mathrm{H}\partial_{\eta}z &
298 : * \mathrm{H}\partial_{\zeta}z \\
299 : * \partial_{\xi}z &
300 : * \partial_{\eta}z &
301 : * \partial_{\zeta}z \\
302 : * \end{bmatrix}
303 : * \label{eq:jacobian_centered_wedge}
304 : * \end{align}
305 : *
306 : * A common factor that shows up in the inverse Jacobian is:
307 : *
308 : * \begin{align}
309 : * T:= \frac{S(\zeta)}{(\partial_{\zeta}z)\rho^3}
310 : * \end{align}
311 : *
312 : * The inverse Jacobian then is:
313 : * \f{align}
314 : * J^{-1} =
315 : * \frac{1}{z}\begin{bmatrix}
316 : * \Xi'^{-1} & 0 & -\Xi\Xi'^{-1} \\
317 : * 0 & \mathrm{H}'^{-1} & -\mathrm{H}\mathrm{H}'^{-1} \\
318 : * T\Xi & T\mathrm{H} & T + F(\partial_{\zeta}z)^{-1}/\sqrt 3 \\
319 : * \end{bmatrix}
320 : * \f}
321 : *
322 : * ### Changing the radial distribution of the gridpoints
323 : * By default, Wedge linearly distributes its gridpoints in the radial
324 : * direction. An exponential distribution of gridpoints can be obtained by
325 : * linearly interpolating in the logarithm of the radius in order to obtain
326 : * a relatively higher resolution at smaller radii. Since this is a radial
327 : * rescaling of Wedge, this option is only supported for fully spherical
328 : * wedges with `sphericity_inner_` = `sphericity_outer_` = 1.
329 : *
330 : * The linear interpolation done for a logarithmic radial distribution
331 : * (`radial_distribution_` is \ref domain::CoordinateMaps::Distribution
332 : * "domain::CoordinateMaps::Distribution::Logarithmic") is:
333 : *
334 : * \begin{align}
335 : * \ln r = \frac{1-\zeta}{2}\ln R_{inner} + \frac{1+\zeta}{2}\ln R_{outer}
336 : * \end{align}
337 : *
338 : * The map then is:
339 : *
340 : * \begin{align}
341 : * \vec{x}(\xi,\eta,\zeta) =
342 : * \frac{\sqrt{R_{inner}^{1-\zeta}R_{outer}^{1+\zeta}}}{\rho}
343 : * \begin{bmatrix}
344 : * \Xi \\
345 : * \mathrm{H} \\
346 : * 1 \\
347 : * \end{bmatrix}
348 : * \end{align}
349 : *
350 : * We can rewrite this map to take on the same form as the map for the linear
351 : * radial distribution, where we set
352 : *
353 : * \begin{align}
354 : * F(\zeta) &= 0 \\
355 : * S(\zeta) &= \sqrt{R_{inner}^{1-\zeta}R_{outer}^{1+\zeta}} \\
356 : * \end{align}
357 : *
358 : * Which gives us
359 : *
360 : * \begin{align}
361 : * \vec{x}(\xi,\eta,\zeta) =
362 : * \frac{S(\zeta)}{\rho}
363 : * \begin{bmatrix}
364 : * \Xi \\
365 : * \mathrm{H} \\
366 : * 1 \\
367 : * \end{bmatrix}
368 : * \end{align}
369 : *
370 : * The Jacobian then is still Eq. ($\ref{eq:jacobian_centered_wedge}$) but
371 : * where $F(\zeta)$ and $S(\zeta)$ are the quantities defined here for the
372 : * logarithmic distribution.
373 : *
374 : * Alternatively, an inverse radial distribution (`radial_distribution_` is
375 : * \ref domain::CoordinateMaps::Distribution
376 : * "domain::CoordinateMaps::Distribution::Inverse") can be chosen where the
377 : * linear interpolation is:
378 : *
379 : * \begin{align}
380 : * \frac{1}{r} =
381 : * \frac{R_\mathrm{inner} + R_\mathrm{outer}}
382 : * {2 R_\mathrm{inner}R_\mathrm{outer}} +
383 : * \frac{R_\mathrm{inner} - R_\mathrm{outer}}
384 : * {2R_\mathrm{inner} R_\mathrm{outer}} \zeta
385 : * \end{align}
386 : *
387 : * Which can be rewritten as:
388 : *
389 : * \begin{align}
390 : * \frac{1}{r} = \frac{1-\zeta}{2R_{inner}} + \frac{1+\zeta}{2R_{outer}}
391 : * \end{align}
392 : *
393 : * The map likewise takes the form:
394 : *
395 : * \begin{align}
396 : * \vec{x}(\xi,\eta,\zeta) =
397 : * \frac{S(\zeta)}{\rho}
398 : * \begin{bmatrix}
399 : * \Xi \\
400 : * \mathrm{H} \\
401 : * 1 \\
402 : * \end{bmatrix}
403 : * \end{align}
404 : *
405 : * Where
406 : *
407 : * \begin{align}
408 : * F(\zeta) &= 0 \\
409 : * S(\zeta) &=
410 : * \frac{2R_{inner}R_{outer}}
411 : * {(1 + \zeta)R_{inner} + (1 - \zeta)R_{outer}}
412 : * \end{align}
413 : *
414 : * Again, the Jacobian is still Eq. ($\ref{eq:jacobian_centered_wedge}$) but
415 : * where $F(\zeta)$ and $S(\zeta)$ are the quantities defined here for the
416 : * inverse distribution.
417 : *
418 : * ### Changing the opening angles
419 : * Consider the following map on \f$\xi \in [-1,1]\f$, which maps this interval
420 : * onto a parameterized curve that extends one fourth of a circle.
421 : *
422 : * \begin{align}
423 : * \vec{\Gamma}(\xi) =
424 : * \frac{R}{\sqrt{1+\xi^2}}
425 : * \begin{bmatrix}
426 : * 1 \\
427 : * \xi \\
428 : * \end{bmatrix}.
429 : * \label{eq:quarter_circle}
430 : * \end{align}
431 : *
432 : * It is convenient to compute the polar coordinate $\theta$ of the mapped
433 : * point as a function of $\xi$:
434 : *
435 : * \begin{align}
436 : * \theta(\xi) = \tan^{-1}\left(\frac{\Gamma_y(\xi)}{\Gamma_x(\xi)}\right).
437 : * \label{eq:polar_coord}
438 : * \end{align}
439 : *
440 : * The *opening angle* of the map is defined to be:
441 : *
442 : * \begin{align}
443 : * \Delta \theta = \theta(1) - \theta(-1),
444 : * \label{eq:define_opening_angle}
445 : * \end{align}
446 : *
447 : * We can see that with $\xi=\pm 1$, we have $\Gamma_x = R/\sqrt{2}$ and
448 : * $\Gamma_y=\pm R/\sqrt{2}$, giving us
449 : * $\theta(1) = \pi/4$ and $\theta(-1) = -\pi/4$. This wedge has an opening
450 : * angle $\pi/2$ radians, as expected.
451 : *
452 : * On the other hand, the following map has an opening angle of $\theta_O$:
453 : *
454 : * \begin{align}
455 : * \vec{\Gamma}(\xi) =
456 : * \frac{R}{\sqrt{1+\tan^2{(\theta_O/2)}\xi^2}}
457 : * \begin{bmatrix}
458 : * 1 \\
459 : * \tan{(\theta_O/2)}\xi \\
460 : * \end{bmatrix}.
461 : * \end{align}
462 : *
463 : * Let us also consider the generalized map
464 : *
465 : * \begin{align}
466 : * \vec{\Gamma}(\xi) =
467 : * \frac{R}{\sqrt{1+\Xi^2}}
468 : * \begin{bmatrix}
469 : * 1 \\
470 : * \Xi \\
471 : * \end{bmatrix},
472 : * \end{align}
473 : *
474 : * where $\Xi(\xi)$ is a function of $\xi$. $\theta(\xi)$ can then be written as
475 : *
476 : * \begin{align}
477 : * \theta(\xi) = \tan^{-1}(\Xi).
478 : * \label{eq:theta}
479 : * \end{align}
480 : *
481 : * For the map $\Xi(\xi) = \tan(\pi\xi/4)$, Eq. ($\ref{eq:theta}$) yields
482 : * $\theta(\xi) = \pi\xi/4$ and $\Delta\theta = \pi/2$. Note that this choice of
483 : * $\Xi(\xi)$ is equivalent to a reparameterization of the previous map given in
484 : * Eq. ($\ref{eq:quarter_circle}$). The reparameterization of the curve
485 : * $\vec{\Gamma}(\xi)$ via the tangent map yields an empirically superior
486 : * gridpoint distribution in practice. That this reparameterization should have
487 : * this property can be motivated by an observation of the following:
488 : *
489 : * \begin{align}
490 : * \frac{\mathrm{d}\tan^{-1}\Xi}{\mathrm{d}\xi}
491 : * = \frac{1}{1+\Xi^2}\frac{\mathrm{d}\Xi}{\mathrm{d}\xi}
492 : * = \frac{\pi}{4}.
493 : * \end{align}
494 : *
495 : * In other words, this parameterization has the property that the logical
496 : * coordinate $\xi$ subtends the angle $\theta$ at a constant rate. In general,
497 : * we say that a curve $\vec{\Gamma}(\xi)$ is parameterized *equiangularly* if
498 : *
499 : * \begin{align}
500 : * \frac{\mathrm{d}\theta}{\mathrm{d}\xi} = \text{const}.
501 : * \end{align}
502 : *
503 : * As for the map
504 : *
505 : * \begin{align}
506 : * \Xi(\xi) =
507 : * \tan{(\theta_O/2)}\frac{\tan{(\theta_D \xi/2)}}{\tan{(\theta_D/2)}},
508 : * \end{align}
509 : *
510 : * this choice of $\Xi(\xi)$ results in a $\vec{\Gamma}(\xi)$ with opening
511 : * angle $\theta_O$, which is equiangularly distributed if
512 : * $\theta_O = \theta_D$. In the Wedge map, the argument
513 : * `with_adapted_equiangular_map` controls whether to set
514 : * $\theta_O = \theta_D$ (the `true` case) or to set $\theta_D = \pi/2$
515 : * (the `false` case). When working with a 3D Wedge, the opening angles for the
516 : * Wedge can be separately controlled for both the $\xi$ and $\eta$ directions,
517 : * but `with_adapted_equiangular_map` will apply to both directions.
518 : * Additionally in the 3D case, it is not possible to set
519 : * `with_equiangular_map_` to `true` for all of the six wedges of a sphere
520 : * unless every opening angle is $\pi/2$. In the
521 : * \ref ::domain::creators::BinaryCompactObject "BinaryCompactObject" domain,
522 : * the outer $+y$, $-y$, $+z$, and $-z$ `Wedge`s are allowed to have a
523 : * user-specified opening angle in the $\xi$-direction, with a corresponding
524 : * $\theta_D$ equal to this opening angle, while in the $\eta$-direction the
525 : * opening angle is set to $\pi/2$. The two end cap `Wedge`s in the $+x$ and
526 : * $-x$ directions have angular dimensions and gridpoint distributions
527 : * determined by the other four `Wedge`s, as the six `Wedge`s must conforming
528 : * have gridpoint distributions at the $\xi = \pm1$, $\eta = \pm 1$ boundaries.
529 : *
530 : * ### Wedge with a Focal Offset
531 : * \image html FocalOffset.jpg "Wedges without and with a focal offset"
532 : *
533 : * In the case of the rectangular
534 : * \ref ::domain::creators::BinaryCompactObject "BinaryCompactObject" domain,
535 : * it becomes desirable to offset the center of the spherical excision surface
536 : * relative to the center of the cubical surface surrounding it. To enable the
537 : * offsetting of the central excision, the Wedge map must be generalized
538 : * according to the *focal lifting* method, which we will now discuss.
539 : *
540 : * We consider the problem of creating parameterized volumes from parameterized
541 : * surfaces. Consider a parameterized surface $\vec{\sigma}_{parent}(\xi,\eta)$,
542 : * also referred to as the *parent surface*. We define *focal lifting* as the
543 : * projection of this parent surface into a three-dimensional parameterized
544 : * volume $\vec{x}(\xi,\eta, \zeta)$ with respect to some *focus* $\vec{x}_0$
545 : * and *lifting scale factor* $\Lambda(\xi,\eta,\zeta)$. The resulting volume
546 : * is then said to be a *focally lifted* volume. These volume maps can be cast
547 : * into the following form:
548 : *
549 : * \begin{align}
550 : * \vec{x} - \vec{x}_0 = \Lambda(\vec{\sigma}_{parent}-\vec{x}_0),
551 : * \label{eq:focal_lifting}
552 : * \end{align}
553 : *
554 : * which makes apparent how the mapped point $\vec{x}(\xi,\eta,\zeta)$ is
555 : * obtained. The parametric equations for the generalized 3D Wedge maps can all
556 : * be written in the above form, which we will refer to as
557 : * *focally lifted form*. In the case of the 3D Wedge map with no focal offset,
558 : * we have:
559 : *
560 : * \begin{align}
561 : * \vec{x}_0 &= 0 \\
562 : * \Lambda &= \left\{\frac{F(\zeta)}{\sqrt{3}} +
563 : * \frac{S(\zeta)}{\rho} \right\} \\
564 : * \vec{\sigma}_{parent} &= \begin{bmatrix} \Xi, \mathrm{H}, 1 \end{bmatrix}^T
565 : * \end{align}
566 : *
567 : * The above map can be thought of as constructing a wedge from a biunit cube
568 : * centered at the origin. Points on the parent surface are scaled by a factor
569 : * of $\Lambda(\xi,\eta,\zeta)$ to obtain the corresponding point in the
570 : * volume. When generalizing the map to have a focus shifted from the origin
571 : * (obtained by setting `focal_offset_` to be non-zero), we scale the original
572 : * parent surface $\vec{\sigma}_{parent} = [\Xi, \mathrm{H},1]^T$ by a factor
573 : * $L$, and let the focus $\vec{x_0}$ shift away from the origin. The
574 : * generalized wedge map is then given by:
575 : *
576 : * \begin{align}
577 : * \vec{x} - \vec{x}_0 =
578 : * \left\{\frac{F(\zeta)}{L\sqrt 3} +
579 : * \frac{S(\zeta)}{L\rho}\right\}
580 : * \begin{bmatrix}
581 : * L\Xi - x_0 \\
582 : * L\mathrm{H} - y_0 \\
583 : * L-z_0 \\
584 : * \end{bmatrix}
585 : * \end{align}
586 : *
587 : * where we are now defining $\rho$ to be
588 : *
589 : * \begin{align}
590 : * \rho = \sqrt{(\Xi - x_0/L)^2 + (\mathrm{H} - y_0/L)^2 + (1 - z_0/L)^2}.
591 : * \label{eq:generalized_rho}
592 : * \end{align}
593 : *
594 : * This map is often written as:
595 : *
596 : * \begin{align}
597 : * \vec{x} - \vec{x}_0 =
598 : * \left\{\frac{F(\zeta)}{\sqrt{3}} +
599 : * \frac{S(\zeta)}{\rho}\right\}(\vec{\sigma}_0 - \vec{x}_0/L),
600 : * \label{eq:focally_lifted_map_with_s_and_f_factors}
601 : * \end{align}
602 : *
603 : * where $\vec{\sigma}_0 = [\Xi, \mathrm{H},1]^T$, as the parent surface
604 : * $\vec{\sigma}_{parent}$ is now $L\vec{\sigma}_0$. We give the quantity in
605 : * braces the name $z_{\Lambda} = L\Lambda$, *generalized z*. With this
606 : * definition, we can rewrite
607 : * Eq. ($\ref{eq:focally_lifted_map_with_s_and_f_factors}$) in the simpler form,
608 : *
609 : * \begin{align}
610 : * \vec{x} - \vec{x}_0 = z_{\Lambda}(\vec{\sigma}_0 - \vec{x}_0/L).
611 : * \label{eq:focally_lifted_map_with_generalized_z_coef}
612 : * \end{align}
613 : *
614 : * \note In the offset case, the frustum factor $F(\zeta)$ and sphere factor
615 : * $S(\zeta)$ (Eqs. ($\ref{eq:frustum_factor}$) and ($\ref{eq:sphere_factor}$))
616 : * for a linear radial distribution are no longer defined by the general $F_0$,
617 : * $F_1$, $S_0$, and $S_1$ given by Eqs.
618 : * ($\ref{eq:frustum_zero_linear}$), ($\ref{eq:frustum_rate_linear}$),
619 : * ($\ref{eq:sphere_zero_linear}$), and ($\ref{eq:sphere_rate_linear}$). In the
620 : * offset case, the inner surface must be spherical $(s_{inner} = 1)$ and the
621 : * outer surface can only be spherical or flat
622 : * $(s_{outer} = 0 \textrm{ or } s_{outer} = 1)$. In the case where
623 : * $s_{outer} = 0$, $L/\sqrt{3}$ is taken to be $R_{outer}$.
624 : *
625 : * The map can be inverted by first solving for \f$z_{\Lambda}\f$ in terms of
626 : * the target coordinates. We make use of the fact that the parent surface
627 : * $\vec{\sigma}_{parent}$ has a constant normal vector $\hat{n} = \hat{z}$.
628 : *
629 : * \begin{align}
630 : * z_{\Lambda} = \frac{(\vec{x} - \vec{x}_0)\cdot\hat{n}}
631 : * {(\vec{\sigma}_0-\vec{x}_0/L)\cdot\hat{n}}.
632 : * \end{align}
633 : *
634 : * In other words, when $\hat{n} = \hat{z}$,
635 : *
636 : * \begin{align}
637 : * z_{\Lambda} = \left\{\frac{F(\zeta)}{\sqrt{3}} +
638 : * \frac{S(\zeta)}{\rho}\right\}
639 : * = \frac{z - z_0}{1 - z_0/L}
640 : * \end{align}
641 : *
642 : * Moving all the known quantities in
643 : * Eq. ($\ref{eq:focally_lifted_map_with_generalized_z_coef}$) to the left hand
644 : * side results in the following expression that solves for the source
645 : * coordinates $\xi$ and $\eta$ in terms of the target coordinates:
646 : *
647 : * \begin{align}
648 : * \frac{\vec{x} - \vec{x}_0}{z_{\Lambda}} + \frac{\vec{x}_0}{L}
649 : * = \vec{\sigma}_0(\xi,\eta)
650 : * = \begin{bmatrix}
651 : * \Xi \\
652 : * \mathrm{H} \\
653 : * 1 \\
654 : * \end{bmatrix},
655 : * \end{align}
656 : *
657 : * Note that $|\vec{\sigma}_0 - \vec{x}_0/L| = \sqrt{(\Xi - x_0/L)^2 +
658 : * (\mathrm{H} - y_0/L)^2 + (1 - z_0/L)^2} = \rho$, indicating that an
659 : * expression for $\rho$ in terms of the target coordinates can be computed via
660 : * taking the magnitude of both sides of
661 : * Eq. ($\ref{eq:focally_lifted_map_with_generalized_z_coef}$):
662 : *
663 : * \begin{align}
664 : * |\vec{x} - \vec{x}_0| = z_{\Lambda}|\vec{\sigma}_0 - \vec{x}_0/L|
665 : * = z_{\Lambda}\rho.
666 : * \end{align}
667 : *
668 : * The quantity $\rho$ is then given by:
669 : *
670 : * \begin{align}
671 : * \rho = \frac{|\vec{x} - \vec{x}_0|}{z_{\Lambda}}.
672 : * \end{align}
673 : *
674 : * With $\rho$ computed, the radial source coordinate $\zeta$ can be computed
675 : * from
676 : *
677 : * \begin{align}
678 : * z_{\Lambda} = \left\{\frac{F(\zeta)}{\sqrt{3}} +
679 : * \frac{S(\zeta)}{\rho} \right\}
680 : * = \left\{\frac{F_0}{\sqrt{3}} + \frac{S_0}{\rho} +
681 : * \frac{F_1\zeta}{\sqrt{3}} + \frac{S_1\zeta}{\rho}\right\},
682 : * \end{align}
683 : *
684 : * which gives
685 : *
686 : * \begin{align}
687 : * \zeta = \frac{z_{\Lambda} -
688 : * \left(\frac{F_0}{\sqrt{3}} + \frac{S_0}{\rho}\right)}
689 : * {\left(\frac{F_1}{\sqrt{3}} + \frac{S_1}{\rho}\right)}.
690 : * \end{align}
691 : *
692 : * To compute the Jacobian, it is useful to first note that $\rho$
693 : * (Eq. ($\ref{eq:generalized_rho}$)) is the magnitude of the vector
694 : *
695 : * \begin{align}
696 : * \vec{\rho} = \vec{\sigma}_0 - \vec{x}_0/L
697 : * = \begin{bmatrix}
698 : * \Xi - x_0/L \\
699 : * \mathrm{H} - y_0/L \\
700 : * 1 - z_0/L
701 : * \end{bmatrix}
702 : * \end{align}
703 : *
704 : * and that we can express the target coordinates in
705 : * Eq. ($\ref{eq:focally_lifted_map_with_generalized_z_coef}$) in terms of the
706 : * components of $\vec{\rho}$:
707 : *
708 : * \begin{align}
709 : * x &= z_{\Lambda}\rho_x + x_0 \\
710 : * y &= z_{\Lambda}\rho_y + y_0 \\
711 : * z &= z_{\Lambda}\rho_z + z_0
712 : * \end{align}
713 : *
714 : * Some common terms used in the Jacobian are the derivatives of $z_{\Lambda}$
715 : * with respect to the source coordinates:
716 : *
717 : * \begin{align}
718 : * \partial_{\xi}z_{\Lambda} &=
719 : * \frac{-S(\zeta)\Xi'\rho_x}{\rho^3} \\
720 : * \partial_{\eta}z_{\Lambda} &=
721 : * \frac{-S(\zeta)\mathrm{H}'\rho_y}{\rho^3} \\
722 : * \partial_{\zeta}z_{\Lambda} &=
723 : * \frac{F'(\zeta)}{\sqrt{3}} + \frac{S'(\zeta)}{\rho}
724 : * \end{align}
725 : *
726 : * The Jacobian then is:
727 : *
728 : * \begin{align}
729 : * J =
730 : * \begin{bmatrix}
731 : * \Xi'z_{\Lambda} + \rho_x\partial_{\xi}z_{\Lambda} &
732 : * \rho_x\partial_{\eta}z_{\Lambda} &
733 : * \rho_x\partial_{\zeta}z_{\Lambda} \\
734 : * \rho_y\partial_{\xi}z_{\Lambda} &
735 : * \mathrm{H}'z_{\Lambda} + \rho_y\partial_{\eta}z_{\Lambda} &
736 : * \rho_y\partial_{\zeta}z_{\Lambda} \\
737 : * \rho_z\partial_{\xi}z_{\Lambda} &
738 : * \rho_z\partial_{\eta}z_{\Lambda} &
739 : * \rho_z\partial_{\zeta}z_{\Lambda} \\
740 : * \end{bmatrix}
741 : * \end{align}
742 : *
743 : * A common factor that shows up in this inverse Jacobian is:
744 : *
745 : * \begin{align}
746 : * T:= \frac{S(\zeta)}{(\partial_{\zeta}z_{\Lambda})\rho^3}
747 : * \end{align}
748 : *
749 : * And the inverse Jacobian is then:
750 : *
751 : * \begin{align}
752 : * J^{-1} =
753 : * \frac{1}{z_{\Lambda}}\begin{bmatrix}
754 : * \Xi'^{-1} & 0 & -\rho_x(\Xi'\rho_z)^{-1} \\
755 : * 0 & \mathrm{H}'^{-1} & -\rho_y(\mathrm{H}'\rho_z)^{-1} \\
756 : * T\rho_x & T\rho_y &
757 : * T\rho_z + F(\partial_{\zeta}z_{\Lambda}\rho_z)^{-1}/\sqrt{3}
758 : * \end{bmatrix}
759 : * \end{align}
760 : *
761 : * ### Offsetting a Rotated Wedge
762 : * The default Wedge map is oriented in the $+z$ direction, so the
763 : * construction of a Wedge oriented along a different direction requires an
764 : * additional OrientationMap $R$ to be passed to `orientation_of_wedge`. When
765 : * offsetting a rotated Wedge, the coordinates passed as parameters to
766 : * `focal_offset` are in the coordinate frame in which the Wedge is rotated
767 : * (the target frame). However, the focal lifting procedure (shown in
768 : * Eq. ($\ref{eq:focal_lifting}$)) is done in the default frame in which the
769 : * Wedge is facing the $+z$ direction, so the focal offset $\vec{x}_0$ is first
770 : * hit by the inverse rotation $R^{-1}$ and then the rotated focus
771 : * $R^{-1}\vec{x}_0$ is used internally as the focus for the $+z$ Wedge. When
772 : * the focal lifting calculation has completed, the rotation of the $+z$ Wedge
773 : * into the desired orientation by $R$ also rotates the focus into the desired
774 : * location. When performing the inverse operation, the focus is similarly
775 : * rotated into the default frame, where the inversion is performed.
776 : *
777 : * ### Interaction between opening angles and focal offsets
778 : * When a Wedge is created with a non-zero focal offset, the resulting shape
779 : * can take on a variety of possible angular sizes, depending on where the
780 : * focus is placed relative to the default centered location. The reader might
781 : * note that the angular size of a Wedge can also be controlled by passing an
782 : * argument to the `opening_angles` parameter in the Wedge constructor. While
783 : * both of these methods allow the angular size of a Wedge to be changed, the
784 : * user is prevented from employing both of them at the same time. In
785 : * particular, when the the offset is set to some non-zero value, the
786 : * `opening_angles_` member variable is set to $\pi/2$. Note that the
787 : * `opening_angles_` member being set to $\pi/2$ does not imply the
788 : * resulting Wedge will have an angular size of $\pi/2$. On the contrary, the
789 : * Wedge will have the angular size that is determined by the application of
790 : * the focal lifting method on the parent surface, which is the upper $+z$ face
791 : * of a cube that is centered at the origin.
792 : *
793 : * Because `opening_angles_` is set to $\pi/2$ when there is a non-zero focal
794 : * offset, when there is a non-zero focal offset and `with_equiangular_map_` is
795 : * `true`, $\Xi$ is given by Eq. ($\ref{eq:equiangular_xi_pi_over_2}$) and
796 : * $\mathrm{H}$ by Eq. ($\ref{eq:equiangular_eta_pi_over_2}$), just as it is
797 : * for the case of a centered Wedge with `opening_angles_` of $\pi/2$.
798 : */
799 : template <size_t Dim>
800 1 : class Wedge {
801 : public:
802 0 : static constexpr size_t dim = Dim;
803 0 : enum class WedgeHalves {
804 : /// Use the entire wedge
805 : Both,
806 : /// Use only the upper logical half
807 : UpperOnly,
808 : /// Use only the lower logical half
809 : LowerOnly
810 : };
811 :
812 : /*!
813 : * \brief Constructs a centered wedge (one with no focal offset)
814 : *
815 : * \param radius_inner Distance from the origin to one of the corners which
816 : * lie on the inner surface.
817 : * \param radius_outer Distance from the origin to one of the corners which
818 : * lie on the outer surface.
819 : * \param orientation_of_wedge The orientation of the desired wedge relative
820 : * to the orientation of the default wedge which is a wedge that has its
821 : * curved surfaces pierced by the upper-z axis. The logical $\xi$ and $\eta$
822 : * coordinates point in the cartesian x and y directions, respectively.
823 : * \param sphericity_inner Value between 0 and 1 which determines
824 : * whether the inner surface is flat (value of 0), spherical (value of 1) or
825 : * somewhere in between.
826 : * \param sphericity_outer Value between 0 and 1 which determines
827 : * whether the outer surface is flat (value of 0), spherical (value of 1) or
828 : * somewhere in between.
829 : * \param with_equiangular_map Determines whether to apply a tangent function
830 : * mapping to the logical coordinates (for `true`) or not (for `false`).
831 : * \param halves_to_use Determines whether to construct a full wedge or only
832 : * half a wedge. If constructing only half a wedge, the resulting shape has a
833 : * face normal to the x direction (assuming default OrientationMap). If
834 : * constructing half a wedge, an intermediate affine map is applied to the
835 : * logical xi coordinate such that the interval [-1,1] is mapped to the
836 : * corresponding logical half of the wedge. For example, if `UpperOnly` is
837 : * specified, [-1,1] is mapped to [0,1], and if `LowerOnly` is specified,
838 : * [-1,1] is mapped to [-1,0]. The case of `Both` means a full wedge, with no
839 : * intermediate map applied. In all cases, the logical points returned by the
840 : * inverse map will lie in the range [-1,1] in each dimension. Half wedges are
841 : * currently only useful in constructing domains for binary systems.
842 : * \param radial_distribution Determines how to distribute gridpoints along
843 : * the radial direction. For wedges that are not exactly spherical, only
844 : * `Distribution::Linear` is currently supported.
845 : * \param opening_angles Determines the angular size of the wedge. The default
846 : * value is $\pi/2$, which corresponds to a wedge size of $\pi/2$. For this
847 : * setting, four Wedges can be put together to cover $2\pi$ in angle along a
848 : * great circle. This option is meant to be used with the equiangular map
849 : * option turned on.
850 : * \param with_adapted_equiangular_map Determines whether to adapt the
851 : * point distribution in the wedge to match its physical angular size. When
852 : * `true`, angular distances are proportional to logical distances. Note
853 : * that it is not possible to use adapted maps in every Wedge of a Sphere
854 : * unless each Wedge has the same size along both angular directions.
855 : */
856 1 : Wedge(double radius_inner, double radius_outer, double sphericity_inner,
857 : double sphericity_outer, OrientationMap<Dim> orientation_of_wedge,
858 : bool with_equiangular_map,
859 : WedgeHalves halves_to_use = WedgeHalves::Both,
860 : Distribution radial_distribution = Distribution::Linear,
861 : const std::array<double, Dim - 1>& opening_angles =
862 : make_array<Dim - 1>(M_PI_2),
863 : bool with_adapted_equiangular_map = true);
864 :
865 : /*!
866 : * \brief Constructs a wedge with a focal offset
867 : *
868 : * \details Can construct an offset Wedge with a spherical inner surface and
869 : * either a spherical or a flat outer surface. If `radius_outer` has a value,
870 : * a spherical Wedge will be constructed, and if not, a flat one will be
871 : * constructed.
872 : *
873 : * Note that because the focal offset is what determines the angular size of
874 : * the Wedge, opening angles cannot be used with offset Wedges.
875 : *
876 : * In the event that `focal_offset` happens to be zero, the Wedge's member
877 : * variables and behavior will be set up to be equivalent to that of a
878 : * centered Wedge:
879 : * - `cube_half_length` will be ignored
880 : * - if `radius_outer` is `std::nullopt`, the outer radius of the Wedge will
881 : * be set to $\sqrt{\mathrm{Dim}}L$, where $L$ is the `cube_half_length`
882 : * - the opening angles ($\theta_O$) and opening angles distribution
883 : * ($\theta_D$) used will be $\pi/2$
884 : *
885 : * \param radius_inner Distance from the origin to one of the corners which
886 : * lie on the inner surface.
887 : * \param radius_outer If this has a value, it creates a spherical Wedge
888 : * (inner and outer sphericity are 1) where this is the distance from the
889 : * origin to one of the corners that lie on the inner surface. If this is
890 : * `std::nullopt`, it creates a Wedge with a flat outer surface
891 : * (inner sphericity is 1 and outer sphericity is 0). In the event that
892 : * `radius_outer == std::nullopt` **and** `focal_offset` is zero,
893 : * the outer radius will instead be set to $\sqrt{\mathrm{Dim}}L$, where $L$
894 : * is the `cube_half_length`. The outer radius is given a value in this
895 : * circumstance so that it can be handled as a centered Wedge (one with no
896 : * offset).
897 : * \param orientation_of_wedge The orientation of the desired wedge relative
898 : * to the orientation of the default wedge which is a wedge that has its
899 : * curved surfaces pierced by the upper-z axis. The logical $\xi$ and $\eta$
900 : * coordinates point in the cartesian x and y directions, respectively.
901 : * \param cube_half_length Half the length of the parent surface (see Wedge
902 : * documentation for more details). If `focal_offset` is zero, this
903 : * parameter has no effect and is ignored so that the Wedge can be handled
904 : * as a centered Wedge (one with no offset).
905 : * \param focal_offset The target frame coordinates of the focus from which
906 : * the Wedge is focally lifted.
907 : * \param with_equiangular_map Determines whether to apply a tangent function
908 : * mapping to the logical coordinates (for `true`) or not (for `false`).
909 : * \param halves_to_use Determines whether to construct a full wedge or only
910 : * half a wedge. If constructing only half a wedge, the resulting shape has a
911 : * face normal to the x direction (assuming default OrientationMap). If
912 : * constructing half a wedge, an intermediate affine map is applied to the
913 : * logical xi coordinate such that the interval [-1,1] is mapped to the
914 : * corresponding logical half of the wedge. For example, if `UpperOnly` is
915 : * specified, [-1,1] is mapped to [0,1], and if `LowerOnly` is specified,
916 : * [-1,1] is mapped to [-1,0]. The case of `Both` means a full wedge, with no
917 : * intermediate map applied. In all cases, the logical points returned by the
918 : * inverse map will lie in the range [-1,1] in each dimension. Half wedges are
919 : * currently only useful in constructing domains for binary systems.
920 : * \param radial_distribution Determines how to distribute gridpoints along
921 : * the radial direction. For wedges that are not exactly spherical, only
922 : * `Distribution::Linear` is currently supported.
923 : */
924 1 : Wedge(double radius_inner, std::optional<double> radius_outer,
925 : double cube_half_length, std::array<double, Dim> focal_offset,
926 : OrientationMap<Dim> orientation_of_wedge, bool with_equiangular_map,
927 : WedgeHalves halves_to_use = WedgeHalves::Both,
928 : Distribution radial_distribution = Distribution::Linear);
929 :
930 0 : Wedge() = default;
931 0 : ~Wedge() = default;
932 0 : Wedge(Wedge&&) = default;
933 0 : Wedge(const Wedge&) = default;
934 0 : Wedge& operator=(const Wedge&) = default;
935 0 : Wedge& operator=(Wedge&&) = default;
936 :
937 : template <typename T>
938 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
939 : const std::array<T, Dim>& source_coords) const;
940 :
941 : /// For a \f$+z\f$-oriented `Wedge`, returns invalid if \f$z<=0\f$
942 : /// or if \f$(x,y,z)\f$ is on or outside the cone defined
943 : /// by \f$(x^2/z^2 + y^2/z^2+1)^{1/2} = -S/F\f$, where
944 : /// \f$S = \frac{1}{2}(s_1 r_1 - s_0 r_0)\f$ and
945 : /// \f$F = \frac{1}{2\sqrt{3}}((1-s_1) r_1 - (1-s_0) r_0)\f$.
946 : /// Here \f$s_0,s_1\f$ and \f$r_0,r_1\f$ are the specified sphericities
947 : /// and radii of the inner and outer \f$z\f$ surfaces. The map is singular on
948 : /// the cone and on the xy plane.
949 : /// The inverse function is only callable with doubles because the inverse
950 : /// might fail if called for a point out of range, and it is unclear
951 : /// what should happen if the inverse were to succeed for some points in a
952 : /// DataVector but fail for other points.
953 1 : std::optional<std::array<double, Dim>> inverse(
954 : const std::array<double, Dim>& target_coords) const;
955 :
956 : template <typename T>
957 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
958 : const std::array<T, Dim>& source_coords) const;
959 :
960 : template <typename T>
961 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
962 : const std::array<T, Dim>& source_coords) const;
963 :
964 : // NOLINTNEXTLINE(google-runtime-references)
965 0 : void pup(PUP::er& p);
966 :
967 0 : static constexpr bool is_identity() { return false; }
968 :
969 : private:
970 : // maps between 2D and 3D choices for coordinate axis orientations
971 0 : static constexpr size_t radial_coord =
972 : detail::WedgeCoordOrientation<Dim>::radial_coord;
973 0 : static constexpr size_t polar_coord =
974 : detail::WedgeCoordOrientation<Dim>::polar_coord;
975 0 : static constexpr size_t azimuth_coord =
976 : detail::WedgeCoordOrientation<Dim>::azimuth_coord;
977 :
978 : /*!
979 : * \brief Factors out the calculation of \f$\Xi(\xi)\f$ and $\mathrm{H}$
980 : *
981 : * \details The **equidistant** parametrization
982 : * (when `with_equiangular_map_ == false`) of the logical coordinates is
983 : *
984 : * \f{align*}{
985 : * \Xi(\xi) = \xi.
986 : * \f}
987 : *
988 : * The **equiangular** reparametrization
989 : * (when `with_equiangular_map_ == true`) of the logical coordinates is
990 : *
991 : * \f{align*}{
992 : * \Xi(\xi) =
993 : * \tan{(\theta_O/2)}\frac{\tan{(\theta_D \xi/2)}}{\tan{(\theta_D/2)}},
994 : * \f}
995 : *
996 : * where $\theta_O$ (element of `opening_angles_`) and $\theta_D$
997 : * (element of `opening_angles_distribution_`) are described in the Wedge
998 : * class documentation.
999 : *
1000 : * When `focal_offset_` is nonzero, the **equiangular** reparametrization
1001 : * is instead
1002 : *
1003 : * \f{align*}{
1004 : * \Xi(\xi) = \tan{(\pi/4)}\xi
1005 : * \f}
1006 : *
1007 : * \tparam FuncIsXi whether the logical cooridnate `lowercase_xi_or_eta` is
1008 : * $\xi$ (polar coordinate) or $\eta$ (azimuthal coordinate)
1009 : * \param lowercase_xi_or_eta the logical coordinate $\xi$ or $\eta$ to map
1010 : */
1011 : template <bool FuncIsXi, typename T>
1012 1 : tt::remove_cvref_wrap_t<T> get_cap_angular_function(
1013 : const T& lowercase_xi_or_eta) const;
1014 :
1015 : /*!
1016 : * \brief Factors out the calculation of \f$\Xi'(\xi)\f$ and $\mathrm{H}'$
1017 : *
1018 : * \details Computes the derivatives of the quantities defined in
1019 : * `get_cap_angular_function()`.
1020 : *
1021 : * \tparam FuncIsXi whether the logical cooridnate `lowercase_xi_or_eta` is
1022 : * $\xi$ (polar coordinate) or $\eta$ (azimuthal coordinate)
1023 : * \param lowercase_xi_or_eta the logical coordinate $\xi$ or $\eta$ to map
1024 : */
1025 : template <bool FuncIsXi, typename T>
1026 1 : tt::remove_cvref_wrap_t<T> get_deriv_cap_angular_function(
1027 : const T& lowercase_xi_or_eta) const;
1028 :
1029 : /*!
1030 : * \brief Factors out the calculation of $\vec{\rho}$
1031 : *
1032 : * \details Computes
1033 : * \f{align*}{
1034 : * \vec{\rho} = [\Xi-x_0/L, \mathrm{H}-y_0/L, 1-z_0/L]^T
1035 : * \f}
1036 : *
1037 : * where \f$\Xi\f$ and $\mathrm{H}$ are the logical coordinate maps defined in
1038 : * `get_cap_angular_function()` and the Wedge class documentation,
1039 : * \f$\vec{x_0} = [x_0, y_0, z_0]^T\f$ is the result of applying the inverse
1040 : * map of the `orientation_of_wedge_` on the `focal_offset_`, and $L$ is the
1041 : * `cube_half_length_`.
1042 : *
1043 : * \param rotated_focus the result of applying the inverse map of the
1044 : * `orientation_of_wedge_` on the `focal_offset_`
1045 : * \param cap the function(s) \f$\Xi\f$ (and $\mathrm{H}$ in 3D)
1046 : */
1047 : template <typename T>
1048 1 : std::array<tt::remove_cvref_wrap_t<T>, Dim> get_rho_vec(
1049 : const std::array<double, Dim>& rotated_focus,
1050 : const std::array<tt::remove_cvref_wrap_t<T>, Dim - 1>& cap) const;
1051 :
1052 : /*!
1053 : * \brief Factors out the calculation of $1/\rho$
1054 : *
1055 : * \details Computes $1/\rho$ where
1056 : *
1057 : * \f{align*}{
1058 : * \rho = \sqrt{(\Xi - x_0/L)^2 + (\mathrm{H} - y_0/L)^2 + (1 - z_0/L)^2}.
1059 : * \f}
1060 : *
1061 : * Here, \f$\Xi\f$ and $\mathrm{H}$ are the logical coordinate maps defined in
1062 : * `get_cap_angular_function()` and the Wedge class documentation,
1063 : * \f$\vec{x_0} = [x_0, y_0, z_0]^T\f$ is the result of applying the inverse
1064 : * map of the `orientation_of_wedge_` on the `focal_offset_`, and $L$ is the
1065 : * `cube_half_length_`.
1066 : *
1067 : * \param rotated_focus the result of applying the inverse map of the
1068 : * `orientation_of_wedge_` on the `focal_offset_`
1069 : * \param cap the function(s) \f$\Xi\f$ (and $\mathrm{H}$ in 3D)
1070 : */
1071 : template <typename T>
1072 1 : tt::remove_cvref_wrap_t<T> get_one_over_rho(
1073 : const std::array<double, Dim>& rotated_focus,
1074 : const std::array<tt::remove_cvref_wrap_t<T>, Dim - 1>& cap) const;
1075 :
1076 : /*!
1077 : * \brief Factors out the calculation of $S(\zeta)$ needed for the map and the
1078 : * Jacobian
1079 : *
1080 : * \details The value of $S(\zeta)$ is computed differently for different
1081 : * radial distributions.
1082 : *
1083 : * For a **linear** radial distribution:
1084 : *
1085 : * \f{align*}{
1086 : * S(\zeta) = S_0 + S_1\zeta
1087 : * \f}
1088 : *
1089 : * where $S_0$ and $S_1$ are defined as
1090 : *
1091 : * \f{align*}{
1092 : * S_0 &=
1093 : * \frac{1}{2} \big\{
1094 : * s_{outer}R_{outer} + s_{inner}R_{inner}
1095 : * \big\} \\
1096 : * S_1 &= \partial_{\zeta}S
1097 : * = \frac{1}{2} \big\{ s_{outer}R_{outer} - s_{inner}R_{inner}\big\}
1098 : * \f}
1099 : *
1100 : * and are stored in `sphere_zero_` and `sphere_rate_`, respectively.
1101 : *
1102 : * For a **logarithmic** radial distribution:
1103 : *
1104 : * \f{align*}{
1105 : * S(\zeta) = \exp{(S_0 + S_1\zeta)}
1106 : * \f}
1107 : *
1108 : * where $S_0$ and $S_1$ are defined as
1109 : *
1110 : * \f{align*}{
1111 : * S_0 &= \frac{1}{2} \ln(R_{outer}R_{inner}) \\
1112 : * S_1 &= \frac{1}{2} \ln(R_{outer}/R_{inner})
1113 : * \f}
1114 : *
1115 : * With these definitions of $S_0$ and $S_1$, we can rewrite the expression
1116 : * for $S(\zeta)$ as:
1117 : *
1118 : * \f{align*}{
1119 : * S(\zeta) &= \sqrt{R_{inner}^{1-\zeta}R_{outer}^{1+\zeta}}
1120 : * \f}
1121 : *
1122 : * As with the linear distribution, $S_0$ and $S_1$ are stored in
1123 : * `sphere_zero_` and `sphere_rate_`, respectively.
1124 : *
1125 : * For an **inverse** radial distribution:
1126 : *
1127 : * \f{align*}{
1128 : * S(\zeta) =
1129 : * \frac{2R_{inner}R_{outer}}
1130 : * {(1 + \zeta)R_{inner} + (1 - \zeta)R_{outer}}
1131 : * \f}
1132 : *
1133 : * In this case, `sphere_zero_` and `sphere_rate_` will simply be `NaN`.
1134 : *
1135 : * See Wedge for more details on these quantities.
1136 : *
1137 : * \param zeta the radial source coordinate
1138 : */
1139 : template <typename T>
1140 1 : tt::remove_cvref_wrap_t<T> get_s_factor(const T& zeta) const;
1141 : /*!
1142 : * \brief Factors out the calculation of $S'(\zeta)$ needed for the Jacobian
1143 : *
1144 : * \details The value of $S'(\zeta)$ is computed differently for different
1145 : * radial distributions.
1146 : *
1147 : * For a **linear** radial distribution:
1148 : *
1149 : * \f{align*}{
1150 : * S'(\zeta) =
1151 : * \frac{1}{2} \big\{ s_{outer}R_{outer} - s_{inner}R_{inner}\big\}
1152 : * \f}
1153 : *
1154 : * For a **logarithmic** radial distribution:
1155 : *
1156 : * \f{align*}{
1157 : * S'(\zeta) = \frac{1}{2} S(\zeta)\ln(R_{outer}/R_{inner})
1158 : * \f}
1159 : *
1160 : * where $S(\zeta)$ is defined in `get_s_factor()`.
1161 : *
1162 : * For an **inverse** radial distribution:
1163 : *
1164 : * \f{align*}{
1165 : * S'(\zeta) =
1166 : * \frac{2(R_{inner} R_{outer}^2 - R_{inner}^2 R_{outer})}
1167 : * {(R_{inner} + R_{outer} + \zeta(R_{inner} - R_{outer}))^2}
1168 : * \f}
1169 : *
1170 : * See Wedge and `get_s_factor()` for more details on these quantities.
1171 : *
1172 : * \param zeta the radial source coordinate
1173 : * \param s_factor $S(\zeta)$ (see `get_s_factor()`)
1174 : */
1175 : template <typename T>
1176 1 : tt::remove_cvref_wrap_t<T> get_s_factor_deriv(const T& zeta,
1177 : const T& s_factor) const;
1178 :
1179 : /*!
1180 : * \brief Factors out the calculation of $z_{\Lambda}$ needed for the map and
1181 : * the Jacobian
1182 : *
1183 : * \details The value of $z_{\Lambda}$ is computed differently for different
1184 : * radial distributions.
1185 : *
1186 : * For a **linear** radial distribution:
1187 : *
1188 : * \f{align*}{
1189 : * z_{\Lambda} = \frac{F(\zeta)}{\sqrt 3} + \frac{S(\zeta)}{\rho}
1190 : * \f}
1191 : *
1192 : * For a **logarithmic** or **inverse** radial distribution:
1193 : *
1194 : * \f{align*}{
1195 : * z_{\Lambda} = \frac{S(\zeta)}{\rho}
1196 : * \f}
1197 : *
1198 : * See Wedge and `get_s_factor()` for more details on these quantities.
1199 : *
1200 : * \param zeta the radial source coordinate
1201 : * \param one_over_rho one over $\rho$ where
1202 : * $\rho = |\vec{\sigma}_0 - \vec{x}_0/L| = \sqrt{(\Xi - x_0/L)^2 +
1203 : * (\mathrm{H} - y_0/L)^2 + (1 - z_0/L)^2}$ (see Wedge)
1204 : * \param s_factor $S(\zeta)$ (see `get_s_factor()`)
1205 : */
1206 : template <typename T>
1207 1 : tt::remove_cvref_wrap_t<T> get_generalized_z(const T& zeta,
1208 : const T& one_over_rho,
1209 : const T& s_factor) const;
1210 : template <typename T>
1211 0 : tt::remove_cvref_wrap_t<T> get_generalized_z(const T& zeta,
1212 : const T& one_over_rho) const;
1213 : /*!
1214 : * \brief Factors out the calculation of $\partial_i z_{\Lambda}$ needed for
1215 : * the Jacobian
1216 : *
1217 : * \details For **all** radial distributions:
1218 : *
1219 : * \f{align*}{
1220 : * \partial_{\xi} z_{\Lambda} &=
1221 : * \frac{-S(\zeta)\Xi'\rho_x}{\rho^3} \\
1222 : * \partial_{\eta} z_{\Lambda} &=
1223 : * \frac{-S(\zeta)\mathrm{H}'\rho_y}{\rho^3} \\
1224 : * \partial_{\zeta} z_{\Lambda} &=
1225 : * \frac{F'(\zeta)}{\sqrt 3} + \frac{S'(\zeta)}{\rho}
1226 : * \f}
1227 : *
1228 : * However, $\partial_{\zeta} z_{\Lambda}$ reduces to
1229 : *
1230 : * \f{align*}{
1231 : * \partial_{\zeta} z_{\Lambda} &= \frac{S'(\zeta)}{\rho}
1232 : * \f}
1233 : *
1234 : * for **logarithmic** and **inverse** radial distributions because
1235 : * $F(\zeta) = 0$.
1236 : *
1237 : * See Wedge and `get_s_factor()` for more details on these quantities.
1238 : *
1239 : * \param zeta the radial source coordinate
1240 : * \param one_over_rho one over $\rho$ where
1241 : * $\rho = |\vec{\sigma}_0 - \vec{x}_0/L| = \sqrt{(\Xi - x_0/L)^2 +
1242 : * (\mathrm{H} - y_0/L)^2 + (1 - z_0/L)^2}$ (see Wedge)
1243 : * \param s_factor $S(\zeta)$ (see `get_s_factor()`)
1244 : * \param cap_deriv $\Xi'$ and $\mathrm{H}'$ (see Wedge)
1245 : * \param rho_vec $\vec{\rho} = [\Xi-x_0/L, \mathrm{H}-y_0/L, 1-z_0/L]^T$
1246 : * (see Wedge)
1247 : */
1248 : template <typename T>
1249 1 : std::array<tt::remove_cvref_wrap_t<T>, Dim> get_d_generalized_z(
1250 : const T& zeta, const T& one_over_rho, const T& s_factor,
1251 : const std::array<tt::remove_cvref_wrap_t<T>, Dim - 1>& cap_deriv,
1252 : const std::array<tt::remove_cvref_wrap_t<T>, Dim>& rho_vec) const;
1253 :
1254 : template <size_t LocalDim>
1255 : // NOLINTNEXTLINE(readability-redundant-declaration)
1256 0 : friend bool operator==(const Wedge<LocalDim>& lhs,
1257 : const Wedge<LocalDim>& rhs);
1258 :
1259 : /// Distance from the origin to one of the corners which lie on the inner
1260 : /// surface.
1261 1 : double radius_inner_{std::numeric_limits<double>::signaling_NaN()};
1262 : /// If this contains a value, it is the distance from the `focal_offset` to
1263 : /// one of the corners that lie on the outer surface. Set to `std::nullopt`
1264 : /// when `focal_offset` is nonzero and the outer surface is flat, because
1265 : /// there is no single outer radius like there is for a centered Wedge or a
1266 : /// spherical offset Wedge.
1267 1 : std::optional<double> radius_outer_ = std::nullopt;
1268 : /// Value between 0 and 1 which determines whether the inner surface is flat
1269 : /// (value of 0), spherical (value of 1) or somewhere in between. If
1270 : /// `focal_offset` is nonzero, `sphericity_inner` must be `1.0`.
1271 1 : double sphericity_inner_{std::numeric_limits<double>::signaling_NaN()};
1272 : /// Value between 0 and 1 which determines whether the outer surface is flat
1273 : /// (value of 0), spherical (value of 1) or somewhere in between. If
1274 : /// `focal_offset` is nonzero, `sphericity_outer` must be `0.0` or `1.0`.
1275 1 : double sphericity_outer_{std::numeric_limits<double>::signaling_NaN()};
1276 : /// Half the length of the parent surface (see Wedge documentation for more
1277 : /// details). This parameter has no effect and is set to `std::nullopt` when
1278 : /// `focal_offset` is zero.
1279 1 : std::optional<double> cube_half_length_ = std::nullopt;
1280 : /// The target frame coordinates of the focus from which the Wedge is focally
1281 : /// lifted.
1282 1 : std::array<double, Dim> focal_offset_{
1283 : make_array<Dim>(std::numeric_limits<double>::signaling_NaN())};
1284 : /// The orientation of the desired wedge relative to the orientation of the
1285 : /// default wedge which is a wedge that has its curved surfaces pierced by the
1286 : /// upper-z axis. The logical $\xi$ and $\eta$ coordinates point in the
1287 : /// cartesian x and y directions, respectively.
1288 1 : OrientationMap<Dim> orientation_of_wedge_ =
1289 : OrientationMap<Dim>::create_aligned();
1290 : /// Determines whether to apply a tangent function mapping to the logical
1291 : /// coordinates (for `true`) or not (for `false`).
1292 1 : bool with_equiangular_map_ = false;
1293 : /// Determines whether to construct a full wedge or only half a wedge (see
1294 : /// Wedge documentation for more details)
1295 1 : WedgeHalves halves_to_use_ = WedgeHalves::Both;
1296 : /// Determines how to distribute gridpoints along the radial direction. For
1297 : /// wedges that are not exactly spherical, only `Distribution::Linear` is
1298 : /// currently supported.
1299 1 : Distribution radial_distribution_ = Distribution::Linear;
1300 : /// $F_0 / \sqrt{3}$ (see Wedge documentation)
1301 1 : double scaled_frustum_zero_{std::numeric_limits<double>::signaling_NaN()};
1302 : /// $S_0$ (see Wedge documentation)
1303 1 : double sphere_zero_{std::numeric_limits<double>::signaling_NaN()};
1304 : /// $F_1 / \sqrt{3}$ (see Wedge documentation)
1305 1 : double scaled_frustum_rate_{std::numeric_limits<double>::signaling_NaN()};
1306 : /// $S_1$ (see Wedge documentation)
1307 1 : double sphere_rate_{std::numeric_limits<double>::signaling_NaN()};
1308 : /// $\theta_O$ (see Wedge documentation). Set to `std::nullopt` when
1309 : /// `focal_offset_` is nonzero.
1310 1 : std::optional<std::array<double, Dim - 1>> opening_angles_ = std::nullopt;
1311 : /// $\theta_D$ (see Wedge documentation). Set to `std::nullopt` when
1312 : /// `focal_offset_` is nonzero.
1313 1 : std::optional<std::array<double, Dim - 1>> opening_angles_distribution_ =
1314 : std::nullopt;
1315 : };
1316 :
1317 : template <size_t Dim>
1318 0 : bool operator!=(const Wedge<Dim>& lhs, const Wedge<Dim>& rhs);
1319 : } // namespace domain::CoordinateMaps
|
