Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 : #include <memory>
10 : #include <optional>
11 : #include <ostream>
12 : #include <pup.h>
13 :
14 : #include "DataStructures/DataVector.hpp"
15 : #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
16 :
17 : namespace domain::CoordinateMaps::ShapeMapTransitionFunctions {
18 : /*!
19 : * \brief A transition function $G(r,\theta,\phi)$ with $f(r,\theta,\phi)$ that
20 : * falls off linearly from an inner surface of a wedge to an outer surface of a
21 : * wedge. Meant to be used in `domain::CoordinateMaps::Wedge` blocks.
22 : *
23 : * \details The functional form of this transition is
24 : *
25 : * \begin{equation}
26 : * G(r,\theta,\phi) = \frac{f(r,\theta,\phi)}{r}
27 : * \label{eq:transition_func}
28 : * \end{equation}
29 : *
30 : * where
31 : *
32 : * \begin{equation}
33 : * f(r, \theta, \phi) = \frac{D_{\text{out}}(r, \theta, \phi) -
34 : * r}{D_{\text{out}}(r, \theta, \phi) - D_{\text{in}}(r, \theta, \phi)},
35 : * \label{eq:linear_transition_func}
36 : * \end{equation}
37 : *
38 : * where, in general,
39 : *
40 : * \begin{equation}
41 : * D(r, \theta, \phi) = R\left(\frac{1 -
42 : * s}{\sqrt{3}\max(|\sin{\theta}\cos{\phi}|,|\sin{\theta}\sin{\phi}|,
43 : * |\cos{\theta}|)} + s\right).
44 : * \label{eq:distance}
45 : * \end{equation}
46 : *
47 : * Here, $s$ is the sphericity of the surface which goes from 0 (flat) to 1
48 : * (spherical), $R$ is the radius of the spherical surface, $\text{out}$ is the
49 : * outer surface, and $\text{in}$ is the inner surface. If the sphericity is 1,
50 : * then the surface is a sphere so $D = R$. If the sphericity is 0, then the
51 : * surface is a cube. This cube is circumscribed by a sphere of radius $R$. See
52 : * `domain::CoordinateMaps::Wedge` for more of an explanation of these boundary
53 : * surfaces and their sphericities.
54 : *
55 : * \note If \p reverse is true, then the functional form of the transition is
56 : * actually
57 : * \begin{equation}
58 : * f(r, \theta, \phi) = 1 - \frac{D_{\text{out}}(r, \theta, \phi) -
59 : * r}{D_{\text{out}}(r, \theta, \phi) - D_{\text{in}}(r, \theta, \phi)} =
60 : * \frac{r - D_{\text{in}}(r, \theta, \phi)}
61 : * {D_{\text{out}}(r, \theta, \phi) - D_{\text{in}}(r, \theta, \phi)}.
62 : * \label{eq:linear_transition_func_reverse}
63 : * \end{equation}
64 : *
65 : * \parblock
66 : *
67 : * \note Because the shape map distorts only radii and does not affect angles,
68 : * $D$ is not a function of $r$ so we have that $D(r,\theta,\phi) =
69 : * D(\theta, \phi)$.
70 : *
71 : * \endparblock
72 : *
73 : * The transition function is also defined within $D_{\text{in}}(r, \theta,
74 : * \phi)$ if the axis is `Axis::Interior`. Within $D_{\text{in}}(r, \theta,
75 : * \phi)$, the transition function is
76 : *
77 : * \begin{equation}
78 : * \label{eq:interior_transition_func}
79 : * G(r,\theta,\phi) = \frac{r^2}{D_{\text{in}}(r, \theta, \phi)^3} =
80 : * \frac{r^3}{R_0^3}
81 : * \end{equation}
82 : *
83 : * and none of the vector formulation is necessary. This is chosen to match Eq.
84 : * $\ref{eq:transition_func}$ at $D_{\text{in}}(r, \theta, \phi)$ and go to 0 at
85 : * $r=0$. The second equality is explained by one of the following bullets.
86 : *
87 : * There are several assumptions made for this mapping:
88 : *
89 : * - The coordinates $r, \theta, \phi$ are assumed to be from the center of the
90 : * inner surface, not the center of the computational domain.
91 : * - The wedges are concentric. (see the constructor) This is also enforced by
92 : * the center of the inner surface being within the outer surface.
93 : * - If the centers of the inner and outer surface are different, then the inner
94 : * sphericity must be 1 (i.e. $D_{\text{in}} = R$).
95 : * - If the axis is `Axis::Interior`, then the inner sphericity must be 1 and
96 : * \p reverse must be false.
97 : * - The $\max$ in the denominator of $\ref{eq:distance}$ can be simplified a
98 : * bit to $\max(|x|, |y|, |z|)/r$. It was written the other way in
99 : * $\ref{eq:distance}$ to emphasize that $D$ has no radial dependence.
100 : *
101 : * ### Vector formulation
102 : *
103 : * \image html WedgeTransition.jpg "Useful vectors for Wedge transition"
104 : *
105 : * It can be more convenient to work with vectors rather than just distances for
106 : * the following equations. Therefore, we define the projection center $\vec P$
107 : * which points from the center of the outer surface to the center of the inner
108 : * surface. If the centers of the inner and outer surface are the same, then
109 : * $\vec P = 0$. Then, we define $\vec x$ as the vector from the center of the
110 : * outer surface to the coordinate we are interested in. We define $\vec x_0$ as
111 : * the vector from the center of the outer surface to the point along ray $\vec
112 : * x - \vec P$ which intersects the inner surface. Similarly, $\vec x_1$
113 : * intersects the outer surface along the same ray as $\vec x - \vec P$. Any
114 : * cube has a side length of $2L$.
115 : *
116 : * In this way, we can reformulate Eq. $\ref{eq:linear_transition_func}$ as
117 : *
118 : * \begin{equation}\label{eq:linear_transition_func_vec}
119 : * f(\vec x) = \frac{|\vec x_1 - \vec P| - r}{|\vec x_1 - \vec P| - |\vec
120 : * x_0 - \vec P|}
121 : * \end{equation}
122 : *
123 : * where
124 : *
125 : * \begin{equation}
126 : * r = |\vec x - P|
127 : * \end{equation}
128 : *
129 : * #### Computing vector to inner surface
130 : *
131 : * Similar to Eq. $\ref{eq:distance}$ we can define $\vec x_0$ as
132 : *
133 : * \begin{equation}\label{eq:x_0_vector}
134 : * \vec x_0 = R_0\left( \frac{r(1-s_0)}{\sqrt{3}\max{(x_i - P_i)}} +
135 : * s\right) \hat r + \vec P
136 : * \end{equation}
137 : *
138 : * with
139 : *
140 : * \begin{equation}
141 : * \hat r = \frac{\vec x - \vec P}{|\vec x - \vec P|}.
142 : * \end{equation}
143 : *
144 : * and also
145 : *
146 : * \begin{equation}\label{eq:x_0_norm}
147 : * |\vec x_0 - \vec P| = R_0\left( \frac{r(1-s_0)}{\sqrt{3}\max{(x_i -
148 : * P_i)}} + s_0\right)
149 : * \end{equation}
150 : *
151 : * ### Computing vector to outer surface
152 : *
153 : * Then, we define
154 : *
155 : * \begin{equation}\label{eq:x_1}
156 : * \vec x_1 = \lambda (\vec x - \vec P) + \vec P
157 : * \end{equation}
158 : *
159 : * where $\lambda$ is a scalar to be determined. We can see that this is valid
160 : * because $\vec x - \vec P$ is centered on $\vec P$ and lies along the ray so
161 : * $\lambda$ will scale this vector to be the correct length.
162 : *
163 : * Computing $\lambda$ is a bit complicated because the outer surface can either
164 : * be a cube, a sphere, or something in between and the way to solve for
165 : * $\lambda$ in the cube case is different than in the sphere case. Therefore we
166 : * define $\lambda_c$ and $\lambda_s$ as the scalar factors for the cube and
167 : * sphere case, respectively. This allows us to define $\lambda$ as
168 : *
169 : * \begin{equation}\label{eq:lambda}
170 : * \lambda = (1 - s_1) \lambda_c + s_1 \lambda_s
171 : * \end{equation}
172 : *
173 : * #### Computing cube lambda
174 : *
175 : * To compute $\lambda_c$ for the cube case, we notice that (in the image)
176 : *
177 : * \begin{align}
178 : * L &= \vec x_1 \cdot \hat z \\
179 : * &= \left(\lambda_c^{+z}(\vec x - \vec P) + \vec P\right) \cdot \hat z\\
180 : * &= \lambda_c^{+z}(x_z - P_z) + P_z
181 : * \end{align}
182 : *
183 : * If we generalize this to all six unit vectors and solve for $\lambda_c^{\pm
184 : * k}$ we get
185 : *
186 : * \begin{equation}\label{eq:lambda_c_pm}
187 : * \lambda_c^{\pm k} = \frac{\pm L - \vec P \cdot \hat x_k}{(\vec x -
188 : * \vec P) \cdot \hat x_k}
189 : * \end{equation}
190 : *
191 : * for $k \in \{x, y, z\}$. This quantity $\lambda_c^{\pm k}$ represents the
192 : * factor needed to scale the vector $\vec x - \vec P$ to intersect with the
193 : * infinite planes that contain the six faces of the cube.
194 : *
195 : * Some $\lambda_c^{\pm k}$ will be negative as they point to the opposite
196 : * planes. We can discard those. Then for the $\lambda_c^{\pm k}$ that are
197 : * positive, we want the one that is smallest in magnitude. This can be written
198 : * as
199 : *
200 : * \begin{equation}\label{eq:lambda_c}
201 : * \lambda_c = \min_{\substack{\lambda_c^{\pm k} > 0}} (\lambda_c^{\pm k})
202 : * \end{equation}
203 : *
204 : * #### Computing sphere lambda
205 : *
206 : * To compute $\lambda_s$ for the sphere case, again we notice that (in the
207 : * image)
208 : *
209 : * \begin{align}
210 : * R_1^2 &= |\vec x_1|^2 \\
211 : * &= |\lambda_s (\vec x - \vec P) + \vec P|^2 \\
212 : * &= \lambda_s^2 |\vec x - \vec P|^2 + 2 \lambda_s (\vec x - \vec P) \cdot
213 : * \vec P + |\vec P|^2
214 : * \end{align}
215 : *
216 : * This is simply a quadratic equation for $\lambda_s$:
217 : *
218 : * \begin{equation}\label{eq:lambda_s}
219 : * 0 = \lambda_s^2 |\vec x - \vec P|^2 + 2 \lambda_s (\vec x - \vec P) \cdot
220 : * \vec P + |\vec P|^2 - R_1^2
221 : * \end{equation}
222 : *
223 : * We explicitly don't write the typical form of the solution of a quadratic
224 : * equation because that form behaviors poorly numerically if the two terms in
225 : * the numerator are close to each other (bad numerical roundoff). Therefore, we
226 : * just resolve to compute the solution numerically, and say that we have
227 : * $\lambda_s$.
228 : *
229 : * We have now solved for $\vec x_1$.
230 : *
231 : * \begin{equation}
232 : * \vec x_1 = \left((1 - s_1) \lambda_c + s_1 \lambda_s \right) (\vec x -
233 : * \vec P) + \vec P
234 : * \end{equation}
235 : *
236 : * given the solutions for $\lambda_c$ in Eqn. $\ref{eq:lambda_c}$ and
237 : * $\lambda_s$ in Eqn. $\ref{eq:lambda_s}$.
238 : *
239 : * ## Original radius divided by mapped radius
240 : *
241 : * Given an already mapped point $\tilde{\vec x}$ and the radial distortion
242 : * $\Sigma(\theta, \phi) = \sum_{\ell,m}\lambda_{\ell m}Y_{\ell m}(\theta,
243 : * \phi)$, we can figure out the ratio of the original radius $r$ to the mapped
244 : * $\tilde{r} = |\tilde{\vec x}|$ by solving
245 : *
246 : * \begin{equation}
247 : * \frac{r}{\tilde{r}} =
248 : * \frac{1}{1-G(r)\Sigma(\theta,\phi)}
249 : * \end{equation}
250 : *
251 : * If the axis is `Axis::Interior`, this formula becomes (after simplification
252 : * and because the inner surface must be spherical)
253 : *
254 : * \begin{equation}
255 : * 0 = \left(\frac{r}{\tilde{r}}\right)^3 - \left(\frac{r}{\tilde{r}}\right)
256 : * \frac{R_0^3}{\tilde{r}^2\Sigma(\theta,\phi)} +
257 : * \frac{R_0^3}{\tilde{r}^2\Sigma(\theta,\phi)}
258 : * \end{equation}
259 : *
260 : * This is a special case of a cubic root, so we use the method outlined in
261 : * \cite NumericalRecipes on pg 228.
262 : *
263 : * For other axes, after plugging in $f(r,\theta,\phi)/r$ and solving, we get
264 : *
265 : * \begin{equation}
266 : * \frac{r}{\tilde{r}} = \frac{1 + \frac{|\vec x_1 - \vec P|\Sigma(\theta,
267 : * \phi)}{\tilde{r}(|\vec x_1 - \vec P| - |\vec x_0 - \vec P|)}}{1 +
268 : * \frac{\Sigma(\theta, \phi)}{|\vec x_1 - \vec P| - |\vec x_0 - \vec P|}}
269 : * \label{eq:r_over_rtil}
270 : * \end{equation}
271 : *
272 : * \note Since $\vec x_0 - \vec P$ and $\vec x_1 - \vec P$ are not functions of
273 : * the radius, we can use $\tilde{\vec x}$ in Eq. $\ref{eq:x_0_vector}$ instead
274 : * of $\vec x$.
275 : *
276 : * \parblock
277 : *
278 : * \note If \p reverse is true, then the value multiplying $\Sigma$ in the
279 : * numerator is now $-|\vec x_0 - \vec P|$ and in the denominator $\Sigma$ picks
280 : * up a minus sign factor.
281 : *
282 : * \endparblock
283 : *
284 : * ## Gradient
285 : *
286 : * The cartesian gradient of the transition function is
287 : *
288 : * \begin{equation}
289 : * \frac{\partial G}{\partial x_i} = \frac{1}{r}\frac{\partial f}{\partial
290 : * \xi^i} - \frac{f \xi_i}{r^3}
291 : * \end{equation}
292 : *
293 : * where
294 : *
295 : * \begin{equation}
296 : * \frac{\partial f}{\partial x_i} = \frac{\frac{\partial
297 : * |\vec x_1 - \vec P|}{\partial x_i} - \frac{x_i - P_i}{r}}{|\vec x_1 - \vec P|
298 : * - |\vec x_0 - \vec P|} - \frac{\left(|\vec x_1 - \vec P| -
299 : * r\right)\left(\frac{\partial |\vec x_1 - \vec P|}{\partial x_i} -
300 : * \frac{\partial |\vec x_0 - \vec P|}{\partial x_i} \right)}{\left(|\vec x_1 -
301 : * \vec P| - |\vec x_0 - \vec P|\right)^2}.
302 : * \end{equation}
303 : *
304 : * \note If \p reverse is true, the gradient picks up an overall factor of -1.0.
305 : *
306 : * If the axis is `Axis::Interior`, then the gradient is
307 : *
308 : * \begin{equation}
309 : * \frac{\partial G}{\partial x_i} = \frac{2(x_i-P_i)}{R_0^3}
310 : * \end{equation}
311 : *
312 : * For the other axes, we need to compute the gradients of $\vec x_0$ and $\vec
313 : * x_1$.
314 : *
315 : * ### Gradient of vector to inner surface
316 : *
317 : * If $\vec P \neq 0$, then we require the inner surface to be a sphere, which
318 : * means the gradient is trivially zero. Therefore the following is only for
319 : * when $\vec P = 0$.
320 : *
321 : * The gradient of $|\vec x_0 - \vec P|$ depends on the result of the $\max$ in
322 : * the denominator of $\ref{eq:x_0_norm}$. To simplify the expression, let $i
323 : * \in \{0,1,2\}$ correspond to the $\max(x_i - P_i)$. Then we can write the
324 : * gradient of as
325 : *
326 : * \begin{align}
327 : * \frac{\partial |\vec x_0 - \vec P|}{\partial x_i} &= -\text{sgn}(x_i -
328 : * P_i)\frac{R(1-s)\left(r^2 - (x_i - P_i)^2\right)}
329 : * {r (x_i - P_i)^2 \sqrt{3}} \\
330 : * \frac{\partial |\vec x_0 - \vec P|}{\partial x_{i+1}} &=
331 : * \frac{R(1-s)(x_{i+1} - P_{i+1})}{r |x_i - P_i| \sqrt{3}} \\
332 : * \frac{\partial |\vec x_0 - \vec P|}{\partial x_{i+2}} &=
333 : * \frac{R(1-s)(x_{i+2} - P_{i+2})}{r |x_i - P_i| \sqrt{3}}
334 : * \end{align}
335 : *
336 : * where $i+1$ and $i+2$ are understood to be $\mod 3$. Also since we don't
337 : * allow points at the center of the inner surfaces, we can be assured that for
338 : * whichever $i$ is max, that $x_i - P_i$ won't be zero.
339 : *
340 : * ### Gradient of vector to outer surface
341 : *
342 : * This gradient is much more complicated because of how we have to define
343 : * $\lambda$. To start off with, we need to compute
344 : *
345 : * \begin{align}\label{eq:x_1_grad}
346 : * \frac{\partial}{\partial x_i}|\vec x_1 - \vec P| &=
347 : * \frac{\partial}{\partial x_i} (\lambda |\vec x - \vec P|) \\
348 : * &= \frac{\lambda (x_i - P_i)}{|\vec x - \vec P|} + \frac{\partial
349 : * \lambda}{\partial x_i}|\vec x - \vec P|
350 : * \end{align}
351 : *
352 : * Since we know $\lambda$ and $\vec x$, that just leaves
353 : *
354 : * \begin{equation}
355 : * \frac{\partial\lambda}{\partial x_i} = \frac{\partial\lambda_c}{\partial
356 : * x_i}(1-s_1) + s_1\frac{\partial\lambda_s}{\partial x_i}
357 : * \end{equation}
358 : *
359 : * #### Gradient of cube lambda
360 : *
361 : * Taking the gradient of Eqn. $\ref{eq:lambda_c_pm}$ we notice that $L$, $\vec
362 : * P$, and $\hat x_k$ are constant, so we only have to worry about $\vec x$
363 : *
364 : * \begin{equation}\label{eq:tmp_lambda_c_grad}
365 : * \frac{\partial}{\partial x_i}\lambda_c^{\pm k} = -\frac{(\pm L -
366 : * \vec P \cdot \hat x_k)(\partial/\partial x_i(\vec x - \vec P) \cdot \hat
367 : * x_k)}{((\vec x - \vec P) \cdot \hat x_k)^2}
368 : * \end{equation}
369 : *
370 : * The gradient in the numerator is simple to take so we end up with
371 : *
372 : * \begin{equation}
373 : * \frac{\partial}{\partial x_i}\lambda_c^{\pm k} = -\frac{(\pm L -
374 : * \vec P \cdot \hat x_k)}{( x_k - P_k)^2} \delta_{ik}
375 : * \end{equation}
376 : *
377 : * #### Gradient of sphere lambda
378 : *
379 : * Taking the gradient of Eqn. $\ref{eq:lambda_s}$ (and using the fact that
380 : * $\vec P$ and $R$ are constant), we get
381 : *
382 : * \begin{equation}
383 : * 0 = \lambda_s \frac{\partial \lambda_s}{\partial x_i}|\vec x - \vec P|^2
384 : * + \lambda_s^2 (\vec x - \vec P) \cdot \frac{\partial}{\partial x_i}(\vec x -
385 : * \vec P) + \frac{\partial\lambda_s}{\partial x_i}(\vec x - \vec P) \cdot \vec
386 : * P + \lambda_s \vec P \cdot \frac{\partial}{\partial x_i}(\vec x - \vec P)
387 : * \end{equation}
388 : *
389 : * Taking the simple gradients and moving things around, we get
390 : *
391 : * \begin{equation}
392 : * \frac{\partial\lambda_s}{\partial x_i} = - \frac{\lambda_s^2 (x_i - P_i)
393 : * + \lambda_s P_i}{\lambda_s |\vec x - \vec P|^2 + (\vec x - \vec P) \cdot \vec
394 : * P}
395 : * \end{equation}
396 : *
397 : * where $\lambda_s$ can be computed from Eqn. $\ref{eq:lambda_s}$.
398 : *
399 : * So now we can put it all together and we get
400 : *
401 : * \begin{equation}
402 : * \frac{\partial\lambda}{\partial x_i} = -\frac{(\pm L - \vec P \cdot
403 : * \hat x_k)}{( x_k - P_k)^2} \delta_{ik}(1 - s_1)
404 : * + s_1 \frac{\lambda_s^2 (x_i - P_i) + \lambda_s P_i}{\lambda_s |\vec x -
405 : * \vec P|^2 + (\vec x - \vec P) \cdot \vec P} \end{equation}
406 : *
407 : * which can then be substituted into Eqn. $\ref{eq:x_1_grad}$ to get the
408 : * gradient of $|\vec x_1 - \vec P|$.
409 : */
410 1 : class Wedge final : public ShapeMapTransitionFunction {
411 0 : struct Surface {
412 0 : std::array<double, 3> center{};
413 0 : double radius{};
414 0 : double sphericity{};
415 0 : double half_cube_length{};
416 :
417 0 : Surface() = default;
418 0 : Surface(const std::array<double, 3>& center_in, double radius_in,
419 : double sphericity_in);
420 :
421 0 : void pup(PUP::er& p);
422 : };
423 :
424 : public:
425 : /*!
426 : * \brief Class to represent the direction of the wedge relative to the outer
427 : * center.
428 : *
429 : * \details \p Interior means within inner surface.
430 : */
431 0 : enum class Axis : int {
432 : Interior = 4,
433 : PlusZ = 3,
434 : MinusZ = -3,
435 : PlusY = 2,
436 : MinusY = -2,
437 : PlusX = 1,
438 : MinusX = -1,
439 : None = 0,
440 : };
441 :
442 0 : friend std::ostream& operator<<(std::ostream& os, Axis axis);
443 :
444 : private:
445 0 : static size_t axis_index(Axis axis);
446 0 : static double axis_sgn(Axis axis);
447 :
448 : public:
449 0 : explicit Wedge() = default;
450 :
451 : /*!
452 : * \brief Construct a Wedge transition for a wedge block in a given direction.
453 : *
454 : * \details Many concentric wedges can be part of the same falloff from 1 at
455 : * the inner boundary of the innermost wedge, to 0 at the outer boundary of
456 : * the outermost wedge.
457 : *
458 : * \note If \p inner_center and \p outer_center are different, then
459 : * \p inner_sphericity must be 1.0. if the \p axis is `Axis::Interior`, then
460 : * \p inner_sphericity must be 1.0 and \p reverse must be false.
461 : *
462 : * \param inner_center Center of the inner surface
463 : * \param inner_radius Inner radius of innermost wedge
464 : * \param inner_sphericity Sphericity of innermost surface of innermost wedge
465 : * \param outer_center Center of the outer surface
466 : * \param outer_radius Outermost radius of outermost wedge
467 : * \param outer_sphericity Sphericity of outermost surface of outermost wedge
468 : * \param axis The direction that this wedge is in.
469 : * \param reverse If true, the function $f$ will be 0 at the inner
470 : * boundary and 1 at the outer boundary (useful for deforming star surfaces).
471 : * Otherwise, it will be 1 at the inner boundary and 0 at the outer boundary
472 : * (useful for deforming black hole excision surfaces).
473 : */
474 1 : Wedge(const std::array<double, 3>& inner_center, double inner_radius,
475 : double inner_sphericity, const std::array<double, 3>& outer_center,
476 : double outer_radius, double outer_sphericity, Axis axis,
477 : bool reverse = false);
478 :
479 1 : double operator()(
480 : const std::array<double, 3>& source_coords,
481 : const std::optional<size_t>& one_over_radius_power) const override;
482 1 : DataVector operator()(
483 : const std::array<DataVector, 3>& source_coords,
484 : const std::optional<size_t>& one_over_radius_power) const override;
485 :
486 1 : std::optional<double> original_radius_over_radius(
487 : const std::array<double, 3>& target_coords,
488 : double radial_distortion) const override;
489 :
490 1 : std::array<double, 3> gradient(
491 : const std::array<double, 3>& source_coords) const override;
492 1 : std::array<DataVector, 3> gradient(
493 : const std::array<DataVector, 3>& source_coords) const override;
494 :
495 0 : WRAPPED_PUPable_decl_template(Wedge);
496 0 : explicit Wedge(CkMigrateMessage* msg);
497 0 : void pup(PUP::er& p) override;
498 :
499 0 : std::unique_ptr<ShapeMapTransitionFunction> get_clone() const override {
500 : return std::make_unique<Wedge>(*this);
501 : }
502 :
503 0 : bool operator==(const ShapeMapTransitionFunction& other) const override;
504 0 : bool operator!=(const ShapeMapTransitionFunction& other) const override;
505 :
506 : private:
507 : template <typename T>
508 0 : T call_impl(const std::array<T, 3>& source_coords,
509 : const std::optional<size_t>& one_over_radius_power) const;
510 :
511 : template <typename T>
512 0 : std::array<T, 3> gradient_impl(const std::array<T, 3>& source_coords) const;
513 :
514 : // This is x_0 - P in the docs
515 : template <typename T>
516 0 : std::array<T, 3> compute_inner_surface_vector(
517 : const std::array<T, 3>& centered_coords,
518 : const T& centered_coords_magnitude,
519 : const std::optional<Axis>& potential_axis) const;
520 :
521 : // This is x_1 - P in the docs
522 : template <typename T>
523 0 : std::array<T, 3> compute_outer_surface_vector(
524 : const std::array<T, 3>& centered_coords, const T& lambda) const;
525 :
526 : template <typename T>
527 0 : T lambda_cube(const std::array<T, 3>& centered_coords,
528 : const std::optional<Axis>& potential_axis) const;
529 :
530 : template <typename T>
531 0 : T lambda_sphere(const std::array<T, 3>& centered_coords,
532 : const T& centered_coords_magnitude) const;
533 :
534 : template <typename T>
535 0 : T compute_lambda(const std::array<T, 3>& centered_coords,
536 : const T& centered_coords_magnitude,
537 : const std::optional<Axis>& potential_axis) const;
538 :
539 : template <typename T>
540 0 : std::array<T, 3> lambda_cube_gradient(
541 : const T& lambda_cube, const std::array<T, 3>& centered_coords) const;
542 :
543 : template <typename T>
544 0 : std::array<T, 3> lambda_sphere_gradient(
545 : const T& lambda_sphere, const std::array<T, 3>& centered_coords,
546 : const T& centered_coords_magnitude) const;
547 :
548 : template <typename T>
549 0 : std::array<T, 3> compute_lambda_gradient(
550 : const std::array<T, 3>& centered_coords,
551 : const T& centered_coords_magnitude) const;
552 :
553 : template <typename T>
554 0 : void check_distances(const T& inner_distance, const T& outer_distance,
555 : const T& centered_coords_magnitude,
556 : const std::array<T, 3>& source_coords,
557 : bool check_bounds) const;
558 :
559 0 : Surface inner_surface_;
560 0 : Surface outer_surface_;
561 0 : std::array<double, 3> projection_center_{};
562 0 : Axis axis_{};
563 0 : bool reverse_{false};
564 :
565 0 : static constexpr double eps_ = std::numeric_limits<double>::epsilon() * 100;
566 : };
567 : } // namespace domain::CoordinateMaps::ShapeMapTransitionFunctions
|
