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