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