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 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 1 : class FlatEndcap {
155 : public:
156 0 : FlatEndcap(const std::array<double, 3>& center, double radius);
157 :
158 0 : FlatEndcap() = default;
159 0 : ~FlatEndcap() = default;
160 0 : FlatEndcap(FlatEndcap&&) = default;
161 0 : FlatEndcap(const FlatEndcap&) = default;
162 0 : FlatEndcap& operator=(const FlatEndcap&) = default;
163 0 : FlatEndcap& operator=(FlatEndcap&&) = default;
164 :
165 : template <typename T>
166 0 : 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;
170 :
171 0 : std::optional<std::array<double, 3>> inverse(
172 : const std::array<double, 3>& target_coords, double sigma_in) const;
173 :
174 : template <typename T>
175 0 : void jacobian(const gsl::not_null<
176 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
177 : jacobian_out,
178 : const std::array<T, 3>& source_coords) const;
179 :
180 : template <typename T>
181 0 : void inv_jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3,
182 : Frame::NoFrame>*>
183 : inv_jacobian_out,
184 : const std::array<T, 3>& source_coords) const;
185 :
186 : template <typename T>
187 0 : void sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,
188 : const std::array<T, 3>& source_coords) const;
189 :
190 : template <typename T>
191 0 : void deriv_sigma(
192 : const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
193 : deriv_sigma_out,
194 : const std::array<T, 3>& source_coords) const;
195 :
196 : template <typename T>
197 0 : void dxbar_dsigma(
198 : const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
199 : dxbar_dsigma_out,
200 : const std::array<T, 3>& source_coords) const;
201 :
202 0 : std::optional<double> lambda_tilde(
203 : const std::array<double, 3>& parent_mapped_target_coords,
204 : const std::array<double, 3>& projection_point,
205 : bool source_is_between_focus_and_target) const;
206 :
207 : template <typename T>
208 0 : void deriv_lambda_tilde(
209 : const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
210 : deriv_lambda_tilde_out,
211 : const std::array<T, 3>& target_coords, const T& lambda_tilde,
212 : const std::array<double, 3>& projection_point) const;
213 :
214 : // NOLINTNEXTLINE(google-runtime-references)
215 0 : void pup(PUP::er& p);
216 :
217 0 : static bool is_identity() { return false; }
218 :
219 : private:
220 0 : friend bool operator==(const FlatEndcap& lhs, const FlatEndcap& rhs);
221 0 : std::array<double, 3> center_{};
222 0 : double radius_{std::numeric_limits<double>::signaling_NaN()};
223 : };
224 0 : bool operator!=(const FlatEndcap& lhs, const FlatEndcap& rhs);
225 : } // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps
|