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