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