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 
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>
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  * ### Jacobian
103  *
104  * Differentiating Eq. (1) above yields
105  *
106  * \f{align}
107  * \frac{\partial x_1^i}{\partial x_0^j} = \lambda \delta^i_j +
108  * (x_0^i - P^i) \frac{\partial \lambda}{\partial x_0^j}.
109  * \f}
110  *
111  * and differentiating Eq. (2) and then inserting Eq. (1) yields
112  *
113  * \f{align}
114  * \frac{\partial\lambda}{\partial x_0^j} &=
115  * \lambda^2 \frac{C_j - x_1^j}{|x_1^i - P^i|^2
116  * + (x_1^i - P^i)(P_i - C_{i})}.
117  * \f}
118  *
119  * The Jacobian can be found by differentiating Eq. (3) above and using the
120  * chain rule, recognizing that \f$x_1^i\f$ depends on \f$\bar{x}^i\f$ only
121  * through its dependence on \f$x_0^i\f$; this is because there is no explicit
122  * dependence on \f$\bar{x}^i\f$ in Eq. (1) (which determines \f$x_1^i\f$)
123  * or Eq. (2) (which determines \f$\lambda\f$).
124  *
125  * \f{align}
126  * \frac{\partial x^i}{\partial \bar{x}^j} &=
127  * \sigma \frac{\partial x_1^i}{\partial x_0^k}
128  * \frac{\partial x_0^k}{\partial \bar{x}^j} +
129  * (1-\sigma)\frac{\partial x_0^i}{\partial \bar{x}^j}
130  * + \frac{\partial \sigma}{\partial \bar{x}^j} (x_1^i - x_0^i),
131  * \nonumber \\
132  * &= (1-\sigma+\lambda\sigma) \frac{\partial x_0^i}{\partial \bar{x}^j} +
133  * \sigma (x_0^i - P^i) \frac{\partial \lambda}{\partial x_0^k}
134  * \frac{\partial x_0^k}{\partial \bar{x}^j}
135  * + \frac{\partial \sigma}{\partial \bar{x}^j} (x_1^i - x_0^i),
136  * \f}
137  * where in the last line we have substituted Eq. (4).
138  *
139  * The class defined by the `InnerMap` template parameter
140  * provides the function `deriv_sigma`, which returns
141  * \f$\partial \sigma/\partial \bar{x}^j\f$, and the function `jacobian`,
142  * which returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
143  *
144  * ### Inverse map.
145  *
146  * Given \f$x^i\f$, we wish to compute \f$\bar{x}^i\f$.
147  *
148  * We first find the coordinates \f$x_0^i\f$ that lie on the 2-surface
149  * and are defined such that \f$P^i\f$, \f$x_0^i\f$,
150  * and \f$x^i\f$ are co-linear.
151  * See the right panel of the above figure.
152  * \f$x_0^i\f$ is determined by the equation
153  *
154  * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda},\f}
155  *
156  * where \f$\tilde{\lambda}\f$ is a scalar factor that depends on
157  * \f$x^i\f$ and is determined by the class defined by the `InnerMap` template
158  * parameter. `InnerMap`
159  * provides a function `lambda_tilde` that takes \f$x^i\f$ and
160  * \f$P^i\f$ as arguments and returns \f$\tilde{\lambda}\f$ (or
161  * a default-constructed `std::optional` if the appropriate
162  * \f$\tilde{\lambda}\f$ cannot be
163  * found; this default-constructed value indicates that the point \f$x^i\f$ is
164  * outside the range of the map).
165  *
166  * Now consider the coordinates \f$x_1^i\f$ that lie on the sphere and
167  * are defined such that \f$P^i\f$, \f$x_1^i\f$, and \f$x^i\f$ are
168  * co-linear. See the right panel of the figure.
169  * \f$x_1^i\f$ is determined by the equation
170  *
171  * \f{align} x_1^i = P^i + (x^i - P^i) \bar{\lambda},\f}
172  *
173  * where \f$\bar{\lambda}\f$ is a scalar factor that depends on \f$x^i\f$ and
174  * is the solution of a quadratic
175  * that is derived by demanding that \f$x_1^i\f$ lies on the sphere:
176  *
177  * \f{align}
178  * |P^i + (x^i - P^i) \bar{\lambda} - C^i |^2 - R^2 = 0.
179  * \f}
180  *
181  * Eq. (9) is a quadratic equation that
182  * takes the usual form \f$a\bar{\lambda}^2+b\bar{\lambda}+c=0\f$,
183  * with
184  *
185  * \f{align*}
186  * a &= |x^i-P^i|^2,\\
187  * b &= 2(x^i-P^i)(P^j-C^j)\delta_{ij},\\
188  * c &= |P^i-C^i|^2 - R^2.
189  * \f}
190  *
191  * Note that we don't actually need to compute \f$x_1^i\f$. Instead, we
192  * can determine \f$\sigma\f$ by the relation (obtained by solving Eq. (3)
193  * for \f$\sigma\f$ and then inserting Eqs. (7) and (8) to eliminate
194  * \f$x_1^i\f$ and \f$x_0^i\f$)
195  *
196  * \f{align}
197  * \sigma = \frac{\tilde{\lambda}-1}{\tilde{\lambda}-\bar{\lambda}}.
198  * \f}
199  *
200  * The denominator of Eq. (10) is nonzero for nonsingular maps:
201  * From Eqs. (7) and (8), \f$\bar{\lambda}=\tilde{\lambda}\f$ means
202  * that \f$x_1^i=x_0^i\f$, which means that
203  * the \f$\sigma=0\f$ and \f$\sigma=1\f$ surfaces intersect, i.e.
204  * the map is singular.
205  *
206  * Once we have \f$x_0^i\f$ and \f$\sigma\f$, the point
207  * \f$(\bar{x},\bar{y},\bar{z})\f$ is uniquely determined by `InnerMap`.
208  * The `inverse` function of `InnerMap` takes \f$x_0^i\f$ and \f$\sigma\f$
209  * as arguments, and returns
210  * \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed `std::optional`
211  * if \f$x_0^i\f$ or
212  * \f$\sigma\f$ are outside the range of the map.
213  *
214  * #### Root polishing
215  *
216  * The inverse function described above will sometimes have errors that
217  * are noticeably larger than roundoff. Therefore we apply a single
218  * Newton-Raphson iteration to refine the result of the inverse map:
219  * Suppose we are given \f$x^i\f$, and we have computed \f$\bar{x}^i\f$
220  * by the above procedure. We then correct \f$\bar{x}^i\f$ by adding
221  *
222  * \f{align}
223  * \delta \bar{x}^i = \left(x^j - x^j(\bar{x})\right)
224  * \frac{\partial \bar{x}^i}{\partial x^j},
225  * \f}
226  *
227  * where \f$x^j(\bar{x})\f$ is the result of applying the forward map
228  * to \f$\bar{x}^i\f$ and \f$\partial \bar{x}^i/\partial x^j\f$ is the
229  * inverse jacobian.
230  *
231  * ### Inverse jacobian
232  *
233  * We write the inverse Jacobian as
234  *
235  * \f{align}
236  * \frac{\partial \bar{x}^i}{\partial x^j} =
237  * \frac{\partial \bar{x}^i}{\partial x_0^k}
238  * \frac{\partial x_0^k}{\partial x^j}
239  * + \frac{\partial \bar{x}^i}{\partial \sigma}
240  * \frac{\partial \sigma}{\partial x^j},
241  * \f}
242  *
243  * where we have recognized that \f$\bar{x}^i\f$ depends both on
244  * \f$x_0^k\f$ (the corresponding point on the 2-surface) and on
245  * \f$\sigma\f$ (encoding the distance away from the 2-surface).
246  *
247  * We now evaluate Eq. (12). The `InnerMap` class provides a function
248  * `inv_jacobian` that returns \f$\partial \bar{x}^i/\partial x_0^k\f$
249  * (where \f$\sigma\f$ is held fixed), and a function `dxbar_dsigma`
250  * that returns \f$\partial \bar{x}^i/\partial \sigma\f$ (where
251  * \f$x_0^i\f$ is held fixed). The factor \f$\partial x_0^j/\partial
252  * x^i\f$ can be computed by differentiating Eq. (7), which yields
253  *
254  * \f{align}
255  * \frac{\partial x_0^j}{\partial x^i} &= \tilde{\lambda} \delta_i^j
256  * + \frac{x_0^j-P^j}{\tilde{\lambda}}
257  * \frac{\partial\tilde{\lambda}}{\partial x^i},
258  * \f}
259  *
260  * where \f$\partial \tilde{\lambda}/\partial x^i\f$ is provided by
261  * the `deriv_lambda_tilde` function of `InnerMap`. Note that for
262  * nonsingular maps there is no worry that \f$\tilde{\lambda}\f$ is
263  * zero in the denominator of the second term of Eq. (13); if
264  * \f$\tilde{\lambda}=0\f$ then \f$x_0^i=P^i\f$ by Eq. (7), and therefore
265  * the map is singular.
266  *
267  * To evaluate the remaining unknown factor in Eq. (12),
268  * \f$\partial \sigma/\partial x^j\f$,
269  * note that [combining Eqs. (1), (7), and (8)]
270  * \f$\bar{\lambda}=\tilde{\lambda}\lambda\f$.
271  * Therefore Eq. (10) is equivalent to
272  *
273  * \f{align}
274  * \sigma &= \frac{\tilde{\lambda}-1}{\tilde{\lambda}(1-\lambda)}.
275  * \f}
276  *
277  * Differentiating this expression yields
278  *
279  * \f{align}
280  * \frac{\partial \sigma}{\partial x^i} &=
281  * \frac{\partial \sigma}{\partial \lambda}
282  * \frac{\partial \lambda}{\partial x_0^j}
283  * \frac{\partial x_0^j}{\partial x^i}
284  * + \frac{\partial \sigma}{\partial \tilde\lambda}
285  * \frac{\partial \tilde\lambda}{\partial x^i}\\
286  * &=
287  * \frac{\sigma}{1-\lambda}
288  * \frac{\partial \lambda}{\partial x_0^j}
289  * \frac{\partial x_0^j}{\partial x^i}
290  * +
291  * \frac{1}{\tilde{\lambda}^2(1-\lambda)}
292  * \frac{\partial \tilde\lambda}{\partial x^i},
293  * \f}
294  *
295  * where the second factor in the first term can be evaluated using
296  * Eq. (5), the third factor in the first term can be evaluated using
297  * Eq. (13), and the second factor in the second term is provided by
298  * `InnerMap`s function `deriv_lambda_tilde`.
299  *
300  */
301 template <typename InnerMap>
302 class FocallyLiftedMap {
303  public:
304  static constexpr size_t dim = 3;
306  const std::array<double, 3>& proj_center, double radius,
307  InnerMap inner_map) noexcept;
308 
309  FocallyLiftedMap() = default;
310  ~FocallyLiftedMap() = default;
311  FocallyLiftedMap(FocallyLiftedMap&&) = default;
312  FocallyLiftedMap(const FocallyLiftedMap&) = default;
313  FocallyLiftedMap& operator=(const FocallyLiftedMap&) = default;
314  FocallyLiftedMap& operator=(FocallyLiftedMap&&) = default;
315 
316  template <typename T>
318  const std::array<T, 3>& source_coords) const noexcept;
319 
321  const std::array<double, 3>& target_coords) const noexcept;
322 
323  template <typename T>
324  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
325  const std::array<T, 3>& source_coords) const noexcept;
326 
327  template <typename T>
328  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
329  const std::array<T, 3>& source_coords) const noexcept;
330 
331  // clang-tidy: google runtime references
332  void pup(PUP::er& p) noexcept; // NOLINT
333 
334  static bool is_identity() noexcept { return false; }
335 
336  private:
337  friend bool operator==
338  <InnerMap>(const FocallyLiftedMap<InnerMap>& lhs,
339  const FocallyLiftedMap<InnerMap>& rhs) noexcept;
340  std::array<double, 3> center_{}, proj_center_{};
342  InnerMap inner_map_;
343 };
344 template <typename InnerMap>
345 bool operator!=(const FocallyLiftedMap<InnerMap>& lhs,
346  const FocallyLiftedMap<InnerMap>& rhs) noexcept;
347 } // 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.cpp:14
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
limits
TypeAliases.hpp
optional