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