FocallyLiftedMap.hpp
1 // 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 bool operator==(const FocallyLiftedMap<InnerMap>& lhs,
27  const FocallyLiftedMap<InnerMap>& rhs) noexcept;
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 class FocallyLiftedMap {
314  public:
315  static constexpr size_t dim = 3;
316  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,
319  InnerMap inner_map) noexcept;
320 
321  FocallyLiftedMap() = default;
322  ~FocallyLiftedMap() = default;
323  FocallyLiftedMap(FocallyLiftedMap&&) = default;
324  FocallyLiftedMap(const FocallyLiftedMap&) = default;
325  FocallyLiftedMap& operator=(const FocallyLiftedMap&) = default;
326  FocallyLiftedMap& operator=(FocallyLiftedMap&&) = default;
327 
328  template <typename T>
329  std::array<tt::remove_cvref_wrap_t<T>, 3> operator()(
330  const std::array<T, 3>& source_coords) const noexcept;
331 
332  /// The inverse function is only callable with doubles because the inverse
333  /// might fail if called for a point out of range, and it is unclear
334  /// what should happen if the inverse were to succeed for some points in a
335  /// DataVector but fail for other points.
336  std::optional<std::array<double, 3>> inverse(
337  const std::array<double, 3>& target_coords) const noexcept;
338 
339  template <typename T>
340  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
341  const std::array<T, 3>& source_coords) const noexcept;
342 
343  template <typename T>
344  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
345  const std::array<T, 3>& source_coords) const noexcept;
346 
347  // clang-tidy: google runtime references
348  void pup(PUP::er& p) noexcept; // NOLINT
349 
350  static bool is_identity() noexcept { return false; }
351 
352  private:
353  friend bool operator==
354  <InnerMap>(const FocallyLiftedMap<InnerMap>& lhs,
355  const FocallyLiftedMap<InnerMap>& rhs) noexcept;
356  std::array<double, 3> center_{}, proj_center_{};
357  double radius_{std::numeric_limits<double>::signaling_NaN()};
358  bool source_is_between_focus_and_target_;
359  InnerMap inner_map_;
360 };
361 template <typename InnerMap>
362 bool operator!=(const FocallyLiftedMap<InnerMap>& lhs,
363  const FocallyLiftedMap<InnerMap>& rhs) noexcept;
364 } // namespace domain::CoordinateMaps
std::rel_ops::operator!=
T operator!=(T... args)
domain::CoordinateMaps::FocallyLiftedMap
Map from to the volume contained between a 2D surface and the surface of a 2-sphere.
Definition: FocallyLiftedMap.hpp:23
Frame::NoFrame
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
cstddef
array
domain::CoordinateMaps
Contains all coordinate maps.
Definition: Affine.hpp:23
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
domain::CoordinateMaps::FocallyLiftedMap::inverse
std::optional< std::array< double, 3 > > inverse(const std::array< double, 3 > &target_coords) const noexcept
The inverse function is only callable with doubles because the inverse might fail if called for a poi...
limits
TypeAliases.hpp
optional
© Copyright 2017 - 2021 SXS Collaboration, MIT License