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