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 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
14 :
15 : /// \cond
16 : namespace PUP {
17 : class er;
18 : } // namespace PUP
19 : /// \endcond
20 :
21 : /// Contains FocallyLiftedInnerMaps
22 : namespace domain::CoordinateMaps::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 1 : class FlatSide {
203 : public:
204 0 : FlatSide(const std::array<double, 3>& center, const double inner_radius,
205 : const double outer_radius);
206 :
207 0 : FlatSide() = default;
208 0 : ~FlatSide() = default;
209 0 : FlatSide(FlatSide&&) = default;
210 0 : FlatSide(const FlatSide&) = default;
211 0 : FlatSide& operator=(const FlatSide&) = default;
212 0 : FlatSide& operator=(FlatSide&&) = default;
213 :
214 : template <typename T>
215 0 : 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;
219 :
220 0 : std::optional<std::array<double, 3>> inverse(
221 : const std::array<double, 3>& target_coords, double sigma_in) const;
222 :
223 : template <typename T>
224 0 : void jacobian(const gsl::not_null<
225 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
226 : jacobian_out,
227 : const std::array<T, 3>& source_coords) const;
228 :
229 : template <typename T>
230 0 : void inv_jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3,
231 : Frame::NoFrame>*>
232 : inv_jacobian_out,
233 : const std::array<T, 3>& source_coords) const;
234 :
235 : template <typename T>
236 0 : void sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,
237 : const std::array<T, 3>& source_coords) const;
238 :
239 : template <typename T>
240 0 : void deriv_sigma(
241 : const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
242 : deriv_sigma_out,
243 : const std::array<T, 3>& source_coords) const;
244 :
245 : template <typename T>
246 0 : void dxbar_dsigma(
247 : const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
248 : dxbar_dsigma_out,
249 : const std::array<T, 3>& source_coords) const;
250 :
251 0 : std::optional<double> lambda_tilde(
252 : const std::array<double, 3>& parent_mapped_target_coords,
253 : const std::array<double, 3>& projection_point,
254 : bool source_is_between_focus_and_target) const;
255 :
256 : template <typename T>
257 0 : void deriv_lambda_tilde(
258 : const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
259 : deriv_lambda_tilde_out,
260 : const std::array<T, 3>& target_coords, const T& lambda_tilde,
261 : const std::array<double, 3>& projection_point) const;
262 :
263 : // NOLINTNEXTLINE(google-runtime-references)
264 0 : void pup(PUP::er& p);
265 :
266 0 : static bool is_identity() { return false; }
267 :
268 : private:
269 0 : friend bool operator==(const FlatSide& lhs, const FlatSide& rhs);
270 0 : std::array<double, 3> center_{};
271 0 : double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
272 0 : double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
273 : };
274 0 : bool operator!=(const FlatSide& lhs, const FlatSide& rhs);
275 : } // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps
|