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