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