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