FocallyLiftedFlatSide.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 /// \cond
16 namespace PUP {
17 class er;
18 } // namespace PUP
19 /// \endcond
20 
21 /// Contains FocallyLiftedInnerMaps
23 /*!
24  * \brief A FocallyLiftedInnerMap that maps a 3D unit right cylinder
25  * to a volume that connects a 2D annulus to a spherical
26  * surface.
27  *
28  * \details The domain of the map is a 3D unit right cylinder with
29  * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
30  * \f$-1\leq\bar{z}\leq 1\f$ and \f$1\leq \bar{x}^2+\bar{y}^2 \leq
31  * 4\f$. The range of the map has coordinates \f$(x,y,z)\f$.
32  *
33  * Consider a 2D annulus in 3D space
34  * oriented normal to the \f$z\f$ axis. The inner and outer radii
35  * of the annulus are \f$R_\mathrm{in}\f$ and \f$R_\mathrm{out}\f$, and the
36  * (3D) center of the annulus is \f$C^i\f$.
37  * `FlatSide` provides the following functions:
38  *
39  * ### forward_map()
40  * `forward_map()` maps \f$(\bar{x},\bar{y},\bar{z}=-1)\f$ to the interior
41  * of the annulus. The arguments to `forward_map()`
42  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
43  * `forward_map()` returns \f$x_0^i\f$,
44  * the 3D coordinates on the annulus, which are given by
45  *
46  * \f{align}
47  * x_0^0 &= \left(R_\mathrm{in}+(R_\mathrm{out}-R_\mathrm{in})
48  * (\bar{\rho}-1)\right)
49  * \frac{\bar{x}}{\bar{\rho}} + C^0, \label{eq:forward_map_x}\\
50  * x_0^1 &= \left(R_\mathrm{in}+(R_\mathrm{out}-R_\mathrm{in})
51  * (\bar{\rho}-1)\right)
52  * \frac{\bar{y}}{\bar{\rho}} + C^1,\\
53  * x_0^2 &= C^2 \label{eq:forward_map_z},
54  * \f}
55  *
56  * where
57  *
58  * \f{align} \bar{\rho} = \sqrt{\bar{x}^2+\bar{y}^2}.\label{eq:rhobar}\f}
59  *
60  * ### sigma
61  *
62  * \f$\sigma\f$ is a function that is zero on the sphere
63  * \f$x^i=x_0^i\f$ and unity at \f$\bar{z}=+1\f$ (corresponding to the
64  * upper surface of the FocallyLiftedMap). We define
65  *
66  * \f{align}
67  * \sigma &= \frac{\bar{z}+1}{2}.
68  * \f}
69  *
70  * ### deriv_sigma
71  *
72  * `deriv_sigma` returns
73  *
74  * \f{align}
75  * \frac{\partial \sigma}{\partial \bar{x}^j} &= (0,0,1/2).
76  * \label{eq:deriv_sigma}
77  * \f}
78  *
79  * ### jacobian
80  *
81  * `jacobian` returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
82  * The arguments to `jacobian`
83  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
84  *
85  * Differentiating
86  * Eqs. (\f$\ref{eq:forward_map_x}\f$--\f$\ref{eq:forward_map_z}\f$)
87  * above yields
88  *
89  * \f{align*}
90  * \frac{\partial x_0^0}{\partial \bar{x}} &=
91  * R_\mathrm{out}-R_\mathrm{in} + (2 R_\mathrm{in}-R_\mathrm{out})
92  * \frac{\bar{y}^2}{\bar{\rho}^3},\\
93  * \frac{\partial x_0^0}{\partial \bar{y}} &=
94  * -(2 R_\mathrm{in}-R_\mathrm{out})
95  * \frac{\bar{x}\bar{y}}{\bar{\rho}^3},\\
96  * \frac{\partial x_0^1}{\partial \bar{x}} &=
97  * -(2 R_\mathrm{in}-R_\mathrm{out})
98  * \frac{\bar{x}\bar{y}}{\bar{\rho}^3},\\
99  * \frac{\partial x_0^1}{\partial \bar{y}} &=
100  * R_\mathrm{out}-R_\mathrm{in} + (2 R_\mathrm{in}-R_\mathrm{out})
101  * \frac{\bar{x}^2}{\bar{\rho}^3},\\
102  * \f}
103  * and all other components are zero.
104  *
105  * ### inverse
106  *
107  * `inverse` takes \f$x_0^i\f$ and \f$\sigma\f$ as arguments, and
108  * returns \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed
109  * `std::optional<std::array<double, 3>>` if
110  * \f$x_0^i\f$ or \f$\sigma\f$ are outside the range of the map.
111  *
112  * Let
113  * \f{align}
114  * \rho = \sqrt{(x_0^0-C^0)^2+(x_0^1-C^1)^2}.
115  * \label{eq:rho}
116  * \f}
117  *
118  * Then
119  * \f{align}
120  * \bar{x} &= \frac{x_0^0-C^0}{\rho}
121  * \frac{\rho+R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
122  * \bar{y} &= \frac{x_0^1-C^1}{\rho}
123  * \frac{\rho+R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
124  * \bar{z} &= 2\sigma - 1.
125  * \f}
126  *
127  * Note that \f$\rho\f$ in Eq. (\f$\ref{eq:rho}\f$) can be written
128  * \f{align}
129  * \rho = R_\mathrm{in}+(R_\mathrm{out}-R_\mathrm{in})(\bar{\rho}-1),
130  * \label{eq:rho_from_rhobar}
131  * \f}
132  * where \f$\bar{\rho}\f$ is given by Eq. (\f$\ref{eq:rhobar}\f$).
133  *
134  * If \f$\bar{z}\f$ is outside the range \f$[-1,1]\f$ or
135  * if \f$\bar{x}^2+\bar{y}^2\f$ is less than 1 or greater than 4
136  * then we return a default-constructed
137  * `std::optional<std::array<double, 3>>`.
138  *
139  * ### lambda_tilde
140  *
141  * `lambda_tilde` takes as arguments a point \f$x^i\f$ and a projection point
142  * \f$P^i\f$, and computes \f$\tilde{\lambda}\f$, the solution to
143  *
144  * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
145  *
146  * Since \f$x_0^i\f$ must lie on the plane \f$x_0^3=C^3\f$,
147  *
148  * \f{align} \tilde{\lambda} &= \frac{C^3-P^3}{x^3-P^3}.\f}
149  *
150  * If \f$\tilde{\lambda}\f$ is less than unity (indicating that the
151  * supplied point is outside the range of the map), then a
152  * default-constructed `std::optional<double>` is returned.
153  *
154  * ### deriv_lambda_tilde
155  *
156  * `deriv_lambda_tilde` takes as arguments \f$x_0^i\f$, a projection point
157  * \f$P^i\f$, and \f$\tilde{\lambda}\f$, and
158  * returns \f$\partial \tilde{\lambda}/\partial x^i\f$. We have
159  *
160  * \f{align}
161  * \frac{\partial\tilde{\lambda}}{\partial x^3} =
162  * -\frac{C^3-P^3}{(x^3-P^3)^2} = -\frac{\tilde{\lambda}^2}{C^3-P^3},
163  * \f}
164  * and other components are zero.
165  *
166  * ### inv_jacobian
167  *
168  * `inv_jacobian` returns \f$\partial \bar{x}^i/\partial x_0^k\f$,
169  * where \f$\sigma\f$ is held fixed.
170  * The arguments to `inv_jacobian`
171  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
172  *
173  * The nonzero components are
174  * \f{align}
175  * \frac{\partial \bar{x}}{\partial x_0^0} &=
176  * \frac{1}{R_\mathrm{out}-R_\mathrm{in}} + \frac{\bar{y}^2}{\bar{\rho}^2\rho}
177  * \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
178  * \frac{\partial \bar{x}}{\partial x_0^1} &=
179  * - \frac{\bar{x}\bar{y}}{\bar{\rho}^2\rho}
180  * \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
181  * \frac{\partial \bar{y}}{\partial x_0^1} &=
182  * \frac{1}{R_\mathrm{out}-R_\mathrm{in}} + \frac{\bar{x}^2}{\bar{\rho}^2\rho}
183  * \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},\\
184  * \frac{\partial \bar{y}}{\partial x_0^0} &=
185  * - \frac{\bar{x}\bar{y}}{\bar{\rho}^2\rho}
186  * \frac{R_\mathrm{out}-2 R_\mathrm{in}}{R_\mathrm{out}-R_\mathrm{in}},
187  * \f}
188  * where \f$\rho\f$ is computed from Eq. (\f$\ref{eq:rho_from_rhobar}\f$).
189  *
190  * ### dxbar_dsigma
191  *
192  * `dxbar_dsigma` returns \f$\partial \bar{x}^i/\partial \sigma\f$,
193  * where \f$x_0^i\f$ is held fixed.
194  *
195  * From Eq. (\f$\ref{eq:deriv_sigma}\f$) we have
196  *
197  * \f{align}
198  * \frac{\partial \bar{x}^i}{\partial \sigma} &= (0,0,2).
199  * \f}
200  *
201  */
202 class FlatSide {
203  public:
204  FlatSide(const std::array<double, 3>& center, const double inner_radius,
205  const double outer_radius) noexcept;
206 
207  FlatSide() = default;
208  ~FlatSide() = default;
209  FlatSide(FlatSide&&) = default;
210  FlatSide(const FlatSide&) = default;
211  FlatSide& operator=(const FlatSide&) = default;
212  FlatSide& operator=(FlatSide&&) = default;
213 
214  template <typename T>
215  void forward_map(
216  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
217  target_coords,
218  const std::array<T, 3>& source_coords) const noexcept;
219 
221  const std::array<double, 3>& target_coords,
222  double sigma_in) const noexcept;
223 
224  template <typename T>
225  void jacobian(const gsl::not_null<
226  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
227  jacobian_out,
228  const std::array<T, 3>& source_coords) const noexcept;
229 
230  template <typename T>
231  void inv_jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3,
232  Frame::NoFrame>*>
233  inv_jacobian_out,
234  const std::array<T, 3>& source_coords) const noexcept;
235 
236  template <typename T>
237  void sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,
238  const std::array<T, 3>& source_coords) const noexcept;
239 
240  template <typename T>
241  void deriv_sigma(
242  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
243  deriv_sigma_out,
244  const std::array<T, 3>& source_coords) const noexcept;
245 
246  template <typename T>
247  void dxbar_dsigma(
248  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
249  dxbar_dsigma_out,
250  const std::array<T, 3>& source_coords) const noexcept;
251 
252  std::optional<double> lambda_tilde(
253  const std::array<double, 3>& parent_mapped_target_coords,
254  const std::array<double, 3>& projection_point,
255  bool source_is_between_focus_and_target) const noexcept;
256 
257  template <typename T>
258  void deriv_lambda_tilde(
259  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
260  deriv_lambda_tilde_out,
261  const std::array<T, 3>& target_coords, const T& lambda_tilde,
262  const std::array<double, 3>& projection_point) const noexcept;
263 
264  // clang-tidy: google runtime references
265  void pup(PUP::er& p) noexcept; // NOLINT
266 
267  static bool is_identity() noexcept { return false; }
268 
269  private:
270  friend bool operator==(const FlatSide& lhs, const FlatSide& rhs) noexcept;
271  std::array<double, 3> center_{};
272  double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
273  double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
274 };
275 bool operator!=(const FlatSide& lhs, const FlatSide& rhs) noexcept;
276 } // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps
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
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
domain::CoordinateMaps::FocallyLiftedInnerMaps::FlatSide
A FocallyLiftedInnerMap that maps a 3D unit right cylinder to a volume that connects a 2D annulus to ...
Definition: FocallyLiftedFlatSide.hpp:202
limits
Gsl.hpp
TypeAliases.hpp
optional
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13