FocallyLiftedMapHelpers.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/Gsl.hpp"
13 #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
14 
15 /// Holds helper functions for use with
16 /// domain::CoordinateMaps::FocallyLiftedMap.
18 
19 /*!
20  * \brief Finds how long to extend a line segment to have it intersect
21  * a point on a 2-sphere.
22  *
23  * \details Consider a 2-sphere with center \f$C^i\f$ and radius \f$R\f$, and
24  * and let \f$P^i\f$ and \f$x_0^i\f$ be two arbitrary 3D points.
25  *
26  * Consider the line passing through \f$P^i\f$ and \f$x_0^i\f$.
27  * If this line intersects the sphere at a point \f$x_1^i\f$, then we can write
28  *
29  * \f[
30  * x_1^i = P^i + (x_0^i-P^i) \lambda,
31  * \f]
32  *
33  * where \f$\lambda\f$ is a scale factor.
34  *
35  * `scale_factor` computes and returns \f$\lambda\f$.
36  *
37  * ### Even more detail:
38  *
39  * To solve for \f$\lambda\f$, we note that \f$x_1^i\f$ is on the surface of
40  * the sphere, so
41  *
42  * \f[
43  * |x_1^i-C^i|^2 = R^2,
44  * \f]
45  *
46  * (where \f$|A^i|^2\f$ means \f$\delta_{ij} A^i A^j\f$),
47  *
48  * or equivalently
49  *
50  * \f[
51  * | P^i-C^i + (x_0^i-P^i)\lambda |^2 = R^2.
52  * \f]
53  *
54  * This is a quadratic equation for \f$\lambda\f$
55  * and it generally has more than one real root.
56  * It takes the usual form \f$a\lambda^2+b\lambda+c=0\f$,
57  * with
58  *
59  * \f{align*}
60  * a &= |x_0^i-P^i|^2,\\
61  * b &= 2(x_0^i-P^i)(P^j-C^j)\delta_{ij},\\
62  * c &= |P^i-C^i|^2 - R^2,
63  * \f}
64  *
65  * So how do we choose between multiple roots? Some of the maps that
66  * use `scale_factor` assume that *for all points*, \f$x_0^i\f$ is
67  * between \f$P^i\f$ and \f$x_1^i\f$. Those maps should set the parameter
68  * `src_is_between_proj_and_target` to true. Other maps assume that
69  * *for all points*, \f$x^i\f$ is always between \f$x_0^i\f$
70  * and \f$P^i\f$. Those maps should set the parameter
71  * `src_is_between_proj_and_target` to false.
72  *
73  * \warning If we ever add maps where
74  * `src_is_between_proj_and_target` can change from point to point,
75  * the logic of `scale_factor` needs to be changed.
76  *
77  * In the arguments to the function below, `src_point` is \f$x_0^i\f$,
78  * `proj_center` is \f$P^i\f$, `sphere_center` is \f$C^i\f$, and
79  * `radius` is \f$R\f$.
80  *
81  */
82 template <typename T>
83 void scale_factor(const gsl::not_null<tt::remove_cvref_wrap_t<T>*>& result,
84  const std::array<T, 3>& src_point,
85  const std::array<double, 3>& proj_center,
86  const std::array<double, 3>& sphere_center, double radius,
87  bool src_is_between_proj_and_target) noexcept;
88 
89 /*!
90  * Solves a problem of the same form as `scale_factor`, but is used
91  * only by the inverse function to compute \f$\tilde{\lambda}\f$ and
92  * \f$\bar{\lambda}\f$. `try_scale_factor` is used in two contexts:
93  *
94  * `try_scale_factor` is used to determine \f$\bar{\lambda}\f$
95  * given \f$x^i\f$. \f$\bar{\lambda}\f$ is defined by
96  * \f{align*} x_1^i = P^i + (x^i - P^i) \bar{\lambda}.\f}
97  *
98  * `try_scale_factor` is used by the `lambda_tilde` functions of some
99  * of the `InnerMap` classes (namely those `InnerMap` classes where
100  * \f$x_0^i\f$ is a spherical surface) to solve for
101  * \f$\tilde{\lambda}\f$ given\f$x^i\f$. \f$\tilde{\lambda}\f$
102  * is defined by
103  * \f{align*} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
104  *
105  * In both of these contexts, the input parameter `src_point` is
106  * \f$x^i\f$, a point that is supposed to be in the range of the
107  * `FocallyLiftedMap`. Because the inverse function can be and is
108  * called for an arbitrary \f$x^i\f$ that might not be in the range
109  * of the `FocallyLiftedMap`, `try_scale_factor` returns a
110  * std::optional, with a default-constructed std::optional if the roots it
111  * finds are not as expected (i.e. if the inverse map was called for
112  * a point not in the range of the map).
113  *
114  * Because `try_scale_factor` can be called in different situations,
115  * it has additional boolean arguments `pick_larger_root` and
116  * `pick_root_greater_than_one` that allow the caller to choose which
117  * root to return.
118  *
119  * `try_scale_factor` is not templated
120  * on type because it is used only by the inverse function, which
121  * works only on doubles.
122  *
123  */
125  const std::array<double, 3>& src_point,
126  const std::array<double, 3>& proj_center,
127  const std::array<double, 3>& sphere_center, double radius,
128  bool pick_larger_root, bool pick_root_greater_than_one) noexcept;
129 
130 /*!
131  * Computes \f$\partial \lambda/\partial x_0^i\f$, where \f$\lambda\f$
132  * is the quantity returned by `scale_factor` and `x_0` is `src_point` in
133  * the `scale_factor` function.
134  *
135  * The formula (see `FocallyLiftedMap`) is
136  * \f{align*}
137  * \frac{\partial\lambda}{\partial x_0^j} &=
138  * \lambda^2 \frac{C_j - x_1^j}{|x_1^i - P^i|^2
139  * + (x_1^i - P^i)(P_i - C_i)}.
140  * \f}
141  *
142  * Note that it takes `intersection_point` and not `src_point` as a
143  * parameter.
144  *
145  * In the arguments to the function below, `intersection_point` is \f$x_1\f$,
146  * `proj_center` is \f$P^i\f$, `sphere_center` is \f$C^i\f$, and
147  * `radius` is \f$R\f$.
148  */
149 template <typename T>
151  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>& result,
152  const std::array<T, 3>& intersection_point,
153  const std::array<double, 3>& proj_center,
154  const std::array<double, 3>& sphere_center, const T& lambda) noexcept;
155 
156 } // namespace domain::CoordinateMaps::FocallyLiftedMapHelpers
domain::CoordinateMaps::FocallyLiftedMapHelpers::d_scale_factor_d_src_point
void d_scale_factor_d_src_point(const gsl::not_null< std::array< tt::remove_cvref_wrap_t< T >, 3 > * > &result, const std::array< T, 3 > &intersection_point, const std::array< double, 3 > &proj_center, const std::array< double, 3 > &sphere_center, const T &lambda) noexcept
cstddef
array
domain::CoordinateMaps::FocallyLiftedMapHelpers::try_scale_factor
std::optional< double > try_scale_factor(const std::array< double, 3 > &src_point, const std::array< double, 3 > &proj_center, const std::array< double, 3 > &sphere_center, double radius, bool pick_larger_root, bool pick_root_greater_than_one) noexcept
domain::CoordinateMaps::FocallyLiftedMapHelpers::scale_factor
void scale_factor(const gsl::not_null< tt::remove_cvref_wrap_t< T > * > &result, const std::array< T, 3 > &src_point, const std::array< double, 3 > &proj_center, const std::array< double, 3 > &sphere_center, double radius, bool src_is_between_proj_and_target) noexcept
Finds how long to extend a line segment to have it intersect a point on a 2-sphere.
limits
Gsl.hpp
TypeAliases.hpp
domain::CoordinateMaps::FocallyLiftedMapHelpers
Holds helper functions for use with domain::CoordinateMaps::FocallyLiftedMap.
Definition: FocallyLiftedMapHelpers.hpp:17
optional
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13