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/TypeTraits/RemoveReferenceWrapper.hpp"
13 :
14 : /// \cond
15 : namespace PUP {
16 : class er;
17 : } // namespace PUP
18 : /// \endcond
19 :
20 : namespace domain::CoordinateMaps {
21 :
22 : template <typename InnerMap>
23 : class FocallyLiftedMap;
24 :
25 : template <typename InnerMap>
26 0 : bool operator==(const FocallyLiftedMap<InnerMap>& lhs,
27 : const FocallyLiftedMap<InnerMap>& rhs);
28 :
29 : /*!
30 : * \ingroup CoordinateMapsGroup
31 : *
32 : * \brief Map from \f$(\bar{x},\bar{y},\bar{z})\f$ to the volume
33 : * contained between a 2D surface and the surface of a 2-sphere.
34 : *
35 : *
36 : * \image html FocallyLifted.svg "2D representation of focally lifted map."
37 : *
38 : * \details We are given the radius \f$R\f$ and
39 : * center \f$C^i\f$ of a sphere, and a projection point \f$P^i\f$. Also,
40 : * through the class defined by the `InnerMap` template parameter,
41 : * we are given the functions
42 : * \f$f^i(\bar{x},\bar{y},\bar{z})\f$ and
43 : * \f$\sigma(\bar{x},\bar{y},\bar{z})\f$ defined below; these
44 : * functions define the mapping from \f$(\bar{x},\bar{y},\bar{z})\f$
45 : * to the 2D surface.
46 : *
47 : * The above figure is a 2D representation of the focally lifted map,
48 : * where the shaded region in \f$\bar{x}^i\f$ coordinates on the left
49 : * is mapped to the shaded region in the \f$x^i\f$ coordinates on the
50 : * right. Shown is an arbitrary point \f$\bar{x}^i\f$ that is mapped
51 : * to \f$x^i\f$; for that point, the corresponding \f$x_0^i\f$ and
52 : * \f$x_1^i\f$ (defined below) are shown on the right, as is the
53 : * projection point \f$P^i\f$. Also shown are the level surfaces
54 : * \f$\sigma=0\f$ and \f$\sigma=1\f$.
55 : *
56 : * The input coordinates are labeled \f$(\bar{x},\bar{y},\bar{z})\f$.
57 : * Let \f$x_0^i = f^i(\bar{x},\bar{y},\bar{z})\f$ be the three coordinates
58 : * of the 2D surface in 3D space.
59 : *
60 : * Now let \f$x_1^i\f$ be a point on the surface of the sphere,
61 : * constructed so that \f$P^i\f$, \f$x_0^i\f$, and \f$x_1^i\f$ are
62 : * co-linear. In particular, \f$x_1^i\f$ is determined by the
63 : * equation
64 : *
65 : * \f{align} x_1^i = P^i + (x_0^i - P^i) \lambda,\f}
66 : *
67 : * where \f$\lambda\f$ is a scalar factor that depends on \f$x_0^i\f$ and
68 : * that can be computed by solving a quadratic equation. This quadratic
69 : * equation is derived by demanding that \f$x_1^i\f$ lies on the sphere:
70 : *
71 : * \f{align}
72 : * |P^i + (x_0^i - P^i) \lambda - C^i |^2 - R^2 = 0,
73 : * \f}
74 : *
75 : * where \f$|A^i|^2\f$ means \f$\delta_{ij} A^i A^j\f$.
76 : *
77 : * The quadratic equation, Eq. (2),
78 : * takes the usual form \f$a\lambda^2+b\lambda+c=0\f$,
79 : * with
80 : *
81 : * \f{align*}
82 : * a &= |x_0^i-P^i|^2,\\
83 : * b &= 2(x_0^i-P^i)(P^j-C^j)\delta_{ij},\\
84 : * c &= |P^i-C^i|^2 - R^2.
85 : * \f}
86 : *
87 : * Now assume that \f$x_0^i\f$ lies on a level surface defined
88 : * by some scalar function \f$\sigma(\bar{x}^i)=0\f$,
89 : * where \f$\sigma\f$ is normalized so that \f$\sigma=1\f$ on the sphere.
90 : * Then, once \f$\lambda\f$ has been computed and \f$x_1^i\f$ has
91 : * been determined, the map is given by
92 : *
93 : * \f{align}x^i = x_0^i + (x_1^i - x_0^i) \sigma(\bar{x}^i).\f}
94 : *
95 : * Note that classes using FocallyLiftedMap will place restrictions on
96 : * \f$P^i\f$, \f$C^i\f$, \f$x_0^i\f$, and \f$R\f$. For example, we
97 : * demand that \f$P^i\f$ does not lie on either the \f$\sigma=0\f$ or
98 : * \f$\sigma=1\f$ surfaces depicted in the right panel of the figure,
99 : * and we demand that the \f$\sigma=0\f$ and \f$\sigma=1\f$ surfaces do
100 : * not intersect; otherwise the map is singular.
101 : *
102 : * Also note that the quadratic Eq. (2) typically has more than one
103 : * root, corresponding to two intersections of the sphere. The
104 : * boolean parameter `source_is_between_focus_and_target` that is
105 : * passed into the constructor of `FocallyLiftedMap` is used to choose the
106 : * appropriate root, or to error if a suitable
107 : * root is not found. `source_is_between_focus_and_target` should be
108 : * true if the source point lies between the projection point
109 : * \f$P^i\f$ and the sphere. `source_is_between_focus_and_target` is
110 : * known only by each particular `CoordinateMap` that uses
111 : * `FocallyLiftedMap`.
112 : *
113 : * ### Jacobian
114 : *
115 : * Differentiating Eq. (1) above yields
116 : *
117 : * \f{align}
118 : * \frac{\partial x_1^i}{\partial x_0^j} = \lambda \delta^i_j +
119 : * (x_0^i - P^i) \frac{\partial \lambda}{\partial x_0^j}.
120 : * \f}
121 : *
122 : * and differentiating Eq. (2) and then inserting Eq. (1) yields
123 : *
124 : * \f{align}
125 : * \frac{\partial\lambda}{\partial x_0^j} &=
126 : * \lambda^2 \frac{C_j - x_1^j}{|x_1^i - P^i|^2
127 : * + (x_1^i - P^i)(P_i - C_{i})}.
128 : * \f}
129 : *
130 : * The Jacobian can be found by differentiating Eq. (3) above and using the
131 : * chain rule, recognizing that \f$x_1^i\f$ depends on \f$\bar{x}^i\f$ only
132 : * through its dependence on \f$x_0^i\f$; this is because there is no explicit
133 : * dependence on \f$\bar{x}^i\f$ in Eq. (1) (which determines \f$x_1^i\f$)
134 : * or Eq. (2) (which determines \f$\lambda\f$).
135 : *
136 : * \f{align}
137 : * \frac{\partial x^i}{\partial \bar{x}^j} &=
138 : * \sigma \frac{\partial x_1^i}{\partial x_0^k}
139 : * \frac{\partial x_0^k}{\partial \bar{x}^j} +
140 : * (1-\sigma)\frac{\partial x_0^i}{\partial \bar{x}^j}
141 : * + \frac{\partial \sigma}{\partial \bar{x}^j} (x_1^i - x_0^i),
142 : * \nonumber \\
143 : * &= (1-\sigma+\lambda\sigma) \frac{\partial x_0^i}{\partial \bar{x}^j} +
144 : * \sigma (x_0^i - P^i) \frac{\partial \lambda}{\partial x_0^k}
145 : * \frac{\partial x_0^k}{\partial \bar{x}^j}
146 : * + \frac{\partial \sigma}{\partial \bar{x}^j} (x_1^i - x_0^i),
147 : * \f}
148 : * where in the last line we have substituted Eq. (4).
149 : *
150 : * The class defined by the `InnerMap` template parameter
151 : * provides the function `deriv_sigma`, which returns
152 : * \f$\partial \sigma/\partial \bar{x}^j\f$, and the function `jacobian`,
153 : * which returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
154 : *
155 : * ### Inverse map.
156 : *
157 : * Given \f$x^i\f$, we wish to compute \f$\bar{x}^i\f$.
158 : *
159 : * We first find the coordinates \f$x_0^i\f$ that lie on the 2-surface
160 : * and are defined such that \f$P^i\f$, \f$x_0^i\f$,
161 : * and \f$x^i\f$ are co-linear.
162 : * See the right panel of the above figure.
163 : * \f$x_0^i\f$ is determined by the equation
164 : *
165 : * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda},\f}
166 : *
167 : * where \f$\tilde{\lambda}\f$ is a scalar factor that depends on
168 : * \f$x^i\f$ and is determined by the class defined by the `InnerMap` template
169 : * parameter. `InnerMap`
170 : * provides a function `lambda_tilde` that takes \f$x^i\f$ and
171 : * \f$P^i\f$ as arguments and returns \f$\tilde{\lambda}\f$ (or
172 : * a default-constructed `std::optional` if the appropriate
173 : * \f$\tilde{\lambda}\f$ cannot be
174 : * found; this default-constructed value indicates that the point \f$x^i\f$ is
175 : * outside the range of the map).
176 : *
177 : * Now consider the coordinates \f$x_1^i\f$ that lie on the sphere and
178 : * are defined such that \f$P^i\f$, \f$x_1^i\f$, and \f$x^i\f$ are
179 : * co-linear. See the right panel of the figure.
180 : * \f$x_1^i\f$ is determined by the equation
181 : *
182 : * \f{align} x_1^i = P^i + (x^i - P^i) \bar{\lambda},\f}
183 : *
184 : * where \f$\bar{\lambda}\f$ is a scalar factor that depends on \f$x^i\f$ and
185 : * is the solution of a quadratic
186 : * that is derived by demanding that \f$x_1^i\f$ lies on the sphere:
187 : *
188 : * \f{align}
189 : * |P^i + (x^i - P^i) \bar{\lambda} - C^i |^2 - R^2 = 0.
190 : * \f}
191 : *
192 : * Eq. (9) is a quadratic equation that
193 : * takes the usual form \f$a\bar{\lambda}^2+b\bar{\lambda}+c=0\f$,
194 : * with
195 : *
196 : * \f{align*}
197 : * a &= |x^i-P^i|^2,\\
198 : * b &= 2(x^i-P^i)(P^j-C^j)\delta_{ij},\\
199 : * c &= |P^i-C^i|^2 - R^2.
200 : * \f}
201 : *
202 : * Note that we don't actually need to compute \f$x_1^i\f$. Instead, we
203 : * can determine \f$\sigma\f$ by the relation (obtained by solving Eq. (3)
204 : * for \f$\sigma\f$ and then inserting Eqs. (7) and (8) to eliminate
205 : * \f$x_1^i\f$ and \f$x_0^i\f$)
206 : *
207 : * \f{align}
208 : * \sigma = \frac{\tilde{\lambda}-1}{\tilde{\lambda}-\bar{\lambda}}.
209 : * \f}
210 : *
211 : * The denominator of Eq. (10) is nonzero for nonsingular maps:
212 : * From Eqs. (7) and (8), \f$\bar{\lambda}=\tilde{\lambda}\f$ means
213 : * that \f$x_1^i=x_0^i\f$, which means that
214 : * the \f$\sigma=0\f$ and \f$\sigma=1\f$ surfaces intersect, i.e.
215 : * the map is singular.
216 : *
217 : * Once we have \f$x_0^i\f$ and \f$\sigma\f$, the point
218 : * \f$(\bar{x},\bar{y},\bar{z})\f$ is uniquely determined by `InnerMap`.
219 : * The `inverse` function of `InnerMap` takes \f$x_0^i\f$ and \f$\sigma\f$
220 : * as arguments, and returns
221 : * \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed `std::optional`
222 : * if \f$x_0^i\f$ or
223 : * \f$\sigma\f$ are outside the range of the map.
224 : *
225 : * #### Root polishing
226 : *
227 : * The inverse function described above will sometimes have errors that
228 : * are noticeably larger than roundoff. Therefore we apply a single
229 : * Newton-Raphson iteration to refine the result of the inverse map:
230 : * Suppose we are given \f$x^i\f$, and we have computed \f$\bar{x}^i\f$
231 : * by the above procedure. We then correct \f$\bar{x}^i\f$ by adding
232 : *
233 : * \f{align}
234 : * \delta \bar{x}^i = \left(x^j - x^j(\bar{x})\right)
235 : * \frac{\partial \bar{x}^i}{\partial x^j},
236 : * \f}
237 : *
238 : * where \f$x^j(\bar{x})\f$ is the result of applying the forward map
239 : * to \f$\bar{x}^i\f$ and \f$\partial \bar{x}^i/\partial x^j\f$ is the
240 : * inverse jacobian.
241 : *
242 : * ### Inverse jacobian
243 : *
244 : * We write the inverse Jacobian as
245 : *
246 : * \f{align}
247 : * \frac{\partial \bar{x}^i}{\partial x^j} =
248 : * \frac{\partial \bar{x}^i}{\partial x_0^k}
249 : * \frac{\partial x_0^k}{\partial x^j}
250 : * + \frac{\partial \bar{x}^i}{\partial \sigma}
251 : * \frac{\partial \sigma}{\partial x^j},
252 : * \f}
253 : *
254 : * where we have recognized that \f$\bar{x}^i\f$ depends both on
255 : * \f$x_0^k\f$ (the corresponding point on the 2-surface) and on
256 : * \f$\sigma\f$ (encoding the distance away from the 2-surface).
257 : *
258 : * We now evaluate Eq. (12). The `InnerMap` class provides a function
259 : * `inv_jacobian` that returns \f$\partial \bar{x}^i/\partial x_0^k\f$
260 : * (where \f$\sigma\f$ is held fixed), and a function `dxbar_dsigma`
261 : * that returns \f$\partial \bar{x}^i/\partial \sigma\f$ (where
262 : * \f$x_0^i\f$ is held fixed). The factor \f$\partial x_0^j/\partial
263 : * x^i\f$ can be computed by differentiating Eq. (7), which yields
264 : *
265 : * \f{align}
266 : * \frac{\partial x_0^j}{\partial x^i} &= \tilde{\lambda} \delta_i^j
267 : * + \frac{x_0^j-P^j}{\tilde{\lambda}}
268 : * \frac{\partial\tilde{\lambda}}{\partial x^i},
269 : * \f}
270 : *
271 : * where \f$\partial \tilde{\lambda}/\partial x^i\f$ is provided by
272 : * the `deriv_lambda_tilde` function of `InnerMap`. Note that for
273 : * nonsingular maps there is no worry that \f$\tilde{\lambda}\f$ is
274 : * zero in the denominator of the second term of Eq. (13); if
275 : * \f$\tilde{\lambda}=0\f$ then \f$x_0^i=P^i\f$ by Eq. (7), and therefore
276 : * the map is singular.
277 : *
278 : * To evaluate the remaining unknown factor in Eq. (12),
279 : * \f$\partial \sigma/\partial x^j\f$,
280 : * note that [combining Eqs. (1), (7), and (8)]
281 : * \f$\bar{\lambda}=\tilde{\lambda}\lambda\f$.
282 : * Therefore Eq. (10) is equivalent to
283 : *
284 : * \f{align}
285 : * \sigma &= \frac{\tilde{\lambda}-1}{\tilde{\lambda}(1-\lambda)}.
286 : * \f}
287 : *
288 : * Differentiating this expression yields
289 : *
290 : * \f{align}
291 : * \frac{\partial \sigma}{\partial x^i} &=
292 : * \frac{\partial \sigma}{\partial \lambda}
293 : * \frac{\partial \lambda}{\partial x_0^j}
294 : * \frac{\partial x_0^j}{\partial x^i}
295 : * + \frac{\partial \sigma}{\partial \tilde\lambda}
296 : * \frac{\partial \tilde\lambda}{\partial x^i}\\
297 : * &=
298 : * \frac{\sigma}{1-\lambda}
299 : * \frac{\partial \lambda}{\partial x_0^j}
300 : * \frac{\partial x_0^j}{\partial x^i}
301 : * +
302 : * \frac{1}{\tilde{\lambda}^2(1-\lambda)}
303 : * \frac{\partial \tilde\lambda}{\partial x^i},
304 : * \f}
305 : *
306 : * where the second factor in the first term can be evaluated using
307 : * Eq. (5), the third factor in the first term can be evaluated using
308 : * Eq. (13), and the second factor in the second term is provided by
309 : * `InnerMap`s function `deriv_lambda_tilde`.
310 : *
311 : */
312 : template <typename InnerMap>
313 1 : class FocallyLiftedMap {
314 : public:
315 0 : static constexpr size_t dim = 3;
316 0 : FocallyLiftedMap(const std::array<double, 3>& center,
317 : const std::array<double, 3>& proj_center, double radius,
318 : bool source_is_between_focus_and_target, InnerMap inner_map);
319 :
320 0 : FocallyLiftedMap() = default;
321 0 : ~FocallyLiftedMap() = default;
322 0 : FocallyLiftedMap(FocallyLiftedMap&&) = default;
323 0 : FocallyLiftedMap(const FocallyLiftedMap&) = default;
324 0 : FocallyLiftedMap& operator=(const FocallyLiftedMap&) = default;
325 0 : FocallyLiftedMap& operator=(FocallyLiftedMap&&) = default;
326 :
327 : template <typename T>
328 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
329 : const std::array<T, 3>& source_coords) const;
330 :
331 : /// The inverse function is only callable with doubles because the inverse
332 : /// might fail if called for a point out of range, and it is unclear
333 : /// what should happen if the inverse were to succeed for some points in a
334 : /// DataVector but fail for other points.
335 1 : std::optional<std::array<double, 3>> inverse(
336 : const std::array<double, 3>& target_coords) const;
337 :
338 : template <typename T>
339 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
340 : const std::array<T, 3>& source_coords) const;
341 :
342 : template <typename T>
343 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
344 : const std::array<T, 3>& source_coords) const;
345 :
346 : // NOLINTNEXTLINE(google-runtime-references)
347 0 : void pup(PUP::er& p);
348 :
349 0 : static bool is_identity() { return false; }
350 :
351 : private:
352 0 : friend bool operator==<InnerMap>(const FocallyLiftedMap<InnerMap>& lhs,
353 : const FocallyLiftedMap<InnerMap>& rhs);
354 0 : std::array<double, 3> center_{}, proj_center_{};
355 0 : double radius_{std::numeric_limits<double>::signaling_NaN()};
356 0 : bool source_is_between_focus_and_target_;
357 0 : InnerMap inner_map_;
358 : };
359 : template <typename InnerMap>
360 0 : bool operator!=(const FocallyLiftedMap<InnerMap>& lhs,
361 : const FocallyLiftedMap<InnerMap>& rhs);
362 : } // namespace domain::CoordinateMaps
|
