Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines the class UniformCylindricalEndcap.
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 cylinder to a volume that connects
29 : * portions of two spherical surfaces.
30 : *
31 : * \image html UniformCylEndcap.svg "A 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$. Let sphere 1 be intersected by a
35 : * plane normal to the \f$z\f$ axis and located at \f$z = z_{\mathrm{P}1}\f$,
36 : * and let sphere 2 be intersected by a plane normal to the \f$z\f$ axis and
37 : * located at \f$z = z_{\mathrm{P}2}\f$.
38 : *
39 : * UniformCylindricalEndcap maps a 3D unit right cylinder (with coordinates
40 : * \f$(\bar{x},\bar{y},\bar{z})\f$ such that \f$-1\leq\bar{z}\leq 1\f$
41 : * and \f$\bar{x}^2+\bar{y}^2 \leq 1\f$) to the shaded area
42 : * in the figure above (with coordinates \f$(x,y,z)\f$). The "bottom"
43 : * of the cylinder \f$\bar{z}=-1\f$ is mapped to the portion of sphere
44 : * 1 that has \f$z \geq z_{\mathrm{P}1}\f$, and on this portion of the
45 : * sphere the angular coordinate \f$\theta_1 = \acos((z-C_1^2)/R_1)\f$
46 : * is uniform in \f$\bar{\rho} = \sqrt{\bar{x}^2+\bar{y}^2}\f$ and the angular
47 : * coordinate \f$\phi_1 = \atan((y-C_1^1)/(x-C_1^0))\f$ is the same as
48 : * \f$\phi = \atan(\bar{y}/\bar{x})\f$.
49 : * Likewise, the "top" of the cylinder
50 : * \f$\bar{z}=+1\f$ is mapped to the portion of sphere 2 that has
51 : * \f$z \geq z_{\mathrm{P}2}\f$, and on this portion of the
52 : * sphere the angular coordinate \f$\theta_2 = \acos((z-C_2^2)/R_2)\f$
53 : * is uniform in \f$\bar\rho\f$ and the angular
54 : * coordinate \f$\phi_2 = \atan((y-C_2^1)/(x-C_2^0))\f$ is the same as
55 : * \f$\phi\f$.
56 : *
57 : * UniformCylindricalEndcap is different from CylindricalEndcap
58 : * because of the distribution of points on the spheres, and because
59 : * for UniformCylindricalEndcap the mapped portion of both Sphere 1
60 : * and Sphere 2 are bounded by planes of constant \f$z\f$, whereas for
61 : * CylindricalEndcap only one of the mapped portions is bounded by a
62 : * plane (except for specially chosen map parameters). Note that
63 : * UniformCylindricalEndcap can be used to construct maps that connect
64 : * an arbitrary number of nested spheres; this is not possible for
65 : * CylindricalEndcap for more than 3 nested spheres because of this
66 : * asymmetry between CylindricalEndcap's two spherical surfaces.
67 : *
68 : * UniformCylindricalEndcap is intended to be composed with `Wedge<2>` maps to
69 : * construct a portion of a cylindrical domain for a binary system.
70 : *
71 : * UniformCylindricalEndcap can be used to construct a domain that is similar
72 : * to, but not identical to, the one described briefly in the Appendix of
73 : * \cite Buchman:2012dw.
74 : * UniformCylindricalEndcap is used to construct the Blocks analogous to
75 : * those labeled 'CA wedge', 'EA wedge', 'CB wedge', 'EE wedge',
76 : * and 'EB wedge' in Figure 20 of that paper.
77 : *
78 : * UniformCylindricalEndcap provides the following functions:
79 : *
80 : * ### operator()
81 : *
82 : * `operator()` maps \f$(\bar{x},\bar{y},\bar{z})\f$ to \f$(x,y,z)\f$
83 : * according to
84 : *
85 : * \f{align}
86 : * x^0 &= C_1^0+\lambda(C_2^0-C_1^0) +
87 : * \cos\phi\left(R_1\sin\theta_1 +
88 : * \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x0} \\
89 : * x^1 &= C_1^1+\lambda(C_2^1-C_1^1) +
90 : * \sin\phi\left(R_1\sin\theta_1 +
91 : * \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x1} \\
92 : * x^2 &= C_1^2+\lambda(C_2^2-C_1^2) +
93 : * R_1\cos\theta_1 +
94 : * \lambda(R_2\cos\theta_2-R_1\cos\theta_1) \label{eq:x2}.
95 : * \f}
96 : *
97 : * Here
98 : * \f{align}
99 : * \lambda &= \frac{\bar{z}+1}{2},\label{eq:lambdafromzbar}\\
100 : * \theta_1 &= \bar{\rho} \theta_{1 \mathrm{max}},\label{eq:deftheta1}\\
101 : * \theta_2 &= \bar{\rho} \theta_{2 \mathrm{max}},\label{eq:deftheta2}\\
102 : * \phi &= \atan(\bar{y}/\bar{x})\label{eq:defphi},
103 : * \f}
104 : * where \f$\theta_{1 \mathrm{max}}\f$ and \f$\theta_{2 \mathrm{max}}\f$
105 : * are defined by
106 : * \f{align}
107 : * \cos(\theta_{1\mathrm{max}}) &= (z_{\mathrm{P}1}-C_1^2)/R_1,\\
108 : * \cos(\theta_{2\mathrm{max}}) &= (z_{\mathrm{P}2}-C_2^2)/R_2,
109 : * \f}
110 : * and
111 : * \f{align}
112 : * \bar{\rho} &= \sqrt{\bar{x}^2+\bar{y}^2}/\bar{R} \label{eq:defrhobar},
113 : * \f}
114 : * where \f$\bar{R}\f$ is the radius of the cylinder in barred
115 : * coordinates, which is always unity.
116 : *
117 : * ### inverse
118 : *
119 : * Given \f$(x,y,z)\f$ we want to find \f$(\bar{x},\bar{y},\bar{z})\f$.
120 : * From Eq. (\f$\ref{eq:x2}\f$) we can write \f$\lambda\f$ as a function
121 : * of \f$\bar\rho\f$:
122 : *
123 : * \f{align}
124 : * \lambda &= \frac{x^2 - C_1^2 - R_1\cos\theta_1}
125 : * {C_2^2-C_1^2 + R_2\cos\theta_2-R_1\cos\theta_1}
126 : * \label{eq:lambda_from_rho}.
127 : * \f}
128 : *
129 : * Then by eliminating \f$\phi\f$ from Eqs. (\f$\ref{eq:x0}\f$) and
130 : * (\f$\ref{eq:x1}\f$) we find that \f$\bar{\rho}\f$ is the solution
131 : * of \f$Q(\bar{\rho})=0\f$, where
132 : *
133 : * \f{align}
134 : * Q(\bar{\rho}) &= \left(x^0-C_1^0-\lambda(C_2^0-C_1^0)\right)^2+
135 : * \left(x^1-C_1^1-\lambda(C_2^1-C_1^1)\right)^2-
136 : * \left((1-\lambda)R_1\sin\theta_1 +
137 : * \lambda R_2\sin\theta_2\right)^2.\label{eq:defQ}
138 : * \f}
139 : * Here \f$\lambda\f$, \f$\theta_1\f$, and \f$\theta_2\f$ are functions
140 : * of \f$\bar{\rho}\f$ through Eqs. (\f$\ref{eq:lambda_from_rho}\f$),
141 : * (\f$\ref{eq:deftheta1}\f$), and (\f$\ref{eq:deftheta2}\f$).
142 : *
143 : * We solve \f$Q(\bar{\rho})=0\f$ numerically; it is a one-dimensional
144 : * root-finding problem.
145 : *
146 : * Once we have determined \f$\bar{\rho}\f$, we then obtain \f$\lambda\f$
147 : * from Eq. (\f$\ref{eq:lambda_from_rho}\f$), and we obtain \f$\phi\f$ from
148 : *
149 : * \f{align}
150 : * \tan\phi &=
151 : * \frac{x^1-C_1^1-\lambda(C_2^1-C_1^1)}{x^0-C_1^0-\lambda(C_2^0-C_1^0)}.
152 : * \f}
153 : *
154 : * Then \f$\bar{z}\f$ is obtained from Eq. (\f$\ref{eq:lambdafromzbar}\f$)
155 : * and \f$\bar{x}\f$ and \f$\bar{y}\f$ are obtained from
156 : *
157 : * \f{align}
158 : * \bar{x} &= \bar{\rho}\bar{R}\cos\phi,\\
159 : * \bar{y} &= \bar{\rho}\bar{R}\sin\phi.
160 : * \f}
161 : *
162 : * #### Considerations when root-finding.
163 : *
164 : * We solve \f$Q(\bar{\rho})=0\f$ numerically for \f$\bar{\rho}\f$,
165 : * where \f$Q(\bar{\rho})\f$ is given by Eq. (\f$\ref{eq:defQ}\f$).
166 : *
167 : * ##### min/max values of \f$\bar{\rho}\f$:
168 : *
169 : * Note that the root we care about must have
170 : * \f$0\leq\lambda\leq 1\f$; therefore from Eq. (\f$\ref{eq:lambda_from_rho}\f$)
171 : * we have
172 : *
173 : * \f{align}
174 : * \bar{\rho}_{\mathrm{min}} &=
175 : * \left\{\begin{array}{ll}
176 : * 0 & \text{for } x^2-C_1^2 \geq R_1, \\
177 : * \displaystyle \frac{1}{\theta_{1 \mathrm{max}}}
178 : * \cos^{-1}\left(\frac{x^2-C_1^2}{R_1}\right) & \text{otherwise}
179 : * \end{array}\right.\label{eq:rhobarmin}\\
180 : * \bar{\rho}_{\mathrm{max}} &=
181 : * \left\{\begin{array}{ll}
182 : * 1 & \text{for } x^2-C_2^2 \leq R_2\cos\theta_{2 \mathrm{max}}, \\
183 : * \displaystyle \frac{1}{\theta_{2 \mathrm{max}}}
184 : * \cos^{-1}\left(\frac{x^2-C_2^2}{R_2}\right) & \text{otherwise}
185 : * \end{array}\right.\label{eq:rhobarmax}
186 : * \f}
187 : *
188 : * so we look for a root only between \f$\bar{\rho}_{\mathrm{min}}\f$
189 : * and \f$\bar{\rho}_{\mathrm{max}}\f$.
190 : *
191 : * ##### Roots within roundoff of endpoints:
192 : *
193 : * Sometimes a root is within roundoff of \f$\bar{\rho}_{\mathrm{min}}\f$
194 : * or \f$\bar{\rho}_{\mathrm{max}}\f$. This tends to happen at points on the
195 : * boundary of the mapped region. In this case, the root might
196 : * not be bracketed by
197 : * \f$[\bar{\rho}_{\mathrm{min}},\bar{\rho}_{\mathrm{max}}]\f$ if the root
198 : * is slightly outside that interval. If we find that
199 : * \f$Q(\bar{\rho}_{\mathrm{min}})\f$ is near zero but has the wrong sign,
200 : * then we slightly expand the interval as follows:
201 : *
202 : * \f{align}
203 : * \bar{\rho}_{\mathrm{min}} \to \bar{\rho}_{\mathrm{min}}
204 : * - 2 \frac{Q(\bar{\rho}_{\mathrm{min}})}{Q'(\bar{\rho}_{\mathrm{min}})},
205 : * \f}
206 : *
207 : * where \f$Q'(\bar{\rho}_{\mathrm{min}})\f$ is the derivative of the function
208 : * in Eq. (\f$\ref{eq:defQ}\f$). Note that without the factor of 2, this is
209 : * a Newton-Raphson step; the factor of 2 is there to overcompensate so that
210 : * the new \f$\bar{\rho}_{\mathrm{min}}\f$ brackets the root.
211 : *
212 : * Similarly, if it is found that \f$Q(\bar{\rho}_{\mathrm{max}})\f$
213 : * is near zero but has the wrong sign so that the root is not
214 : * bracketed, then the same formula is used to expand the interval
215 : * near \f$\bar{\rho} = \bar{\rho}_{\mathrm{max}}\f$ to bracket the
216 : * root.
217 : *
218 : * Note that by differentiating Eqs. (\f$\ref{eq:defQ}\f$) and
219 : * (\f$\ref{eq:lambda_from_rho}\f$), one obtains
220 : *
221 : * \f{align}
222 : * Q'(\bar{\rho}) =& -2 \frac{d\lambda}{d\bar{\rho}}\left[
223 : * \left(x^0-C_1^0-\lambda(C_2^0-C_1^0)\right)(C_2^0-C_1^0)+
224 : * \left(x^1-C_1^1-\lambda(C_2^1-C_1^1)\right)(C_2^1-C_1^1)
225 : * \right]\nonumber \\
226 : * &
227 : * -2 \left((1-\lambda)R_1\sin\theta_1 +
228 : * \lambda R_2\sin\theta_2\right)
229 : * \left[
230 : * \frac{d\lambda}{d\bar{\rho}} (R_2\sin\theta_2-R_1\sin\theta_1)
231 : * +(1-\lambda)R_1\theta_{1 \mathrm{max}}\cos\theta_1
232 : * +\lambda R_2\theta_{2 \mathrm{max}}\cos\theta_2
233 : * \right], \label{eq:defQderiv}
234 : * \f}
235 : *
236 : * where
237 : * \f{align}
238 : * \frac{d\lambda}{d\bar{\rho}} &=
239 : * \frac{(1-\lambda)R_1\theta_{1 \mathrm{max}}\sin\theta_1
240 : * +\lambda R_2\theta_{2 \mathrm{max}}\sin\theta_2}
241 : * {C_2^2-C_1^2 + R_2\cos\theta_2-R_1\cos\theta_1}
242 : * \label{eq:dlambda_drhobar}.
243 : * \f}
244 : *
245 : * ##### Roots within roundoff of \f$\bar{\rho}=0\f$ or \f$\bar{\rho}=1\f$:
246 : *
247 : * For some points on the boundary of the mapped domain, the root will
248 : * be within roundoff of \f$\bar{\rho}=0\f$ or \f$\bar{\rho}=1\f$.
249 : * Here it does not always make sense to expand the range of the map
250 : * if the root fails (by roundoff) to be bracketed, as is done above.
251 : * Furthermore, when \f$\bar{\rho}=0\f$ is a root it turns
252 : * out that both \f$Q(\bar{\rho})=0\f$ and \f$Q'(\bar{\rho})=0\f$ for
253 : * \f$\bar{\rho}=0\f$, so some root-finders (e.g. Newton-Raphson) have
254 : * difficulty converging. Therefore the cases where the root is
255 : * within roundoff of \f$\bar{\rho}=0\f$ or \f$\bar{\rho}=1\f$ are
256 : * treated separately.
257 : *
258 : * These cases are detected by comparing terms in the first-order
259 : * power series of \f$Q(\bar{\rho})=0\f$ when expanded about
260 : * \f$\bar{\rho}=0\f$ or \f$\bar{\rho}=1\f$. When one of these cases is
261 : * recognized, the root is returned as either \f$\bar{\rho}=0\f$ or
262 : * \f$\bar{\rho}=1\f$.
263 : *
264 : * #### Quick rejection of points out of range of the map.
265 : *
266 : * It is expected that `inverse()` will often be passed points
267 : * \f$(x,y,z)\f$ that are out of the range of the map; in this case
268 : * `inverse()` returns a default-constructed
269 : * `std::optional<std::array<double, 3>>`. To avoid the difficulty and
270 : * expense of attempting to solve \f$Q(\bar{\rho})=0\f$ numerically
271 : * for such points (and then having this solution fail), it is useful
272 : * to quickly reject points \f$(x,y,z)\f$ that are outside the range
273 : * of the map.
274 : *
275 : * Any point in the range of the map must be inside or on
276 : * sphere 2, and it must be outside or on sphere 1, so the inverse map
277 : * can immediately return a default-constructed
278 : * `std::optional<std::array<double, 3>>` for a point that does not
279 : * satisfy these conditions.
280 : *
281 : * Likewise, the inverse map can immediately reject any point with
282 : * \f$z < z_{\mathrm{P}1}\f$.
283 : *
284 : * Finally, consider the circle \f$S_1\f$ defining the intersection of sphere 1
285 : * and the plane \f$z = z_{\mathrm{P}1}\f$; this circle has radius
286 : * \f$r_1 = R_1 \sin\theta_{1\mathrm{max}}\f$. Similarly, the circle
287 : * \f$S_2\f$ defining the intersection of sphere 2 and the plane \f$z =
288 : * z_{\mathrm{P}2}\f$ has radius \f$r_2 = R_2
289 : * \sin\theta_{2\mathrm{max}}\f$. Now consider the cone that passes
290 : * through these two circles. A point in the range of the map must be inside
291 : * or on this cone. The cone can be defined parametrically as
292 : *
293 : * \f{align}
294 : * x_c &= C_1^0 + \tilde{\lambda}(C_2^0-C_1^0) +
295 : * \cos\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\
296 : * y_c &= C_1^1 + \tilde{\lambda}(C_2^1-C_1^1),+
297 : * \sin\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\
298 : * z_c &= C_1^2 + R_1 \cos\theta_{1\mathrm{max}} +
299 : * \tilde{\lambda}(C_2^2 + R_2 \cos\theta_{2\mathrm{max}} -
300 : * C_1^2 - R_1 \cos\theta_{1\mathrm{max}}),
301 : * \f}
302 : *
303 : * where \f$(x_c,y_c,z_c)\f$ is a point on the cone, and the two
304 : * parameters defining a point on the cone are the angle \f$\varphi\f$
305 : * around the cone and the parameter \f$\tilde{\lambda}\f$, which is
306 : * defined to be zero on \f$S_1\f$ and unity on \f$S_2\f$.
307 : *
308 : * Given an arbitrary point \f$(x, y, z)\f$, we can determine whether
309 : * or not that point is inside the cone as follows. First determine
310 : *
311 : * \f{align}
312 : * \tilde{\lambda} &= \frac{z - C_1^2 - R_1 \cos\theta_{1\mathrm{max}}}
313 : * {C_2^2+ R_2 \cos\theta_{2\mathrm{max}}-
314 : * C_1^2- R_1 \cos\theta_{1\mathrm{max}}}, \\
315 : * \tilde{x} &= x - C_1^0 - \tilde{\lambda} (C_2^0-C_1^0),\\
316 : * \tilde{y} &= y - C_1^1 - \tilde{\lambda} (C_2^1-C_1^1).\\
317 : * \f}
318 : *
319 : * Then the condition for the point to be inside or on the cone is
320 : * \f{align}
321 : * \sqrt{\tilde{x}^2+\tilde{y}^2} \le r_1 + (r_2-r_1)\tilde{\lambda}.
322 : * \f}
323 : *
324 : * The inverse map can therefore reject any points that do
325 : * not satisfy this criterion.
326 : *
327 : * ### jacobian
328 : *
329 : * One can rewrite Eqs.(\f$\ref{eq:x0}\f$) through (\f$\ref{eq:x2}\f$) as
330 : *
331 : * \f{align}
332 : * x^0 &= \frac{1}{2}\left((1-\bar{z})C_1^0+ (1+\bar{z})C_2^0\right) +
333 : * \frac{\bar{x}}{2}\left(
334 : * (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) +
335 : * (1+\bar{z}) R_2 S(\bar{\rho},\theta_{2 \mathrm{max}})
336 : * \right), \label{eq:x0alt} \\
337 : * x^1 &= \frac{1}{2}\left((1-\bar{z})C_1^1 + (1+\bar{z})C_2^1\right) +
338 : * \frac{\bar{y}}{2}\left(
339 : * (1-\bar{z})R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) +
340 : * (1+\bar{z})R_2 S(\bar{\rho},\theta_{2 \mathrm{max}})
341 : * \right), \label{eq:x1alt} \\
342 : * x^2 &= \frac{1}{2}\left((1-\bar{z})C_1^2 + (1+\bar{z})C_2^2\right) +
343 : * \left(
344 : * (1-\bar{z})R_1 \cos\theta_1 +
345 : * (1+\bar{z})R_2 \cos\theta_2
346 : * \right), \label{eq:x2alt} \\
347 : * \f}
348 : *
349 : * where we have used Eq. (\f$\ref{eq:lambdafromzbar}\f$) to eliminate
350 : * \f$\lambda\f$ in favor of \f$\bar{z}\f$, and where we have defined the
351 : * function
352 : *
353 : * \f{align}
354 : * S(\bar{\rho},a) = \frac{\sin(\bar{\rho} a)}{\bar{\rho}}. \label{eq:Sdef}
355 : * \f}
356 : *
357 : * Note that \f$S(\bar{\rho},a)\f$ is finite as \f$\bar{\rho}\f$
358 : * approaches zero, and in the code we must take care that everything
359 : * remains well-behaved in that limit.
360 : *
361 : * Then differentiating Eqs. (\f$\ref{eq:x0alt}\f$) and (\f$\ref{eq:x1alt}\f$)
362 : * with respect to \f$\bar{x}\f$ and \f$\bar{y}\f$, taking into account the
363 : * dependence of \f$\bar{\rho}\f$ on \f$\bar{x}\f$ and \f$\bar{y}\f$ from Eq.
364 : * (\f$\ref{eq:defrhobar}\f$), we find:
365 : *
366 : * \f{align}
367 : * \frac{\partial x^0}{\partial \bar{x}} &=
368 : * \frac{1}{2}\left(
369 : * (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) +
370 : * (1+\bar{z}) R_2 S(\bar{\rho},\theta_{2 \mathrm{max}})
371 : * \right) +
372 : * \frac{\bar{x}^2}{2\bar{\rho}} \left[
373 : * (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}) +
374 : * (1+\bar{z}) R_2 S'(\bar{\rho},\theta_{2 \mathrm{max}})
375 : * \right], \\
376 : * \frac{\partial x^1}{\partial \bar{y}} &=
377 : * \frac{1}{2}\left(
378 : * (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) +
379 : * (1+\bar{z}) R_2 S(\bar{\rho},\theta_{2 \mathrm{max}})
380 : * \right) +
381 : * \frac{\bar{y}^2}{2\bar{\rho}} \left[
382 : * (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}) +
383 : * (1+\bar{z}) R_2 S'(\bar{\rho},\theta_{2 \mathrm{max}})
384 : * \right], \\
385 : * \frac{\partial x^0}{\partial \bar{y}} &=
386 : * \frac{\bar{x}\bar{y}}{2\bar{\rho}}\left[
387 : * (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}) +
388 : * (1+\bar{z}) R_2 S'(\bar{\rho},\theta_{2 \mathrm{max}})
389 : * \right], \\
390 : * \frac{\partial x^1}{\partial \bar{x}} &=
391 : * \frac{\partial x^0}{\partial \bar{y}},
392 : * \f}
393 : *
394 : * where \f$S'(\bar{\rho},a)\f$ means the derivative of \f$S(\bar{\rho},a)\f$
395 : * with respect to \f$\bar\rho\f$. Note that \f$S'(\bar{\rho},a)/\bar{\rho}\f$
396 : * approaches a constant value as \f$\bar{\rho}\f$ approaches zero.
397 : *
398 : * Differentiating Eq. (\f$\ref{eq:x2alt}\f$) with respect to
399 : * \f$\bar{x}\f$ and \f$\bar{y}\f$ we find
400 : *
401 : * \f{align}
402 : * \frac{\partial x^2}{\partial \bar{x}} &=
403 : * - \frac{\bar{x}}{2}\left[
404 : * (1-\bar{z}) R_1 \theta_{1 \mathrm{max}}
405 : * S(\bar{\rho},\theta_{1 \mathrm{max}}) +
406 : * (1+\bar{z}) R_2 \theta_{2 \mathrm{max}}
407 : * S(\bar{\rho},\theta_{2 \mathrm{max}})-
408 : * \right],\\
409 : * \frac{\partial x^2}{\partial \bar{y}} &=
410 : * - \frac{\bar{y}}{2}\left[
411 : * (1-\bar{z}) R_1 \theta_{1 \mathrm{max}}
412 : * S(\bar{\rho},\theta_{1 \mathrm{max}}) +
413 : * (1+\bar{z}) R_2 \theta_{2 \mathrm{max}}
414 : * S(\bar{\rho},\theta_{2 \mathrm{max}})-
415 : * \right].
416 : * \f}
417 : *
418 : * Differentiating Eqs. (\f$\ref{eq:x0alt}\f$) through (\f$\ref{eq:x2alt}\f$)
419 : * with respect to \f$\bar{z}\f$ yields
420 : *
421 : * \f{align}
422 : * \frac{\partial x^0}{\partial \bar{z}} &=
423 : * \frac{1}{2}\left[
424 : * C_2^0-C_1^0 +
425 : * \bar{x}\left(R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) -
426 : * R_1 S(\bar{\rho},\theta_{1 \mathrm{max}})\right)
427 : * \right],\\
428 : * \frac{\partial x^1}{\partial \bar{z}} &=
429 : * \frac{1}{2}\left[
430 : * C_2^1-C_1^1 +
431 : * \bar{y}\left(R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) -
432 : * R_1 S(\bar{\rho},\theta_{1 \mathrm{max}})\right)
433 : * \right],\\
434 : * \frac{\partial x^2}{\partial \bar{z}} &=
435 : * \frac{1}{2}\left(
436 : * C_2^2-C_1^2 + R_2\cos\theta_2 - R_1\cos\theta_1
437 : * \right).
438 : * \f}
439 : *
440 : * ### inv_jacobian
441 : *
442 : * The inverse Jacobian is computed by numerically inverting the
443 : * Jacobian.
444 : *
445 : * ### Restrictions on map parameters
446 : *
447 : * We demand that Sphere 1 is fully contained inside Sphere 2, and
448 : * that the two spheres have at least some small separation between
449 : * them. In particular, we demand that
450 : * - \f$0.98 R_2 >= R_1 + |C_1-C_2|\f$
451 : *
452 : * where 0.98 is a safety factor. It is possible to construct a valid
453 : * map without this assumption, but the assumption simplifies the
454 : * code, and the expected use cases obey this restriction.
455 : *
456 : * We also demand that \f$z_{\mathrm{P}1} <= z_{\mathrm{P}2} -0.04 R_2\f$, and
457 : * that the z planes in the above figures lie above the centers of the
458 : * corresponding spheres and are not too close to the centers or edges of
459 : * those spheres; specificially, we demand that
460 : * - \f$ 0.075\pi < \theta_{1 \mathrm{max}} < 0.45\pi\f$
461 : * - \f$ 0.075\pi < \theta_{2 \mathrm{max}} < 0.45\pi\f$
462 : *
463 : * Here 0.075 and 0.45 are safety factors. These restrictions are not
464 : * strictly necessary but are made for simplicity and to ensure the
465 : * accuracy of the inverse map (the inverse map becomes less accurate if
466 : * the map parameters are extreme).
467 : *
468 : * Consider the line segment \f$L\f$ that connects a point on the
469 : * circle \f$S_1\f$ (the circle formed by the intersection of sphere 1
470 : * and the plane \f$z=z_{\mathrm{P}1}\f$) with the center of the
471 : * circle \f$S_1\f$. Consider another line segment \f$L'\f$ that
472 : * connects the same point on the circle \f$S_1\f$ with the
473 : * corresponding point on the circle \f$S_2\f$ (the circle formed by
474 : * the intersection of sphere 2 and the plane
475 : * \f$z=z_{\mathrm{P}2}\f$). Now consider the angle between \f$L\f$
476 : * and \f$L'\f$, as measured from the interior of sphere 1, and Let
477 : * \f$\alpha\f$ be the minimum value of this angle over the circle.
478 : * \f$\alpha\f$ is shown in the figure above. If
479 : * \f$\alpha < \theta_{1 \mathrm{max}}\f$, then the line segment \f$L'\f$
480 : * intersects the mapped portion of sphere 1 twice, so the map is
481 : * multiply valued. Similarly, if \f$\alpha < \theta_{2 \mathrm{max}}\f$,
482 : * then the line segment \f$L'\f$ intersects the mapped portion of
483 : * sphere 2 twice, and again the map is multiply valued. Therefore we
484 : * demand that the map parameters are such that
485 : * - \f$\alpha > 1.1 \theta_{1 \mathrm{max}}\f$
486 : * - \f$\alpha > 1.1 \theta_{2 \mathrm{max}}\f$
487 : *
488 : * where 1.1 is a safety factor.
489 : *
490 : * The condition on \f$\alpha\f$ is guaranteed to provide an
491 : * invertible map if \f$C_1^0=C_2^0\f$ and \f$C_1^1=C_2^1\f$.
492 : * However, for \f$C_1^0 \neq C_2^0\f$ or \f$C_1^1\neq C_2^1\f$, even
493 : * if the \f$\alpha\f$ condition is satisfied, it is possible for two
494 : * lines of constant \f$(\bar{x},\bar{y})\f$ (each line has different
495 : * values of \f$(\bar{x},\bar{y})\f$) to pass through the same point
496 : * \f$(x,y,z)\f$ if those lines are not coplanar. This condition is
497 : * difficult to check analytically, so we check it numerically. We
498 : * have found empirically that if \f$Q(\bar{\rho})\f$ from
499 : * Eq. (\f$\ref{eq:defQ}\f$) has only a single root between
500 : * \f$\bar{\rho}_{\mathrm{min}}\f$ and \f$\bar{\rho}_{\mathrm{max}}\f$
501 : * for all points \f$(x,y,z)\f$ on the surface of sphere 1 with
502 : * \f$z\geq z_{\mathrm{P}1}\f$ and with \f$(x-C_1^0)/(y-C_1^1) =
503 : * (C_1^0-C_2^0)/(C_1^1-C_2^1)\f$, then the map is single-valued
504 : * everywhere. We cannot numerically check every point in this
505 : * one-parameter family of points, but we demand that this condition
506 : * is satisfied for a reasonably large number of points in this
507 : * family. This check is not very expensive since it is done only
508 : * once, in the constructor.
509 : *
510 : */
511 1 : class UniformCylindricalEndcap {
512 : public:
513 0 : static constexpr size_t dim = 3;
514 0 : UniformCylindricalEndcap(const std::array<double, 3>& center_one,
515 : const std::array<double, 3>& center_two,
516 : double radius_one, double radius_two,
517 : double z_plane_one, double z_plane_two);
518 0 : UniformCylindricalEndcap() = default;
519 0 : ~UniformCylindricalEndcap() = default;
520 0 : UniformCylindricalEndcap(UniformCylindricalEndcap&&) = default;
521 0 : UniformCylindricalEndcap(const UniformCylindricalEndcap&) = default;
522 0 : UniformCylindricalEndcap& operator=(const UniformCylindricalEndcap&) =
523 : default;
524 0 : UniformCylindricalEndcap& operator=(UniformCylindricalEndcap&&) = default;
525 :
526 : template <typename T>
527 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
528 : const std::array<T, 3>& source_coords) const;
529 :
530 : /// The inverse function is only callable with doubles because the inverse
531 : /// might fail if called for a point out of range, and it is unclear
532 : /// what should happen if the inverse were to succeed for some points in a
533 : /// DataVector but fail for other points.
534 1 : std::optional<std::array<double, 3>> inverse(
535 : const std::array<double, 3>& target_coords) const;
536 :
537 : template <typename T>
538 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
539 : const std::array<T, 3>& source_coords) const;
540 :
541 : template <typename T>
542 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
543 : const std::array<T, 3>& source_coords) const;
544 :
545 : // clang-tidy: google runtime references
546 0 : void pup(PUP::er& p); // NOLINT
547 :
548 0 : static bool is_identity() { return false; }
549 :
550 : private:
551 0 : friend bool operator==(const UniformCylindricalEndcap& lhs,
552 : const UniformCylindricalEndcap& rhs);
553 0 : std::array<double, 3> center_one_{};
554 0 : std::array<double, 3> center_two_{};
555 0 : double radius_one_{std::numeric_limits<double>::signaling_NaN()};
556 0 : double radius_two_{std::numeric_limits<double>::signaling_NaN()};
557 0 : double z_plane_one_{std::numeric_limits<double>::signaling_NaN()};
558 0 : double z_plane_two_{std::numeric_limits<double>::signaling_NaN()};
559 0 : double theta_max_one_{std::numeric_limits<double>::signaling_NaN()};
560 0 : double theta_max_two_{std::numeric_limits<double>::signaling_NaN()};
561 : };
562 :
563 0 : bool operator!=(const UniformCylindricalEndcap& lhs,
564 : const UniformCylindricalEndcap& rhs);
565 :
566 : /// Given parameters for UniformCylindricalEndcap, returns whether
567 : /// the map is invertible for target points on sphere_one.
568 : ///
569 : /// `is_uniform_cylindrical_endcap_invertible_on_sphere_one` is
570 : /// publicly visible because it is useful for unit tests that need to
571 : /// choose valid parameters to pass into the map.
572 1 : bool is_uniform_cylindrical_endcap_invertible_on_sphere_one(
573 : const std::array<double, 3>& center_one,
574 : const std::array<double, 3>& center_two, const double radius_one,
575 : const double radius_two, const double theta_max_one,
576 : const double theta_max_two);
577 :
578 : } // namespace domain::CoordinateMaps
|
