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 <optional>
10 :
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Utilities/Gsl.hpp"
13 : #include "Utilities/MakeArray.hpp"
14 :
15 : /// \cond
16 : namespace PUP {
17 : class er;
18 : } // namespace PUP
19 : /// \endcond
20 :
21 : namespace domain::CoordinateMaps::FocallyLiftedInnerMaps {
22 : /*!
23 : * \brief A FocallyLiftedInnerMap that maps a 3D unit right cylindrical shell
24 : * to a volume that connects portions of two spherical surfaces.
25 : *
26 : * \details The domain of the map is a 3D unit right cylinder with
27 : * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
28 : * \f$-1\leq\bar{z}\leq 1\f$ and \f$1\leq \bar{x}^2+\bar{y}^2 \leq
29 : * 4\f$. The range of the map has coordinates \f$(x,y,z)\f$.
30 : *
31 : * Consider a sphere with center \f$C^i\f$ and radius \f$R\f$ that is
32 : * intersected by two planes normal to the \f$z\f$ axis located at
33 : * \f$z = z_\mathrm{L}\f$ and \f$z = z_\mathrm{U}\f$, with
34 : * \f$z_\mathrm{L} < z_\mathrm{U}\f$.
35 : * `Side` provides the following functions:
36 : *
37 : * ### forward_map()
38 : * `forward_map()` maps \f$(\bar{x},\bar{y},\bar{z})\f$ to a point on the inner
39 : * surface
40 : * \f$\bar{x}^2+\bar{y}^2=1\f$ by dividing \f$\bar{x}\f$ and \f$\bar{y}\f$
41 : * by \f$(1+\sigma)\f$, where \f$\sigma\f$ is the function given by Eq. (7)
42 : * below.
43 : * Then it maps that point to a point on the portion of the sphere with
44 : * \f$z_\mathrm{L} \leq z \leq z_\mathrm{U}\f$.
45 : * `forward_map()` returns
46 : * \f$x_0^i\f$, the 3D coordinates on that sphere, which are given by
47 : *
48 : * \f{align}
49 : * x_0^0 &= R \sin\theta \frac{\bar{x}}{1+\sigma} + C^0,\\
50 : * x_0^1 &= R \sin\theta \frac{\bar{y}}{1+\sigma} + C^1,\\
51 : * x_0^2 &= R \cos\theta + C^2.\\
52 : * \f}
53 : *
54 : * Here
55 : * \f{align}
56 : * \theta = \theta_\mathrm{max} +
57 : * (\theta_\mathrm{min}-\theta_\mathrm{max}) \frac{\bar{z}+1}{2},
58 : * \f}
59 : *
60 : * where
61 : * \f{align}
62 : * \cos(\theta_\mathrm{max}) &= (z_\mathrm{L}-C^2)/R,\\
63 : * \cos(\theta_\mathrm{min}) &= (z_\mathrm{U}-C^2)/R.
64 : * \f}
65 : *
66 : * Note that \f$\theta\f$ decreases with increasing \f$\bar{z}\f$,
67 : * which is the usual convention for a polar angle but might otherwise
68 : * cause confusion.
69 : *
70 : * ### sigma
71 : *
72 : * \f$\sigma\f$ is a function that is zero on the sphere
73 : * \f$x^i=x_0^i\f$ and unity at \f$\bar{x}^2+\bar{y}^2=4\f$
74 : * (corresponding to the upper surface of the FocallyLiftedMap). We define
75 : *
76 : * \f{align}
77 : * \sigma &= \sqrt{\bar{x}^2+\bar{y}^2}-1.
78 : * \f}
79 : *
80 : * ### deriv_sigma
81 : *
82 : * `deriv_sigma` returns
83 : *
84 : * \f{align}
85 : * \frac{\partial \sigma}{\partial \bar{x}^j} &=
86 : * \left(\frac{\bar{x}}{1+\sigma},
87 : * \frac{\bar{y}}{1+\sigma},0\right).
88 : * \f}
89 : *
90 : * ### jacobian
91 : *
92 : * `jacobian` returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
93 : * The arguments to `jacobian` are \f$(\bar{x},\bar{y},\bar{z})\f$.
94 : * Differentiating Eqs.(1--4) above yields
95 : *
96 : * \f{align*}
97 : * \frac{\partial x_0^0}{\partial \bar{x}} &= R \sin\theta
98 : * \frac{\bar{y}^2}{(1+\sigma)^3}, \\
99 : * \frac{\partial x_0^0}{\partial \bar{y}} &= -R \sin\theta
100 : * \frac{\bar{x}\bar{y}}{(1+\sigma)^3}, \\
101 : * \frac{\partial x_0^0}{\partial \bar{z}} &=
102 : * R \cos\theta \frac{\theta_\mathrm{min}-\theta_\mathrm{max}}{2(1+\sigma)}
103 : * \bar{x},\\
104 : * \frac{\partial x_0^1}{\partial \bar{x}} &= -R \sin\theta
105 : * \frac{\bar{x}\bar{y}}{(1+\sigma)^3}, \\
106 : * \frac{\partial x_0^1}{\partial \bar{y}} &= R \sin\theta
107 : * \frac{\bar{x}^2}{(1+\sigma)^3}, \\
108 : * \frac{\partial x_0^1}{\partial \bar{z}} &=
109 : * R \cos\theta \frac{\theta_\mathrm{min}-\theta_\mathrm{max}}{2(1+\sigma)}
110 : * \bar{y},\\
111 : * \frac{\partial x_0^2}{\partial \bar{x}} &= 0,\\
112 : * \frac{\partial x_0^2}{\partial \bar{y}} &= 0,\\
113 : * \frac{\partial x_0^2}{\partial \bar{z}} &=
114 : * - R \sin\theta \frac{\theta_\mathrm{min}-\theta_\mathrm{max}}{2}.
115 : * \f}
116 : *
117 : * ### inverse
118 : *
119 : * `inverse` takes \f$x_0^i\f$ and \f$\sigma\f$ as arguments, and
120 : * returns \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed
121 : * `std::optional<std::array<double, 3>>` if \f$x_0^i\f$ or \f$\sigma\f$ are
122 : * outside the range of the map.
123 : *
124 : * If \f$\sigma\f$ is outside the range \f$[0,1]\f$ then we return
125 : * a default-constructed `std::optional<std::array<double, 3>>`.
126 : *
127 : * To get \f$\bar{z}\f$ we invert Eq. (4):
128 : * \f{align}
129 : * \bar{z} &= 2\frac{\acos\left((x_0^2-C^2)/R\right)-\theta_\mathrm{max}}
130 : * {\theta_\mathrm{min}-\theta_\mathrm{max}} - 1.
131 : * \f}
132 : *
133 : * If \f$\bar{z}\f$ is outside the range \f$[-1,1]\f$ then we return
134 : * a default-constructed `std::optional<std::array<double, 3>>`.
135 : *
136 : * To compute \f$\bar{x}\f$ and \f$\bar{y}\f$, we invert Eqs. (1--3) and
137 : * use \f$\sigma\f$:
138 : *
139 : * \f{align}
140 : * \bar{x} &= \frac{(x_0^0-C^0) (1+\sigma)}{\rho},\\
141 : * \bar{y} &= \frac{(x_0^1-C^1) (1+\sigma)}{\rho},
142 : * \f}
143 : *
144 : * where
145 : *
146 : * \f{align}
147 : * \rho = \sqrt{(x_0^0-C^0)^2+(x_0^1-C^1)^2}.
148 : * \f}
149 : *
150 : * ### lambda_tilde
151 : *
152 : * `lambda_tilde` takes as arguments a point \f$x^i\f$ and a projection point
153 : * \f$P^i\f$, and computes \f$\tilde{\lambda}\f$, the solution to
154 : *
155 : * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
156 : *
157 : * Since \f$x_0^i\f$ must lie on the sphere, \f$\tilde{\lambda}\f$ is the
158 : * solution of the quadratic equation
159 : *
160 : * \f{align}
161 : * |P^i + (x^i - P^i) \tilde{\lambda} - C^i |^2 - R^2 = 0.
162 : * \f}
163 : *
164 : * In solving the quadratic, we choose the larger root if
165 : * \f$x^2>z_\mathrm{P}\f$ and the smaller root otherwise. We demand
166 : * that the root is greater than unity. If there is no such root,
167 : * this means that the point \f$x^i\f$ is not in the range of the map
168 : * so we return a default-constructed `std::optional<double>`.
169 : *
170 : * ### deriv_lambda_tilde
171 : *
172 : * `deriv_lambda_tilde` takes as arguments \f$x_0^i\f$, a projection point
173 : * \f$P^i\f$, and \f$\tilde{\lambda}\f$, and
174 : * returns \f$\partial \tilde{\lambda}/\partial x^i\f$.
175 : * By differentiating Eq. (14), we find
176 : *
177 : * \f{align}
178 : * \frac{\partial\tilde{\lambda}}{\partial x^j} &=
179 : * \tilde{\lambda}^2 \frac{C^j - x_0^j}{
180 : * (x_0^i - P^i)(x_{0i} - C_{i})} \nonumber \\
181 : * &= \tilde{\lambda}^2 \frac{C^j - x_0^j}{|x_0^i - P^i|^2
182 : * + (x_0^i - P^i)(P_i - C_{i})}.
183 : * \f}
184 : *
185 : * ### inv_jacobian
186 : *
187 : * `inv_jacobian` returns \f$\partial \bar{x}^i/\partial x_0^k\f$,
188 : * where \f$\sigma\f$ is held fixed.
189 : * The arguments to `inv_jacobian` are \f$(\bar{x},\bar{y},\bar{z})\f$.
190 : *
191 : * Note from Eqs. (9--12) that \f$\bar{x}\f$ and \f$\bar{y}\f$
192 : * depend only on \f$x_0^0\f$ and \f$x_0^1\f$ but not on \f$x_0^2\f$.
193 : *
194 : * By differentiating Eqs. (9--12), we find
195 : *
196 : * \f{align*}
197 : * \frac{\partial \bar{x}}{\partial x_0^0} &=
198 : * \frac{\bar{y}^2}{(1+\sigma)\rho},\\
199 : * \frac{\partial \bar{x}}{\partial x_0^1} &=
200 : * - \frac{\bar{x}\bar{y}}{(1+\sigma)\rho},\\
201 : * \frac{\partial \bar{x}}{\partial x_0^2} &= 0,\\
202 : * \frac{\partial \bar{y}}{\partial x_0^0} &=
203 : * - \frac{\bar{x}\bar{y}}{(1+\sigma)\rho},\\
204 : * \frac{\partial \bar{y}}{\partial x_0^1} &=
205 : * \frac{\bar{x}^2}{(1+\sigma)\rho},\\
206 : * \frac{\partial \bar{y}}{\partial x_0^2} &= 0,\\
207 : * \frac{\partial \bar{z}}{\partial x_0^0} &= 0,\\
208 : * \frac{\partial \bar{z}}{\partial x_0^1} &= 0,\\
209 : * \frac{\partial \bar{z}}{\partial x_0^2} &=
210 : * -\frac{2}{\rho(\theta_\mathrm{min}-\theta_\mathrm{max})},
211 : * \f}
212 : *
213 : * where
214 : *
215 : * \f[
216 : * \rho = R \sin\theta = R\sin\left(\theta_\mathrm{max} +
217 : * (\theta_\mathrm{min}-\theta_\mathrm{max}) \frac{\bar{z}+1}{2}\right),
218 : * \f]
219 : *
220 : * which is also equal to the quantity in Eq. (12).
221 : *
222 : * ### dxbar_dsigma
223 : *
224 : * `dxbar_dsigma` returns \f$\partial \bar{x}^i/\partial \sigma\f$,
225 : * where \f$x_0^i\f$ is held fixed.
226 : *
227 : * From Eqs. (10) and (11) we have
228 : *
229 : * \f{align}
230 : * \frac{\partial \bar{x}^i}{\partial \sigma} &=
231 : * \left(\frac{\bar{x}}{\sqrt{\bar{x}^2+\bar{y}^2}},
232 : * \frac{\bar{y}}{\sqrt{\bar{x}^2+\bar{y}^2}},0\right).
233 : * \f}
234 : *
235 : */
236 1 : class Side {
237 : public:
238 0 : static constexpr size_t dim = 3;
239 0 : Side(const std::array<double, 3>& center, const double radius,
240 : const double z_lower, const double z_upper);
241 :
242 0 : Side() = default;
243 0 : ~Side() = default;
244 0 : Side(Side&&) = default;
245 0 : Side(const Side&) = default;
246 0 : Side& operator=(const Side&) = default;
247 0 : Side& operator=(Side&&) = default;
248 :
249 : template <typename T>
250 0 : void forward_map(gsl::not_null<std::array<T, 3>*> target_coords,
251 : const std::array<T, 3>& source_coords) const;
252 :
253 0 : std::optional<std::array<double, 3>> inverse(
254 : const std::array<double, 3>& target_coords, double sigma_in) const;
255 :
256 : template <typename T>
257 0 : void jacobian(gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> jacobian_out,
258 : const std::array<T, 3>& source_coords) const;
259 :
260 : template <typename T>
261 0 : void inv_jacobian(
262 : gsl::not_null<tnsr::Ij<T, 3, Frame::NoFrame>*> inv_jacobian_out,
263 : const std::array<T, 3>& source_coords) const;
264 :
265 : template <typename T>
266 0 : void sigma(gsl::not_null<T*> sigma_out,
267 : const std::array<T, 3>& source_coords) const;
268 :
269 : template <typename T>
270 0 : void deriv_sigma(gsl::not_null<std::array<T, 3>*> deriv_sigma_out,
271 : const std::array<T, 3>& source_coords) const;
272 :
273 : template <typename T>
274 0 : void dxbar_dsigma(gsl::not_null<std::array<T, 3>*> dxbar_dsigma_out,
275 : const std::array<T, 3>& source_coords) const;
276 :
277 0 : std::optional<double> lambda_tilde(
278 : const std::array<double, 3>& parent_mapped_target_coords,
279 : const std::array<double, 3>& projection_point,
280 : bool source_is_between_focus_and_target) const;
281 :
282 : template <typename T>
283 0 : void deriv_lambda_tilde(
284 : gsl::not_null<std::array<T, 3>*> deriv_lambda_tilde_out,
285 : const std::array<T, 3>& target_coords, const T& lambda_tilde,
286 : const std::array<double, 3>& projection_point) const;
287 :
288 : // NOLINTNEXTLINE(google-runtime-references)
289 0 : void pup(PUP::er& p);
290 :
291 0 : static bool is_identity() { return false; }
292 :
293 0 : static constexpr bool supports_hessian{false};
294 :
295 : private:
296 0 : friend bool operator==(const Side& lhs, const Side& rhs);
297 0 : std::array<double, 3> center_{
298 : make_array<3>(std::numeric_limits<double>::signaling_NaN())};
299 0 : double radius_{std::numeric_limits<double>::signaling_NaN()};
300 0 : double theta_min_{std::numeric_limits<double>::signaling_NaN()};
301 0 : double theta_max_{std::numeric_limits<double>::signaling_NaN()};
302 : };
303 0 : bool operator!=(const Side& lhs, const Side& rhs);
304 : } // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps
|
