Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines the class UniformCylindricalSide.
6 :
7 : #pragma once
8 :
9 : #include <array>
10 : #include <cstddef>
11 : #include <limits>
12 : #include <optional>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 :
16 : /// \cond
17 : namespace PUP {
18 : class er;
19 : } // namespace PUP
20 : /// \endcond
21 :
22 : namespace domain::CoordinateMaps {
23 :
24 : /*!
25 : * \ingroup CoordinateMapsGroup
26 : *
27 : * \brief Map from 3D unit right cylindrical shell to a volume that connects
28 : * portions of two spherical surfaces.
29 : *
30 : * \image html UniformCylSide.svg "A hollow cylinder maps to the shaded region."
31 : *
32 : * \details Consider two spheres with centers \f$C_1\f$ and \f$C_2\f$,
33 : * and radii \f$R_1\f$ and \f$R_2\f$. Sphere 1 is assumed to be contained
34 : * inside Sphere 2.
35 : * Let sphere 1 be intersected by two
36 : * planes normal to the \f$z\f$ axis and located at
37 : * \f$z = z^{\pm}_{\mathrm{P}1}\f$,
38 : * and let sphere 2 be intersected by two planes normal to the \f$z\f$ axis and
39 : * located at \f$z = z^{\pm}_{\mathrm{P}2}\f$. Here we assume that
40 : * \f$z^{-}_{\mathrm{P}2} \leq z^{-}_{\mathrm{P}1}<
41 : * z^{+}_{\mathrm{P}1} \leq z^{+}_{\mathrm{P}2}\f$.
42 : *
43 : * UniformCylindricalSide maps a 3D unit right cylindrical shell (with
44 : * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
45 : * \f$-1\leq\bar{z}\leq 1\f$ and \f$1 \leq \bar{x}^2+\bar{y}^2 \leq 4\f$, where
46 : * the values of 1 and 2 for the inner and outer cylindrical radii
47 : * are arbitrary choices but are required by UniformCylindricalSide)
48 : * to the shaded area in the figure above (with coordinates
49 : * \f$(x,y,z)\f$). The "inner surface" of the cylindrical shell
50 : * \f$\bar{x}^2+\bar{y}^2=1\f$ is mapped to the portion of sphere 1
51 : * that has \f$z^{-}_{\mathrm{P}1} \leq z \leq z^{+}_{\mathrm{P}1} \f$,
52 : * and on this portion of the sphere the cosine of the polar angular coordinate
53 : * \f$\cos\theta_1 =(z-C_1^z)/R_1\f$ is uniform in \f$\bar{z}\f$,
54 : * and the angular coordinate \f$\phi_1 = \atan((y-C_1^y)/(x-C_1^x))\f$
55 : * is the same as \f$\phi = \atan(\bar{y}/\bar{x})\f$.
56 : * Likewise, the "outer surface" of the cylindrical shell
57 : * \f$\bar{x}^2+\bar{y}^2=4\f$ is mapped to the portion of sphere 2
58 : * that has \f$z^{-}_{\mathrm{P}2} \leq z \leq z^{+}_{\mathrm{P}2}
59 : * \f$, and on this portion of the sphere the cosine of the azimuthal
60 : * angular coordinate
61 : * \f$\cos\theta_2 = (z-C_2^z)/R_2\f$ is uniform in \f$\bar{z}\f$,
62 : * and the angular coordinate \f$\phi_2 =
63 : * \atan((y-C_2^y)/(x-C_2^x))\f$ is the same as \f$\phi\f$.
64 : *
65 : * UniformCylindricalSide is different from CylindricalSide
66 : * because of the distribution of points on the spheres, and because
67 : * for UniformCylindricalSide the mapped portion of both Sphere 1
68 : * and Sphere 2 are bounded by planes of constant \f$z\f$, whereas for
69 : * CylindricalSide only one of the mapped portions is bounded by a
70 : * plane (except for specially chosen map parameters). Note that
71 : * UniformCylindricalSide can be used to construct maps that connect
72 : * an arbitrary number of nested spheres; this is not possible for
73 : * CylindricalSide for more than 3 nested spheres because of this
74 : * asymmetry between CylindricalSide's two spherical surfaces.
75 : *
76 : * Note that the entire region between Sphere 1 and Sphere 2 can be covered
77 : * by a single cylindrical shell (mapped using UniformCylindricalSide) and
78 : * two cylinders (each mapped by UniformCylindricalEndcap).
79 : *
80 : * UniformCylindricalSide is intended to be composed with `Wedge<2>` maps to
81 : * construct a portion of a cylindrical domain for a binary system.
82 : *
83 : * UniformCylindricalSide can be used to construct a domain that is similar
84 : * to, but not identical to, the one described briefly in the Appendix of
85 : * \cite Buchman:2012dw.
86 : * UniformCylindricalSide is used to construct the Blocks analogous to
87 : * those labeled 'CA cylinder', 'EA cylinder', 'CB cylinder', 'EE cylinder',
88 : * and 'EB cylinder' in Figure 20 of that paper.
89 : *
90 : * UniformCylindricalSide provides the following functions:
91 : *
92 : * ## operator()
93 : *
94 : * `operator()` maps \f$(\bar{x},\bar{y},\bar{z})\f$ to \f$(x,y,z)\f$
95 : * according to
96 : *
97 : * \f{align}
98 : * x &= C_1^x+\lambda(C_2^x-C_1^x) +
99 : * \cos\phi\left(R_1\sin\theta_1 +
100 : * \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x0} \\
101 : * y &= C_1^y+\lambda(C_2^y-C_1^y) +
102 : * \sin\phi\left(R_1\sin\theta_1 +
103 : * \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x1} \\
104 : * z &= C_1^z+\lambda(C_2^z-C_1^z) +
105 : * R_1\cos\theta_1 +
106 : * \lambda(R_2\cos\theta_2-R_1\cos\theta_1) \label{eq:x2}.
107 : * \f}
108 : *
109 : * Here
110 : * \f{align}
111 : * \lambda &= \bar{\rho}-1,\label{eq:lambdafromrhobar}\\
112 : * \cos\theta_1 &= \cos\theta_{1 \mathrm{max}} +
113 : * \left(\cos\theta_{1 \mathrm{min}}-\cos\theta_{1 \mathrm{max}}\right)
114 : * \frac{\bar{z}+1}{2}\label{eq:deftheta1}\\
115 : * \cos\theta_2 &= \cos\theta_{2 \mathrm{max}} +
116 : * \left(\cos\theta_{2 \mathrm{min}}-
117 : * \cos\theta_{2 \mathrm{max}}\right)
118 : * \frac{\bar{z}+1}{2}\label{eq:deftheta2}\\
119 : * \phi &= \atan(\bar{y}/\bar{x})\label{eq:defphi},
120 : * \f}
121 : * where \f$\theta_{1 \mathrm{min}}\f$, \f$\theta_{2 \mathrm{min}}\f$,
122 : * \f$\theta_{1 \mathrm{max}}\f$, and \f$\theta_{2 \mathrm{max}}\f$
123 : * are defined by
124 : * \f{align}
125 : * \label{eq:deftheta1min}
126 : * \cos(\theta_{1\mathrm{min}}) &= (z^{+}_{\mathrm{P}1}-C_1^z)/R_1,\\
127 : * \cos(\theta_{1\mathrm{max}}) &= (z^{-}_{\mathrm{P}1}-C_1^z)/R_1,\\
128 : * \cos(\theta_{2\mathrm{min}}) &= (z^{+}_{\mathrm{P}2}-C_2^z)/R_2,\\
129 : * \label{eq:deftheta2max}
130 : * \cos(\theta_{2\mathrm{max}}) &= (z^{-}_{\mathrm{P}2}-C_2^z)/R_2,
131 : * \f}
132 : * and
133 : * \f{align}
134 : * \bar{\rho} &= \sqrt{\bar{x}^2+\bar{y}^2}/\bar{R} \label{eq:defrhobar},
135 : * \f}
136 : * where \f$\bar{R}\f$ is the inner radius of the cylindrical shell in barred
137 : * coordinates, which is always unity.
138 : *
139 : * Note that \f$\theta_{1\mathrm{min}}<\theta_{1\mathrm{max}}\f$ but
140 : * \f$\cos\theta_{1\mathrm{min}}>\cos\theta_{1\mathrm{max}}\f$ (and same
141 : * for sphere 2).
142 : *
143 : * Also note that Eqs. (\f$\ref{eq:deftheta1}\f$) and
144 : * (\f$\ref{eq:deftheta2}\f$) can be simplified using Eqs.
145 : * (\f$\ref{eq:deftheta1min}\f$-\f$\ref{eq:deftheta2max}\f$):
146 : * \f{align}
147 : * R_1\cos\theta_1 &= z^{-}_{\mathrm{P}1}-C_1^z
148 : * +(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})
149 : * \frac{\bar{z}+1}{2}\label{eq:deftheta1alt}\\
150 : * R_2\cos\theta_2 &= z^{-}_{\mathrm{P}2}-C_2^z
151 : * +(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})
152 : * \frac{\bar{z}+1}{2}\label{eq:deftheta2alt}\\
153 : * \f}
154 : *
155 : * ## inverse
156 : *
157 : * Given \f$(x,y,z)\f$ we want to find \f$(\bar{x},\bar{y},\bar{z})\f$.
158 : * From Eqs. (\f$\ref{eq:x2}\f$), (\f$\ref{eq:deftheta1alt}\f$), and
159 : * (\f$\ref{eq:deftheta2alt}\f$) we can write \f$\bar{z}\f$ as a function
160 : * of \f$\lambda\f$:
161 : *
162 : * \f{align}
163 : * \frac{1+\bar{z}}{2} &=
164 : * \frac{z +
165 : * \lambda (z^{-}_{\mathrm{P}1}-z^{-}_{\mathrm{P}2}) - z^{-}_{\mathrm{P}1}}
166 : * {(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})
167 : * + \lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})}
168 : * \label{eq:zbar_from_lambda},
169 : * \f}
170 : * Note that the denominator of
171 : * Eq. (\f$\ref{eq:zbar_from_lambda}\f$) is always positive because
172 : * \f$0\leq\lambda\leq 1\f$, \f$z^{+}_{\mathrm{P}1}>z^{-}_{\mathrm{P}1}\f$,
173 : * and \f$z^{+}_{\mathrm{P}2}>z^{-}_{\mathrm{P}2}\f$.
174 : *
175 : * By eliminating \f$\phi\f$ from Eqs. (\f$\ref{eq:x0}\f$) and
176 : * (\f$\ref{eq:x1}\f$) we find that \f$\lambda\f$ is the solution
177 : * of \f$Q(\lambda)=0\f$, where
178 : *
179 : * \f{align}
180 : * Q(\lambda) &= \left(x-C_1^x-\lambda(C_2^x-C_1^x)\right)^2+
181 : * \left(y-C_1^y-\lambda(C_2^y-C_1^y)\right)^2-
182 : * \left((1-\lambda)R_1\sin\theta_1 +
183 : * \lambda R_2\sin\theta_2\right)^2.\label{eq:defQ}
184 : * \f}
185 : * Here \f$\theta_1\f$ and \f$\theta_2\f$ are functions
186 : * of \f$\bar{z}\f$ through Eqs. (\f$\ref{eq:deftheta1alt}\f$) and
187 : * (\f$\ref{eq:deftheta2alt}\f$), and \f$\bar{z}\f$ is a function of
188 : * \f$\lambda\f$ through Eq. (\f$\ref{eq:zbar_from_lambda}\f$).
189 : *
190 : * We solve \f$Q(\lambda)=0\f$ numerically; it is a one-dimensional
191 : * root-finding problem.
192 : *
193 : * Once we have determined \f$\lambda\f$, we then obtain \f$\bar{z}\f$
194 : * from Eq. (\f$\ref{eq:zbar_from_lambda}\f$), and we obtain \f$\phi\f$ from
195 : *
196 : * \f{align}
197 : * \tan\phi &=
198 : * \frac{y-C_1^y-\lambda(C_2^y-C_1^y)}{x-C_1^x-\lambda(C_2^x-C_1^x)}.
199 : * \f}
200 : *
201 : * Then \f$\bar{\rho}\f$ is obtained from Eq. (\f$\ref{eq:lambdafromrhobar}\f$)
202 : * and \f$\bar{x}\f$ and \f$\bar{y}\f$ are obtained from
203 : *
204 : * \f{align}
205 : * \bar{x} &= \bar{\rho}\bar{R}\cos\phi,\\
206 : * \bar{y} &= \bar{\rho}\bar{R}\sin\phi.
207 : * \f}
208 : *
209 : * ### Considerations when root-finding.
210 : *
211 : * We solve \f$Q(\lambda)=0\f$ numerically for \f$\lambda\f$,
212 : * where \f$Q(\lambda)\f$ is given by Eq. (\f$\ref{eq:defQ}\f$).
213 : *
214 : * #### min/max values of \f$\lambda\f$:
215 : *
216 : * Note that the root we care about must have \f$-1\leq\bar{z}\leq 1\f$;
217 : * therefore from Eq. (\f$\ref{eq:zbar_from_lambda}\f$) we have
218 : *
219 : * \f{align}
220 : * \lambda_{\mathrm{min}} &=
221 : * \mathrm{max}\left\{0,
222 : * \frac{z-z^{+}_{\mathrm{P}1}}
223 : * {(z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1})},
224 : * \frac{z^{-}_{\mathrm{P}1}-z}
225 : * {(z^{-}_{\mathrm{P}1}-z^{-}_{\mathrm{P}2})}
226 : * \right\}\label{eq:lambdamin}
227 : * \f}
228 : * In the case where \f$z^{+}_{\mathrm{P}2}=z^{+}_{\mathrm{P}1}\f$
229 : * we treat the middle term in Eq.(\f$\ref{eq:lambdamin}\f$) as zero since
230 : * in that case \f$z-z^{+}_{\mathrm{P}1}\f$ can never be positive for
231 : * \f$x^2\f$ in the range of the map, and for
232 : * \f$z=z^{+}_{\mathrm{P}2}=z^{+}_{\mathrm{P}1}\f$
233 : * it turns out that
234 : * (\f$\ref{eq:zbar_from_lambda}\f$) places no restriction on
235 : * \f$\lambda_{\mathrm{min}}\f$. For the same reason, if
236 : * \f$z^{-}_{\mathrm{P}2}=z^{-}_{\mathrm{P}1}\f$ we treat
237 : * the last term in Eq.(\f$\ref{eq:lambdamin}\f$) as zero.
238 : *
239 : * We look for a root only between \f$\lambda_{\mathrm{min}}\f$
240 : * and \f$\lambda_{\mathrm{max}}=1\f$.
241 : *
242 : * ##### Roots within roundoff of min or max \f$\lambda\f$
243 : *
244 : * Sometimes a root is within roundoff of \f$\lambda_{\mathrm{min}}\f$.
245 : * In this case, the root might not be bracketed by
246 : * \f$[\lambda_{\mathrm{min}},\lambda_{\mathrm{max}}]\f$ if the root
247 : * is slightly outside that interval by roundoff error. If we find that
248 : * \f$Q(\lambda_{\mathrm{min}})\f$ is near zero but has the wrong sign,
249 : * then we slightly expand the interval as follows:
250 : *
251 : * \f{align}
252 : * \lambda_{\mathrm{min}} \to \lambda_{\mathrm{min}}
253 : * - 2 \frac{Q(\lambda_{\mathrm{min}})}{Q'(\lambda_{\mathrm{min}})},
254 : * \f}
255 : *
256 : * where \f$Q'(\lambda_{\mathrm{min}})\f$ is the derivative of the function
257 : * in Eq. (\f$\ref{eq:defQ}\f$). Note that without the factor of 2, this is
258 : * a Newton-Raphson step; the factor of 2 is there to overcompensate so that
259 : * the new \f$\lambda_{\mathrm{min}}\f$ brackets the root.
260 : *
261 : * Note that by differentiating Eqs. (\f$\ref{eq:defQ}\f$) and
262 : * (\f$\ref{eq:zbar_from_lambda}\f$), one obtains
263 : *
264 : * \f{align}
265 : * Q'(\lambda) =& -2 \left[
266 : * \left(x-C_1^x-\lambda(C_2^x-C_1^x)\right)(C_2^x-C_1^x)+
267 : * \left(y-C_1^y-\lambda(C_2^y-C_1^y)\right)(C_2^y-C_1^y)
268 : * \right]\nonumber \\
269 : * &
270 : * -\left[
271 : * 2(R_2\sin\theta_2-R_1\sin\theta_1)
272 : * -(1-\lambda)\cot\theta_1 (z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})
273 : * \frac{d\bar{z}}{d\lambda} \right. \nonumber \\
274 : * & \left.\qquad
275 : * -\lambda \cot\theta_2 (z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})
276 : * \frac{d\bar{z}}{d\lambda}
277 : * \right]
278 : * \left((1-\lambda)R_1\sin\theta_1 +
279 : * \lambda R_2\sin\theta_2\right), \label{eq:defQderiv}
280 : * \f}
281 : *
282 : * where
283 : * \f{align}
284 : * \frac{d\bar{z}}{d\lambda} &=
285 : * \frac{(1-\bar{z})(z^{-}_{\mathrm{P}1}-z^{-}_{\mathrm{P}2})
286 : * -(1+\bar{z})(z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1})}
287 : * {(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})
288 : * + \lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})}
289 : * \label{eq:dzbar_dlambda}.
290 : * \f}
291 : *
292 : * A root within roundoff of \f$\lambda_{\mathrm{max}}\f$ is treated
293 : * similarly.
294 : *
295 : * #### Special cases:
296 : *
297 : * For some points on the boundary of the mapped domain,
298 : * \f$\lambda_{\mathrm{min}}\f$ will be within roundoff of
299 : * \f$\lambda=1\f$. We check explicitly for this case, and we
300 : * compute the root as exactly \f$\lambda=1\f$.
301 : *
302 : * ### Quick rejection of points out of range of the map.
303 : *
304 : * It is expected that `inverse()` will often be passed points
305 : * \f$(x,y,z)\f$ that are out of the range of the map; in this case
306 : * `inverse()` returns a `std::nullopt`. To avoid the difficulty and
307 : * expense of attempting to solve \f$Q(\lambda)=0\f$ numerically
308 : * for such points (and then having this solution fail), it is useful
309 : * to quickly reject points \f$(x,y,z)\f$ that are outside the range
310 : * of the map.
311 : *
312 : * Any point in the range of the map must be inside or on
313 : * sphere 2, and it must be outside or on sphere 1, so the inverse map
314 : * can immediately return a default-constructed
315 : * `std::optional<std::array<double, 3>>` for a point that does not
316 : * satisfy these conditions.
317 : *
318 : * Likewise, the inverse map can immediately reject any point with
319 : * \f$z < z^{-}_{\mathrm{P}2}\f$ or \f$z > z^{+}_{\mathrm{P}2}\f$.
320 : *
321 : * Finally, for \f$z^{+}_{\mathrm{P}2}\neq z^{+}_{\mathrm{P}1}\f$,
322 : * consider the circle \f$S^{+}_1\f$
323 : * defining the intersection of sphere 1
324 : * and the plane \f$z = z^{+}_{\mathrm{P}1}\f$; this circle has radius
325 : * \f$r_1 = R_1 \sin\theta_{1\mathrm{min}}\f$. Similarly, the circle
326 : * \f$S^{+}_2\f$ defining the intersection of sphere 2 and the plane \f$z =
327 : * z^{+}_{\mathrm{P}2}\f$ has radius \f$r_2 = R_2
328 : * \sin\theta_{2\mathrm{min}}\f$. Now consider the cone that passes
329 : * through these two circles. A point in the range of the map must be outside
330 : * (where "outside" means farther from the \f$z\f$ axis) or on this cone.
331 : * The cone can be defined parametrically as
332 : *
333 : * \f{align}
334 : * x_c &= C_1^x + \tilde{\lambda}(C_2^x-C_1^x) +
335 : * \cos\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\
336 : * y_c &= C_1^y + \tilde{\lambda}(C_2^y-C_1^y),+
337 : * \sin\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\
338 : * z_c &= z^{+}_{\mathrm{P}1} +
339 : * \tilde{\lambda}(z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1}),
340 : * \f}
341 : *
342 : * where \f$(x_c,y_c,z_c)\f$ is a point on the cone, and the two
343 : * parameters defining a point on the cone are the angle \f$\varphi\f$
344 : * around the cone and the parameter \f$\tilde{\lambda}\f$, which is
345 : * defined to be zero on \f$S^{+}_1\f$ and unity on \f$S^{+}_2\f$.
346 : *
347 : * Given an arbitrary point \f$(x, y, z)\f$, we can determine whether
348 : * or not that point is inside the cone as follows. First determine
349 : *
350 : * \f{align}
351 : * \tilde{\lambda} &= \frac{z - z^{+}_{\mathrm{P}1}}
352 : * {z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1}}, \\
353 : * \tilde{x} &= x - C_1^x - \tilde{\lambda} (C_2^x-C_1^x),\\
354 : * \tilde{y} &= y - C_1^y - \tilde{\lambda} (C_2^y-C_1^y).\\
355 : * \f}
356 : *
357 : * Then the condition for the point to be outside or on the cone is
358 : * \f{align}
359 : * \sqrt{\tilde{x}^2+\tilde{y}^2} \ge r_1 + (r_2-r_1)\tilde{\lambda}.
360 : * \f}
361 : *
362 : * The inverse map therefore rejects any points that do
363 : * not satisfy this criterion. The cone criterion makes sense only
364 : * for points with \f$z\geq z^{+}_{\mathrm{P}1}\f$.
365 : *
366 : * For \f$z^{-}_{\mathrm{P}2} \neq z^{-}_{\mathrm{P}1}\f$,
367 : * a similar cone can be constructed for the southern hemisphere. That
368 : * cone passes through
369 : * the circle \f$S^{-}_1\f$
370 : * defining the intersection of sphere 1
371 : * and the plane \f$z = z^{-}_{\mathrm{P}1}\f$ and the circle
372 : * \f$S^{-}_2\f$ defining the intersection of sphere 2 and the plane \f$z =
373 : * z^{-}_{\mathrm{P}2}\f$. The inverse map rejects any point that is inside
374 : * that cone as well, provided that the point has
375 : * \f$z\leq z^{-}_{\mathrm{P}1}\f$. For points with
376 : * \f$z > z^{-}_{\mathrm{P}1}\f$ checking the cone criterion
377 : * does not make sense.
378 : *
379 : * ## jacobian
380 : *
381 : * From Eqs. (\f$\ref{eq:deftheta1alt}\f$) and (\f$\ref{eq:deftheta2alt}\f$)
382 : * we see that \f$\theta_1\f$ and \f$\theta_2\f$ depend on \f$\bar{z}\f$ and
383 : * are independent of \f$\bar{x}\f$ and \f$\bar{y}\f$, and that
384 : * \f{align}
385 : * \frac{\partial (R_1\cos\theta_1)}{\partial\bar{z}}
386 : * &= \frac{1}{2}(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}),
387 : * \label{eq:dcostheta1} \\
388 : * \frac{\partial (R_2\cos\theta_2)}{\partial\bar{z}}
389 : * &= \frac{1}{2}(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}),
390 : * \label{eq:dcostheta2} \\
391 : * \frac{\partial (R_1\sin\theta_1)}{\partial\bar{z}}
392 : * &= -\frac{1}{2}\cot\theta_1 (z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}),
393 : * \label{eq:dsintheta1} \\
394 : * \frac{\partial (R_2\sin\theta_2)}{\partial\bar{z}}
395 : * &= -\frac{1}{2}\cot\theta_2(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}),
396 : * \label{eq:dsintheta2}
397 : * \f}
398 : *
399 : * Also, from Eqs. (\f$\ref{eq:defphi}\f$) and (\f$\ref{eq:defrhobar}\f$)
400 : * we have
401 : * \f{align}
402 : * \frac{\partial\cos\phi}{\partial\bar{x}}
403 : * &= \frac{\bar{y}^2}{\bar{R}^3\bar{\rho}^3},
404 : * \label{eq:dcosphidxbar} \\
405 : * \frac{\partial\cos\phi}{\partial\bar{y}}
406 : * &= -\frac{\bar{x}\bar{y}}{\bar{R}^3\bar{\rho}^3},
407 : * \label{eq:dcosphidybar} \\
408 : * \frac{\partial\sin\phi}{\partial\bar{x}}
409 : * &= -\frac{\bar{x}\bar{y}}{\bar{R}^3\bar{\rho}^3},
410 : * \label{eq:dsinphidxbar} \\
411 : * \frac{\partial\sin\phi}{\partial\bar{y}}
412 : * &= \frac{\bar{x}^2}{\bar{R}^3\bar{\rho}^3},
413 : * \label{eq:dsinphidybar}
414 : * \f}
415 : * and we know that \f$\phi\f$ is independent of \f$\bar{z}\f$.
416 : *
417 : * Finally, from Eqs. (\f$\ref{eq:defrhobar}\f$) and
418 : * (\f$\ref{eq:lambdafromrhobar}\f$) we have
419 : *
420 : * \f{align}
421 : * \frac{\partial\lambda}{\partial\bar{x}}
422 : * &= \frac{\bar{x}}{\bar{R}^2\bar{\rho}},
423 : * \label{eq:dlambdadxbar} \\
424 : * \frac{\partial\lambda}{\partial\bar{y}}
425 : * &= \frac{\bar{y}}{\bar{R}^2\bar{\rho}},
426 : * \label{eq:dlambdadybar}
427 : * \f}
428 : * with no dependence on \f$\bar{z}\f$.
429 : *
430 : * Putting these results together yields
431 : * \f{align}
432 : * \frac{\partial x^0}{\partial \bar{x}} &=
433 : * \frac{\bar{y}^2}{\bar{\rho}^3\bar{R}^3}R_1\sin\theta_1 +
434 : * (R_2\sin\theta_2-R_1\sin\theta_1)
435 : * \frac{\lambda \bar{R}^2\bar\rho^2+\bar{x}^2}{\bar\rho^3\bar{R}^3}
436 : * + \frac{\bar{x}}{\bar\rho\bar{R}^2}(C_2^x-C_1^x),\\
437 : * \frac{\partial x^0}{\partial \bar{y}} &=
438 : * \frac{\bar{x}\bar{y}}{\bar{\rho}^3\bar{R}^3}
439 : * (R_2\sin\theta_2-2 R_1\sin\theta_1)
440 : * + \frac{\bar{y}}{\bar\rho\bar{R}^2}(C_2^x-C_1^x),\\
441 : * \frac{\partial x^0}{\partial \bar{z}} &=
442 : * -\frac{1}{2}\frac{\bar{x}}{\bar\rho\bar{R}}\left[
443 : * \cot\theta_1(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})+
444 : * \cot\theta_2\lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})\right],\\
445 : * \frac{\partial x^1}{\partial \bar{x}} &=
446 : * \frac{\bar{x}\bar{y}}{\bar{\rho}^3\bar{R}^3}
447 : * (R_2\sin\theta_2-2 R_1\sin\theta_1)
448 : * + \frac{\bar{x}}{\bar\rho\bar{R}^2}(C_2^y-C_1^y),\\
449 : * \frac{\partial x^1}{\partial \bar{y}} &=
450 : * \frac{\bar{x}^2}{\bar{\rho}^3\bar{R}^3}R_1\sin\theta_1 +
451 : * (R_2\sin\theta_2-R_1\sin\theta_1)
452 : * \frac{\lambda \bar{R}^2\bar\rho^2+\bar{y}^2}{\bar\rho^3\bar{R}^3}
453 : * + \frac{\bar{y}}{\bar\rho\bar{R}^2}(C_2^y-C_1^y),\\
454 : * \frac{\partial x^1}{\partial \bar{z}} &=
455 : * -\frac{1}{2}\frac{\bar{y}}{\bar\rho\bar{R}}\left[
456 : * \cot\theta_1(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})+
457 : * \cot\theta_2\lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})\right],\\
458 : * \frac{\partial x^2}{\partial \bar{x}} &=
459 : * \frac{\bar{x}}{\bar\rho\bar{R}^2}\left(
460 : * C_2^z-C_1^z + R_2\cos\theta_2-R_1\cos\theta_1\right),\\
461 : * \frac{\partial x^2}{\partial \bar{y}} &=
462 : * \frac{\bar{y}}{\bar\rho\bar{R}^2}\left(
463 : * C_2^z-C_1^z + R_2\cos\theta_2-R_1\cos\theta_1\right),\\
464 : * \frac{\partial x^2}{\partial \bar{z}} &=
465 : * \frac{1}{2}(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})+
466 : * \frac{1}{2}\lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}).
467 : * \f}
468 : *
469 : * ## inv_jacobian
470 : *
471 : * The inverse Jacobian is computed by numerically inverting the
472 : * Jacobian.
473 : *
474 : * ## Restrictions on map parameters
475 : *
476 : * We demand that Sphere 1 is fully contained inside Sphere 2, and
477 : * that the two spheres have at least some small separation between
478 : * them. In particular, we demand that
479 : * \f{align}
480 : * 0.98 R_2 &\geq R_1 + |C_1-C_2|, \label{eq:spherecontained}
481 : * \f}
482 : * where 0.98 is a safety factor. It is possible to construct a valid
483 : * map without this assumption, but the assumption simplifies the
484 : * code, and the expected use cases obey this restriction.
485 : *
486 : * We also demand that \f$R_1 \geq 0.08 R_2\f$. Again, this assumption
487 : * is made for accuracy purposes and might be relaxed.
488 : *
489 : * ### Invertibility condition
490 : *
491 : * Consider the line segment \f$L^+_1\f$ that connects a point on the
492 : * circle \f$S^+_1\f$ (the circle formed by the intersection of sphere 1
493 : * and the plane \f$z=z^+_{\mathrm{P}1}\f$) with the center of the
494 : * circle \f$S^+_1\f$. Consider another line segment \f$L^+_2\f$ that
495 : * connects the same point on the circle \f$S^+_1\f$ with the
496 : * corresponding point on the circle \f$S^+_2\f$ (the circle formed by
497 : * the intersection of sphere 2 and the plane
498 : * \f$z=z^+_{\mathrm{P}2}\f$). Now consider the angle between \f$L^+_1\f$
499 : * and \f$L^+_2\f$, as measured from the interior of sphere 1, and Let
500 : * \f$\alpha^+\f$ be the minimum value of this angle over the circle.
501 : * \f$\alpha^+\f$ is shown in the figure above. If
502 : * \f$\alpha^+ < \theta_{1 \mathrm{min}}\f$, then the line segment \f$L^+_2\f$
503 : * twice intersects the unmapped portion of sphere 1 near the north pole,
504 : * so the map is ill-defined.
505 : * Similarly, if \f$\alpha^+ < \theta_{2 \mathrm{min}}\f$,
506 : * then the line segment \f$L^+_2\f$ twice intersects the mapped portion of
507 : * sphere 2 near the north pole, and again the map is poorly defined.
508 : * Therefore we demand that the map parameters satisfy
509 : * - \f$\alpha^+ > 1.1 \theta_{1 \mathrm{min}}\f$
510 : * - \f$\alpha^+ > 1.1 \theta_{2 \mathrm{min}}\f$
511 : *
512 : * where 1.1 is a safety factor.
513 : *
514 : * Similarly, one can define an angle \f$\alpha^-\f$ for the region
515 : * near the south pole, and we require similar restrictions on that angle.
516 : *
517 : * ### Restrictions on z-planes
518 : *
519 : * We also demand that either
520 : * \f$z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\f$
521 : * or that \f$z^+_{\mathrm{P}1} <= z^+_{\mathrm{P}2} -0.03 R_2\f$.
522 : * Similarly, we demand that either \f$z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\f$
523 : * or \f$z^-_{\mathrm{P}1} >= z^-_{\mathrm{P}2} + 0.03 R_2\f$.
524 : * These restrictions follow expected use cases and avoid extreme distortions.
525 : *
526 : * ### Restrictions for unequal z planes
527 : * For \f$z^+_{\mathrm{P}1} \neq z^+_{\mathrm{P}2}\f$ and
528 : * \f$z^-_{\mathrm{P}1} \neq z^-_{\mathrm{P}2}\f$, we assume the following
529 : * restrictions on other parameters:
530 : *
531 : * We prohibit a tiny Sphere 1 near the edge of Sphere 2 by demanding that
532 : * \f{align}
533 : * C^z_1 - R_1 &\leq C^z_2 + R_2/5,\\
534 : * C^z_1 + R_1 &\geq C^z_2 - R_2/5.
535 : * \f}
536 : * We also demand that the polar axis of Sphere 2 intersects Sphere 1
537 : * somewhere:
538 : * \f{align}
539 : * \sqrt{(C^x_1-C^x_2)^2 + (C^y_1-C^y_2)^2} &\leq R_1.
540 : * \f}
541 : * and we demand that Sphere 1 is not too close to the edge of Sphere 2
542 : * in the \f$x\f$ or \f$y\f$ directions:
543 : * \f{align}
544 : * \sqrt{(C^y_1-C^y_2)^2 + (C^y_1-C^y_2)^2} &\leq \mathrm{max}(0,0.95 R_2-R_1),
545 : * \f}
546 : * where the max avoids problems when \f$0.95 R_2-R_1\f$ is negative
547 : * (which, if it occurs, means that the \f$x\f$ and \f$y\f$ centers of the
548 : * two spheres are equal).
549 : *
550 : * We require that the z planes in the above figures lie above/below
551 : * the centers of the corresponding spheres and are not too close to
552 : * the centers or edges of those spheres; specificially, we demand
553 : * that
554 : * \f{align}
555 : * \label{eq:theta_1_min_res}
556 : * 0.15\pi &< \theta_{1 \mathrm{min}} < 0.4\pi \\
557 : * \label{eq:theta_1_max_res}
558 : * 0.6\pi &< \theta_{1 \mathrm{max}} < 0.85\pi \\
559 : * \label{eq:theta_2_min_res}
560 : * 0.15\pi &< \theta_{2 \mathrm{min}} < 0.4\pi \\
561 : * \label{eq:theta_2_max_res}
562 : * 0.6\pi &< \theta_{2 \mathrm{max}} < 0.85\pi .
563 : * \f}
564 : *
565 : * Here the numerical values are safety factors.
566 : * These restrictions are not strictly necessary but are made for simplicity.
567 : * Increasing the range will make the maps less accurate because the domain
568 : * is more distorted. These parameters
569 : * can be changed provided the unit tests are changed to test the
570 : * appropriate parameter ranges.
571 : *
572 : * ### Restrictions for equal z planes
573 : *
574 : * If \f$z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\f$ or
575 : * \f$z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\f$ we demand that
576 : * \f$C_1^x=C_2^x\f$ and \f$C_1^y=C_2^y\f$, which simplifies the cases
577 : * we need to test and agrees with our expected use cases.
578 : * We also demand
579 : * \f{align}
580 : * z^+_{\mathrm{P}2} &\geq z^-_{\mathrm{P}2} + 0.18 R_2
581 : * \f}
582 : * This condition is necessary because for unequal z planes,
583 : * \f$\theta_{2 \mathrm{min}}\f$ and
584 : * \f$\theta_{2 \mathrm{max}}\f$ are no longer required
585 : * to be on opposite sides of the equator of sphere 2 (see the paragraph below).
586 : * Note that for unequal z planes \f$\theta_{1 \mathrm{min}}\f$ and
587 : * \f$\theta_{1 \mathrm{max}}\f$ are no longer required
588 : * to be on opposite sides of the equator of sphere 1, but the conditions
589 : * in the paragraph below guarantee that
590 : * \f$z^+_{\mathrm{P}1} \geq z^-_{\mathrm{P}1}\f$.
591 : *
592 : * Unlike the case with unequal z planes, we no longer require that the
593 : * z planes in the above figures lie above/below
594 : * the centers of the corresponding spheres, but we still require that
595 : * the z planes are not too close to the edges of those spheres.
596 : * The restrictions are the same as
597 : * Eqs. (\f$\ref{eq:theta_1_min_res}\f$--\f$\ref{eq:theta_2_max_res}\f$)
598 : * except for the following changes:
599 : * If \f$z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\f$,
600 : * then we replace Eq. (\f$\ref{eq:theta_1_min_res}\f$) with
601 : * \f{align}
602 : * \label{eq:equal_plus_theta_1_min_res}
603 : * 0.15\pi &< \theta_{1 \mathrm{min}} < 0.59\pi,
604 : * \f}
605 : * and furthermore, if \f$z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\f$ and
606 : * \f$\theta_{1 \mathrm{min}} > 0.4\pi\f$ we replace
607 : * Eqs. (\f$\ref{eq:theta_1_max_res}\f$--\f$\ref{eq:theta_2_min_res}\f$)
608 : * with
609 : * \f{align}
610 : * \label{eq:equal_plus_high_theta_1_max_res}
611 : * 0.7\pi &< \theta_{1 \mathrm{max}} < 0.85\pi \\
612 : * \label{eq:equal_plus_high_theta_2_min_res}
613 : * 0.25\pi &< \theta_{2 \mathrm{min}} < 0.75\pi,
614 : * \f}
615 : * but if \f$z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\f$ and
616 : * \f$\theta_{1 \mathrm{min}} \leq 0.4\pi\f$ we replace
617 : * Eq. (\f$\ref{eq:theta_2_min_res}\f$)
618 : * with
619 : * \f{align}
620 : * \label{eq:equal_plus_low_theta_2_min_res}
621 : * 0.15\pi &< \theta_{2 \mathrm{min}} < 0.75\pi.
622 : * \f}
623 : *
624 : * Similarly, if \f$z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\f$ we replace
625 : * (\f$\ref{eq:theta_1_max_res}\f$) with
626 : * \f{align}
627 : * \label{eq:equal_minus_theta_1_max_res}
628 : * 0.41\pi &< \theta_{1 \mathrm{max}} < 0.85\pi,
629 : * \f}
630 : * and furthermore, if \f$z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\f$ and
631 : * \f$\theta_{1 \mathrm{max}} < 0.6\pi\f$ we replace
632 : * Eqs. (\f$\ref{eq:theta_1_min_res}\f$) and (\f$\ref{eq:theta_2_max_res}\f$)
633 : * with
634 : * \f{align}
635 : * \label{eq:equal_minus_high_theta_1_min_res}
636 : * 0.15\pi &< \theta_{1 \mathrm{min}} < 0.3\pi \\
637 : * \label{eq:equal_minus_high_theta_2_max_res}
638 : * 0.25\pi &< \theta_{2 \mathrm{max}} < 0.75\pi,
639 : * \f}
640 : * but if \f$z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\f$ and
641 : * \f$\theta_{1 \mathrm{max}} \geq 0.6\pi\f$ we replace
642 : * Eq. (\f$\ref{eq:theta_2_max_res}\f$)
643 : * with
644 : * \f{align}
645 : * \label{eq:equal_minus_low_theta_2_max_res}
646 : * 0.25\pi &< \theta_{2 \mathrm{max}} < 0.85\pi .
647 : * \f}
648 : */
649 1 : class UniformCylindricalSide {
650 : public:
651 0 : static constexpr size_t dim = 3;
652 0 : UniformCylindricalSide(const std::array<double, 3>& center_one,
653 : const std::array<double, 3>& center_two,
654 : double radius_one, double radius_two,
655 : double z_plane_plus_one, double z_plane_minus_one,
656 : double z_plane_plus_two, double z_plane_minus_two);
657 0 : UniformCylindricalSide() = default;
658 0 : ~UniformCylindricalSide() = default;
659 0 : UniformCylindricalSide(UniformCylindricalSide&&) = default;
660 0 : UniformCylindricalSide(const UniformCylindricalSide&) = default;
661 0 : UniformCylindricalSide& operator=(const UniformCylindricalSide&) = default;
662 0 : UniformCylindricalSide& operator=(UniformCylindricalSide&&) = default;
663 :
664 : template <typename T>
665 0 : std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
666 :
667 : /// The inverse function is only callable with doubles because the inverse
668 : /// might fail if called for a point out of range, and it is unclear
669 : /// what should happen if the inverse were to succeed for some points in a
670 : /// DataVector but fail for other points.
671 1 : std::optional<std::array<double, 3>> inverse(
672 : const std::array<double, 3>& target_coords) const;
673 :
674 : template <typename T>
675 0 : tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
676 : const std::array<T, 3>& source_coords) const;
677 :
678 : template <typename T>
679 0 : tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
680 : const std::array<T, 3>& source_coords) const;
681 :
682 : // clang-tidy: google runtime references
683 0 : void pup(PUP::er& p); // NOLINT
684 :
685 0 : static bool is_identity() { return false; }
686 :
687 0 : static constexpr bool supports_hessian{false};
688 :
689 : private:
690 0 : friend bool operator==(const UniformCylindricalSide& lhs,
691 : const UniformCylindricalSide& rhs);
692 0 : std::array<double, 3> center_one_{};
693 0 : std::array<double, 3> center_two_{};
694 0 : double radius_one_{std::numeric_limits<double>::signaling_NaN()};
695 0 : double radius_two_{std::numeric_limits<double>::signaling_NaN()};
696 0 : double z_plane_plus_one_{std::numeric_limits<double>::signaling_NaN()};
697 0 : double z_plane_minus_one_{std::numeric_limits<double>::signaling_NaN()};
698 0 : double z_plane_plus_two_{std::numeric_limits<double>::signaling_NaN()};
699 0 : double z_plane_minus_two_{std::numeric_limits<double>::signaling_NaN()};
700 : };
701 :
702 0 : bool operator!=(const UniformCylindricalSide& lhs,
703 : const UniformCylindricalSide& rhs);
704 : } // namespace domain::CoordinateMaps
|
