FocallyLiftedSide.hpp
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/MakeArray.hpp"
14 #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
15
16 /// \cond
17 namespace PUP {
18 class er;
19 } // namespace PUP
20 /// \endcond
21
23 /*!
24  * \brief A FocallyLiftedInnerMap that maps a 3D unit right cylindrical shell
25  * to a volume that connects portions of two spherical surfaces.
26  *
27  * \details The domain of the map is a 3D unit right cylinder with
28  * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
29  * \f$-1\leq\bar{z}\leq 1\f$ and \f$1\leq \bar{x}^2+\bar{y}^2 \leq 30 * 4\f$. The range of the map has coordinates \f$(x,y,z)\f$.
31  *
32  * Consider a sphere with center \f$C^i\f$ and radius \f$R\f$ that is
33  * intersected by two planes normal to the \f$z\f$ axis located at
34  * \f$z = z_\mathrm{L}\f$ and \f$z = z_\mathrm{U}\f$, with
35  * \f$z_\mathrm{L} < z_\mathrm{U}\f$.
36  * Side provides the following functions:
37  *
38  * ### forward_map()
39  * forward_map() maps \f$(\bar{x},\bar{y},\bar{z})\f$ to a point on the inner
40  * surface
41  * \f$\bar{x}^2+\bar{y}^2=1\f$ by dividing \f$\bar{x}\f$ and \f$\bar{y}\f$
42  * by \f$(1+\sigma)\f$, where \f$\sigma\f$ is the function given by Eq. (7)
43  * below.
44  * Then it maps that point to a point on the portion of the sphere with
45  * \f$z_\mathrm{L} \leq z \leq z_\mathrm{U}\f$.
46  * forward_map() returns
47  * \f$x_0^i\f$, the 3D coordinates on that sphere, which are given by
48  *
49  * \f{align}
50  * x_0^0 &= R \sin\theta \frac{\bar{x}}{1+\sigma} + C^0,\\
51  * x_0^1 &= R \sin\theta \frac{\bar{y}}{1+\sigma} + C^1,\\
52  * x_0^2 &= R \cos\theta + C^2.\\
53  * \f}
54  *
55  * Here
56  * \f{align}
57  * \theta = \theta_\mathrm{max} +
58  * (\theta_\mathrm{min}-\theta_\mathrm{max}) \frac{\bar{z}+1}{2},
59  * \f}
60  *
61  * where
62  * \f{align}
63  * \cos(\theta_\mathrm{max}) &= (z_\mathrm{L}-C^2)/R,\\
64  * \cos(\theta_\mathrm{min}) &= (z_\mathrm{U}-C^2)/R.
65  * \f}
66  *
67  * Note that \f$\theta\f$ decreases with increasing \f$\bar{z}\f$,
68  * which is the usual convention for a polar angle but might otherwise
69  * cause confusion.
70  *
71  * ### sigma
72  *
73  * \f$\sigma\f$ is a function that is zero on the sphere
74  * \f$x^i=x_0^i\f$ and unity at \f$\bar{x}^2+\bar{y}^2=4\f$
75  * (corresponding to the upper surface of the FocallyLiftedMap). We define
76  *
77  * \f{align}
78  * \sigma &= \sqrt{\bar{x}^2+\bar{y}^2}-1.
79  * \f}
80  *
81  * ### deriv_sigma
82  *
83  * deriv_sigma returns
84  *
85  * \f{align}
86  * \frac{\partial \sigma}{\partial \bar{x}^j} &=
87  * \left(\frac{\bar{x}}{1+\sigma},
88  * \frac{\bar{y}}{1+\sigma},0\right).
89  * \f}
90  *
91  * ### jacobian
92  *
93  * jacobian returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
94  * The arguments to jacobian are \f$(\bar{x},\bar{y},\bar{z})\f$.
95  * Differentiating Eqs.(1--4) above yields
96  *
97  * \f{align*}
98  * \frac{\partial x_0^0}{\partial \bar{x}} &= R \sin\theta
99  * \frac{\bar{y}^2}{(1+\sigma)^3}, \\
100  * \frac{\partial x_0^0}{\partial \bar{y}} &= -R \sin\theta
101  * \frac{\bar{x}\bar{y}}{(1+\sigma)^3}, \\
102  * \frac{\partial x_0^0}{\partial \bar{z}} &=
103  * R \cos\theta \frac{\theta_\mathrm{min}-\theta_\mathrm{max}}{2(1+\sigma)}
104  * \bar{x},\\
105  * \frac{\partial x_0^1}{\partial \bar{x}} &= -R \sin\theta
106  * \frac{\bar{x}\bar{y}}{(1+\sigma)^3}, \\
107  * \frac{\partial x_0^1}{\partial \bar{y}} &= R \sin\theta
108  * \frac{\bar{x}^2}{(1+\sigma)^3}, \\
109  * \frac{\partial x_0^1}{\partial \bar{z}} &=
110  * R \cos\theta \frac{\theta_\mathrm{min}-\theta_\mathrm{max}}{2(1+\sigma)}
111  * \bar{y},\\
112  * \frac{\partial x_0^2}{\partial \bar{x}} &= 0,\\
113  * \frac{\partial x_0^2}{\partial \bar{y}} &= 0,\\
114  * \frac{\partial x_0^2}{\partial \bar{z}} &=
115  * - R \sin\theta \frac{\theta_\mathrm{min}-\theta_\mathrm{max}}{2}.
116  * \f}
117  *
118  * ### inverse
119  *
120  * inverse takes \f$x_0^i\f$ and \f$\sigma\f$ as arguments, and
121  * returns \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed
122  * std::optional<std::array<double, 3>> if \f$x_0^i\f$ or \f$\sigma\f$ are
123  * outside the range of the map.
124  *
125  * If \f$\sigma\f$ is outside the range \f$[0,1]\f$ then we return
126  * a default-constructed std::optional<std::array<double, 3>>.
127  *
128  * To get \f$\bar{z}\f$ we invert Eq. (4):
129  * \f{align}
130  * \bar{z} &= 2\frac{\acos\left((x_0^2-C^2)/R\right)-\theta_\mathrm{max}}
131  * {\theta_\mathrm{min}-\theta_\mathrm{max}} - 1.
132  * \f}
133  *
134  * If \f$\bar{z}\f$ is outside the range \f$[-1,1]\f$ then we return
135  * a default-constructed std::optional<std::array<double, 3>>.
136  *
137  * To compute \f$\bar{x}\f$ and \f$\bar{y}\f$, we invert Eqs. (1--3) and
138  * use \f$\sigma\f$:
139  *
140  * \f{align}
141  * \bar{x} &= \frac{(x_0^0-C^0) (1+\sigma)}{\rho},\\
142  * \bar{y} &= \frac{(x_0^1-C^1) (1+\sigma)}{\rho},
143  * \f}
144  *
145  * where
146  *
147  * \f{align}
148  * \rho = \sqrt{(x_0^0-C^0)^2+(x_0^1-C^1)^2}.
149  * \f}
150  *
151  * ### lambda_tilde
152  *
153  * lambda_tilde takes as arguments a point \f$x^i\f$ and a projection point
154  * \f$P^i\f$, and computes \f$\tilde{\lambda}\f$, the solution to
155  *
156  * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
157  *
158  * Since \f$x_0^i\f$ must lie on the sphere, \f$\tilde{\lambda}\f$ is the
159  * solution of the quadratic equation
160  *
161  * \f{align}
162  * |P^i + (x^i - P^i) \tilde{\lambda} - C^i |^2 - R^2 = 0.
163  * \f}
164  *
165  * In solving the quadratic, we choose the larger root if
166  * \f$x^2>z_\mathrm{P}\f$ and the smaller root otherwise. We demand
167  * that the root is greater than unity. If there is no such root,
168  * this means that the point \f$x^i\f$ is not in the range of the map
169  * so we return a default-constructed std::optional<double>.
170  *
171  * ### deriv_lambda_tilde
172  *
173  * deriv_lambda_tilde takes as arguments \f$x_0^i\f$, a projection point
174  * \f$P^i\f$, and \f$\tilde{\lambda}\f$, and
175  * returns \f$\partial \tilde{\lambda}/\partial x^i\f$.
176  * By differentiating Eq. (14), we find
177  *
178  * \f{align}
179  * \frac{\partial\tilde{\lambda}}{\partial x^j} &=
180  * \tilde{\lambda}^2 \frac{C^j - x_0^j}{
181  * (x_0^i - P^i)(x_{0i} - C_{i})} \nonumber \\
182  * &= \tilde{\lambda}^2 \frac{C^j - x_0^j}{|x_0^i - P^i|^2
183  * + (x_0^i - P^i)(P_i - C_{i})}.
184  * \f}
185  *
186  * ### inv_jacobian
187  *
188  * inv_jacobian returns \f$\partial \bar{x}^i/\partial x_0^k\f$,
189  * where \f$\sigma\f$ is held fixed.
190  * The arguments to inv_jacobian are \f$(\bar{x},\bar{y},\bar{z})\f$.
191  *
192  * Note from Eqs. (9--12) that \f$\bar{x}\f$ and \f$\bar{y}\f$
193  * depend only on \f$x_0^0\f$ and \f$x_0^1\f$ but not on \f$x_0^2\f$.
194  *
195  * By differentiating Eqs. (9--12), we find
196  *
197  * \f{align*}
198  * \frac{\partial \bar{x}}{\partial x_0^0} &=
199  * \frac{\bar{y}^2}{(1+\sigma)\rho},\\
200  * \frac{\partial \bar{x}}{\partial x_0^1} &=
201  * - \frac{\bar{x}\bar{y}}{(1+\sigma)\rho},\\
202  * \frac{\partial \bar{x}}{\partial x_0^2} &= 0,\\
203  * \frac{\partial \bar{y}}{\partial x_0^0} &=
204  * - \frac{\bar{x}\bar{y}}{(1+\sigma)\rho},\\
205  * \frac{\partial \bar{y}}{\partial x_0^1} &=
206  * \frac{\bar{x}^2}{(1+\sigma)\rho},\\
207  * \frac{\partial \bar{y}}{\partial x_0^2} &= 0,\\
208  * \frac{\partial \bar{z}}{\partial x_0^0} &= 0,\\
209  * \frac{\partial \bar{z}}{\partial x_0^1} &= 0,\\
210  * \frac{\partial \bar{z}}{\partial x_0^2} &=
211  * -\frac{2}{\rho(\theta_\mathrm{min}-\theta_\mathrm{max})},
212  * \f}
213  *
214  * where
215  *
216  * \f[
217  * \rho = R \sin\theta = R\sin\left(\theta_\mathrm{max} +
218  * (\theta_\mathrm{min}-\theta_\mathrm{max}) \frac{\bar{z}+1}{2}\right),
219  * \f]
220  *
221  * which is also equal to the quantity in Eq. (12).
222  *
223  * ### dxbar_dsigma
224  *
225  * dxbar_dsigma returns \f$\partial \bar{x}^i/\partial \sigma\f$,
226  * where \f$x_0^i\f$ is held fixed.
227  *
228  * From Eqs. (10) and (11) we have
229  *
230  * \f{align}
231  * \frac{\partial \bar{x}^i}{\partial \sigma} &=
232  * \left(\frac{\bar{x}}{\sqrt{\bar{x}^2+\bar{y}^2}},
233  * \frac{\bar{y}}{\sqrt{\bar{x}^2+\bar{y}^2}},0\right).
234  * \f}
235  *
236  */
237 class Side {
238  public:
239  static constexpr size_t dim = 3;
240  Side(const std::array<double, 3>& center, const double radius,
241  const double z_lower, const double z_upper) noexcept;
242
243  Side() = default;
244  ~Side() = default;
245  Side(Side&&) = default;
246  Side(const Side&) = default;
247  Side& operator=(const Side&) = default;
248  Side& operator=(Side&&) = default;
249
250  template <typename T>
251  void forward_map(
252  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
253  target_coords,
254  const std::array<T, 3>& source_coords) const noexcept;
255
257  const std::array<double, 3>& target_coords,
258  double sigma_in) const noexcept;
259
260  template <typename T>
261  void jacobian(const gsl::not_null<
262  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
263  jacobian_out,
264  const std::array<T, 3>& source_coords) const noexcept;
265
266  template <typename T>
267  void inv_jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3,
268  Frame::NoFrame>*>
269  inv_jacobian_out,
270  const std::array<T, 3>& source_coords) const noexcept;
271
272  template <typename T>
273  void sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,
274  const std::array<T, 3>& source_coords) const noexcept;
275
276  template <typename T>
277  void deriv_sigma(
278  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
279  deriv_sigma_out,
280  const std::array<T, 3>& source_coords) const noexcept;
281
282  template <typename T>
283  void dxbar_dsigma(
284  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
285  dxbar_dsigma_out,
286  const std::array<T, 3>& source_coords) const noexcept;
287
288  std::optional<double> lambda_tilde(
289  const std::array<double, 3>& parent_mapped_target_coords,
290  const std::array<double, 3>& projection_point,
291  bool source_is_between_focus_and_target) const noexcept;
292
293  template <typename T>
294  void deriv_lambda_tilde(
295  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
296  deriv_lambda_tilde_out,
297  const std::array<T, 3>& target_coords, const T& lambda_tilde,
298  const std::array<double, 3>& projection_point) const noexcept;
299
300  // clang-tidy: google runtime references
301  void pup(PUP::er& p) noexcept; // NOLINT
302
303  static bool is_identity() noexcept { return false; }
304
305  private:
306  friend bool operator==(const Side& lhs, const Side& rhs) noexcept;
307  std::array<double, 3> center_{
310  double theta_min_{std::numeric_limits<double>::signaling_NaN()};
311  double theta_max_{std::numeric_limits<double>::signaling_NaN()};
312 };
313 bool operator!=(const Side& lhs, const Side& rhs) noexcept;
314 } // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps
MakeArray.hpp
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
domain::CoordinateMaps::FocallyLiftedInnerMaps
Contains FocallyLiftedInnerMaps.
Definition: FocallyLiftedEndcap.hpp:22
cstddef
array
domain::CoordinateMaps::FocallyLiftedInnerMaps::Side
A FocallyLiftedInnerMap that maps a 3D unit right cylindrical shell to a volume that connects portion...
Definition: FocallyLiftedSide.hpp:237
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
limits
Gsl.hpp
TypeAliases.hpp
optional
std::numeric_limits
gsl::not_null
Require a pointer to not be a nullptr