FocallyLiftedFlatEndcap.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 portion of a plane and 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$\bar{x}^2+\bar{y}^2 \leq
31  * 1\f$. The range of the map has coordinates \f$(x,y,z)\f$.
32  *
33  * Consider a 2D circle in 3D space that is normal to the \f$z\f$ axis
34  * and has (3D) center \f$C^i\f$ and radius \f$R\f$. `FlatEndcap`
35  * provides the following functions:
36  *
37  * ### forward_map()
38  * `forward_map()` maps \f$(\bar{x},\bar{y},\bar{z}=-1)\f$ to the interior
39  * of the circle. The arguments to `forward_map()`
40  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
41  * `forward_map()` returns \f$x_0^i\f$,
42  * the 3D coordinates on the circle, which are given by
43  *
44  * \f{align}
45  * x_0^0 &= R \bar{x} + C^0,\\
46  * x_0^1 &= R \bar{y} + C^1,\\
47  * x_0^2 &= C^2.
48  * \f}
49  *
50  * ### sigma
51  *
52  * \f$\sigma\f$ is a function that is zero on the plane
53  * \f$x^i=x_0^i\f$ and unity at \f$\bar{z}=+1\f$ (corresponding to the
54  * upper surface of the FocallyLiftedMap). We define
55  *
56  * \f{align}
57  * \sigma &= \frac{\bar{z}+1}{2}.
58  * \f}
59  *
60  * ### deriv_sigma
61  *
62  * `deriv_sigma` returns
63  *
64  * \f{align}
65  * \frac{\partial \sigma}{\partial \bar{x}^j} &= (0,0,1/2).
66  * \f}
67  *
68  * ### jacobian
69  *
70  * `jacobian` returns \f$\partial x_0^k/\partial \bar{x}^j\f$.
71  * The arguments to `jacobian`
72  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
73  *
74  * Differentiating Eqs.(1--3) above yields
75  *
76  * \f{align*}
77  * \frac{\partial x_0^0}{\partial \bar{x}} &= R,\\
78  * \frac{\partial x_0^1}{\partial \bar{y}} &= R,
79  * \f}
80  * and all other components are zero.
81  *
82  * ### inverse
83  *
84  * `inverse` takes \f$x_0^i\f$ and \f$\sigma\f$ as arguments, and
85  * returns \f$(\bar{x},\bar{y},\bar{z})\f$, or a default-constructed
86  * `std::optional<std::array<double, 3>>` if
87  * \f$x_0^i\f$ or \f$\sigma\f$ are outside the range of the map.
88  * The formula for the inverse is straightforward:
89  *
90  * \f{align}
91  * \bar{x} &= \frac{x_0^0-C^0}{R},\\
92  * \bar{y} &= \frac{x_0^1-C^1}{R},\\
93  * \bar{z} &= 2\sigma - 1.
94  * \f}
95  *
96  * If \f$\bar{z}\f$ is outside the range \f$[-1,1]\f$ or
97  * if \f$\bar{x}^2+\bar{y}^2 > 1\f$ then we return
98  * a default-constructed `std::optional<std::array<double, 3>>`
99  *
100  * ### lambda_tilde
101  *
102  * `lambda_tilde` takes as arguments a point \f$x^i\f$ and a projection point
103  * \f$P^i\f$, and computes \f$\tilde{\lambda}\f$, the solution to
104  *
105  * \f{align} x_0^i = P^i + (x^i - P^i) \tilde{\lambda}.\f}
106  *
107  * Since \f$x_0^i\f$ must lie on the plane \f$x_0^3=C^3\f$,
108  *
109  * \f{align} \tilde{\lambda} &= \frac{C^3-P^3}{x^3-P^3}.\f}
110  *
111  * For `FocallyLiftedInnerMaps::FlatEndcap`, \f$x^i\f$ is always between
112  * \f$P^i\f$ and \f$x_0^i\f$, so \f$\tilde{\lambda}\ge 1\f$.
113  * Therefore a default-constructed `std::optional<double>` is returned
114  * if \f$\tilde{\lambda}\f$ is less than unity (meaning that \f$x^i\f$
115  * is outside the range of the map).
116  *
117  * ### deriv_lambda_tilde
118  *
119  * `deriv_lambda_tilde` takes as arguments \f$x_0^i\f$, a projection point
120  * \f$P^i\f$, and \f$\tilde{\lambda}\f$, and
121  * returns \f$\partial \tilde{\lambda}/\partial x^i\f$. We have
122  *
123  * \f{align}
124  * \frac{\partial\tilde{\lambda}}{\partial x^3} =
125  * -\frac{C^3-P^3}{(x^3-P^3)^2} = -\frac{\tilde{\lambda}^2}{C^3-P^3},
126  * \f}
127  * and other components are zero.
128  *
129  * ### inv_jacobian
130  *
131  * `inv_jacobian` returns \f$\partial \bar{x}^i/\partial x_0^k\f$,
132  * where \f$\sigma\f$ is held fixed.
133  * The arguments to `inv_jacobian`
134  * are \f$(\bar{x},\bar{y},\bar{z})\f$, but \f$\bar{z}\f$ is ignored.
135  *
136  * The nonzero components are
137  * \f{align}
138  * \frac{\partial \bar{x}}{\partial x_0^0} &= \frac{1}{R},\\
139  * \frac{\partial \bar{y}}{\partial x_0^1} &= \frac{1}{R}.
140  * \f}
141  *
142  * ### dxbar_dsigma
143  *
144  * `dxbar_dsigma` returns \f$\partial \bar{x}^i/\partial \sigma\f$,
145  * where \f$x_0^i\f$ is held fixed.
146  *
147  * From Eq. (6) we have
148  *
149  * \f{align}
150  * \frac{\partial \bar{x}^i}{\partial \sigma} &= (0,0,2).
151  * \f}
152  *
153  */
154 class FlatEndcap {
155  public:
156  FlatEndcap(const std::array<double, 3>& center, double radius) noexcept;
157 
158  FlatEndcap() = default;
159  ~FlatEndcap() = default;
160  FlatEndcap(FlatEndcap&&) = default;
161  FlatEndcap(const FlatEndcap&) = default;
162  FlatEndcap& operator=(const FlatEndcap&) = default;
163  FlatEndcap& operator=(FlatEndcap&&) = default;
164 
165  template <typename T>
166  void forward_map(
167  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
168  target_coords,
169  const std::array<T, 3>& source_coords) const noexcept;
170 
172  const std::array<double, 3>& target_coords,
173  double sigma_in) const noexcept;
174 
175  template <typename T>
176  void jacobian(const gsl::not_null<
177  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
178  jacobian_out,
179  const std::array<T, 3>& source_coords) const noexcept;
180 
181  template <typename T>
182  void inv_jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3,
183  Frame::NoFrame>*>
184  inv_jacobian_out,
185  const std::array<T, 3>& source_coords) const noexcept;
186 
187  template <typename T>
188  void sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,
189  const std::array<T, 3>& source_coords) const noexcept;
190 
191  template <typename T>
192  void deriv_sigma(
193  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
194  deriv_sigma_out,
195  const std::array<T, 3>& source_coords) const noexcept;
196 
197  template <typename T>
198  void dxbar_dsigma(
199  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
200  dxbar_dsigma_out,
201  const std::array<T, 3>& source_coords) const noexcept;
202 
203  std::optional<double> lambda_tilde(
204  const std::array<double, 3>& parent_mapped_target_coords,
205  const std::array<double, 3>& projection_point,
206  bool source_is_between_focus_and_target) const noexcept;
207 
208  template <typename T>
209  void deriv_lambda_tilde(
210  const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
211  deriv_lambda_tilde_out,
212  const std::array<T, 3>& target_coords, const T& lambda_tilde,
213  const std::array<double, 3>& projection_point) const noexcept;
214 
215  // clang-tidy: google runtime references
216  void pup(PUP::er& p) noexcept; // NOLINT
217 
218  static bool is_identity() noexcept { return false; }
219 
220  private:
221  friend bool operator==(const FlatEndcap& lhs, const FlatEndcap& rhs) noexcept;
222  std::array<double, 3> center_{};
224 };
225 bool operator!=(const FlatEndcap& lhs, const FlatEndcap& rhs) noexcept;
226 } // 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
domain::CoordinateMaps::FocallyLiftedInnerMaps::FlatEndcap
A FocallyLiftedInnerMap that maps a 3D unit right cylinder to a volume that connects a portion of a p...
Definition: FocallyLiftedFlatEndcap.hpp:154
cstddef
array
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
limits
Gsl.hpp
TypeAliases.hpp
optional
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13