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