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