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 "Domain/CoordinateMaps/Distribution.hpp"
13 : #include "Domain/Structure/OrientationMap.hpp"
14 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp"
15 :
16 : /// \cond
17 : namespace PUP {
18 : class er;
19 : } // namespace PUP
20 : /// \endcond
21 :
22 : namespace domain::CoordinateMaps {
23 :
24 : namespace detail {
25 : // This mapping can be deleted once the 2D and 3D wedges are oriented the same
26 : // (see issue https://github.com/sxs-collaboration/spectre/issues/2988)
27 : template <size_t Dim>
28 : struct WedgeCoordOrientation;
29 : template <>
30 : struct WedgeCoordOrientation<2> {
31 : static constexpr size_t radial_coord = 0;
32 : static constexpr size_t polar_coord = 1;
33 : static constexpr size_t azimuth_coord = 2; // unused
34 : };
35 : template <>
36 : struct WedgeCoordOrientation<3> {
37 : static constexpr size_t radial_coord = 2;
38 : static constexpr size_t polar_coord = 0;
39 : static constexpr size_t azimuth_coord = 1;
40 : };
41 : } // namespace detail
42 :
43 : /*!
44 : * \ingroup CoordinateMapsGroup
45 : *
46 : * \brief Map from a square or cube to a wedge.
47 : * \image html Shell.png "A shell can be constructed out of six wedges."
48 : *
49 : * \details The mapping that goes from a reference cube (in 3D) or square (in
50 : * 2D) to a wedge centered on a coordinate axis covering a volume between an
51 : * inner surface and outer surface. Each surface can be given a curvature
52 : * between flat (a sphericity of 0) or spherical (a sphericity of 1).
53 : *
54 : * In 2D, the first logical coordinate corresponds to the radial coordinate,
55 : * and the second logical coordinates correspond to the angular coordinate. In
56 : * 3D, the first two logical coordinates correspond to the two angular
57 : * coordinates, and the third to the radial coordinate. This difference
58 : * originates from separate implementations for the 2D and 3D map that were
59 : * merged. The 3D implementation can be changed to use the first logical
60 : * coordinate as radial direction, but this requires propagating the change
61 : * through the rest of the domain code (see issue
62 : * https://github.com/sxs-collaboration/spectre/issues/2988).
63 : *
64 : * The following documentation is for the 3D map. The 2D map is obtained by
65 : * setting either of the two angular coordinates to zero (and using \f$\xi\f$
66 : * as radial coordinate).
67 : *
68 : * The Wedge map is constructed by linearly interpolating between a bulged
69 : * face of radius `radius_of_inner_surface` to a bulged face of
70 : * radius `radius_of_outer_surface`, where the radius of each bulged face
71 : * is defined to be the radius of the sphere circumscribing the bulge.
72 : *
73 : * We make a choice here as to whether we wish to use the logical coordinates
74 : * parameterizing these surface as they are, in which case we have the
75 : * equidistant choice of coordinates, or whether to apply a tangent map to them
76 : * which leads us to the equiangular choice of coordinates. In terms of the
77 : * logical coordinates, the equiangular coordinates are:
78 : *
79 : * \f[\textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4)\f]
80 : *
81 : * \f[\textrm{equiangular eta} : \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)\f]
82 : *
83 : * With derivatives:
84 : *
85 : * \f[\Xi'(\xi) = \frac{\pi}{4}(1+\Xi^2)\f]
86 : *
87 : * \f[\mathrm{H}'(\eta) = \frac{\pi}{4}(1+\mathrm{H}^2)\f]
88 : *
89 : * The equidistant coordinates are:
90 : *
91 : * \f[ \textrm{equidistant xi} : \Xi = \xi\f]
92 : *
93 : * \f[ \textrm{equidistant eta} : \mathrm{H} = \eta\f]
94 : *
95 : * with derivatives:
96 : *
97 : * <center>\f$\Xi'(\xi) = 1\f$, and \f$\mathrm{H}'(\eta) = 1\f$</center>
98 : *
99 : * We also define the variable \f$\rho\f$, given by:
100 : *
101 : * \f[\textrm{rho} : \rho = \sqrt{1+\Xi^2+\mathrm{H}^2}\f]
102 : *
103 : * ### The Spherical Face Map
104 : * The surface map for the spherical face of radius \f$R\f$ lying in the
105 : * \f$+z\f$
106 : * direction in either choice of coordinates is then given by:
107 : *
108 : * \f[\vec{\sigma}_{spherical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\f]
109 : * Where
110 : * \f[
111 : * \vec{x}(\xi,\eta) =
112 : * \begin{bmatrix}
113 : * x(\xi,\eta)\\
114 : * y(\xi,\eta)\\
115 : * z(\xi,\eta)\\
116 : * \end{bmatrix} = \frac{R}{\rho}
117 : * \begin{bmatrix}
118 : * \Xi\\
119 : * \mathrm{H}\\
120 : * 1\\
121 : * \end{bmatrix}\f]
122 : *
123 : * ### The Bulged Face Map
124 : * The bulged surface is itself constructed by linearly interpolating between
125 : * a cubical face and a spherical face. The surface map for the cubical face
126 : * of side length \f$2L\f$ lying in the \f$+z\f$ direction is given by:
127 : *
128 : * \f[\vec{\sigma}_{cubical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\f]
129 : * Where
130 : * \f[
131 : * \vec{x}(\xi,\eta) =
132 : * \begin{bmatrix}
133 : * x(\xi,\eta)\\
134 : * y(\xi,\eta)\\
135 : * L\\
136 : * \end{bmatrix} = L
137 : * \begin{bmatrix}
138 : * \Xi\\
139 : * \mathrm{H}\\
140 : * 1\\
141 : * \end{bmatrix}\f]
142 : *
143 : * To construct the bulged map we interpolate between this cubical face map
144 : * and a spherical face map of radius \f$R\f$, with the
145 : * interpolation parameter being \f$s\f$. The surface map for the bulged face
146 : * lying in the \f$+z\f$ direction is then given by:
147 : *
148 : * \f[\vec{\sigma}_{bulged}(\xi,\eta) = {(1-s)L + \frac{sR}{\rho}}
149 : * \begin{bmatrix}
150 : * \Xi\\
151 : * \mathrm{H}\\
152 : * 1\\
153 : * \end{bmatrix}\f]
154 : *
155 : * We constrain L by demanding that the spherical face circumscribe the cube.
156 : * With this condition, we have \f$L = R/\sqrt3\f$.
157 : * \note This differs from the choice in SpEC where it is demanded that the
158 : * surfaces touch at the center, which leads to \f$L = R\f$.
159 : *
160 : * ### The Full Volume Map
161 : * The final map for the wedge which lies along the \f$+z\f$ is obtained by
162 : * interpolating between the two surfaces with the
163 : * interpolation parameter being the logical coordinate \f$\zeta\f$. This
164 : * results in:
165 : *
166 : * \f[\vec{x}(\xi,\eta,\zeta) =
167 : * \frac{1}{2}\left\{(1-\zeta)\Big[(1-s_{inner})\frac{R_{inner}}{\sqrt 3}
168 : * + s_{inner}\frac{R_{inner}}{\rho}\Big] +
169 : * (1+\zeta)\Big[(1-s_{outer})\frac{R_{outer}}{\sqrt 3} +s_{outer}
170 : * \frac{R_{outer}}{\rho}\Big] \right\}\begin{bmatrix}
171 : * \Xi\\
172 : * \mathrm{H}\\
173 : * 1\\
174 : * \end{bmatrix}\f]
175 : *
176 : * We will define the variables \f$F(\zeta)\f$ and \f$S(\zeta)\f$, the frustum
177 : * and sphere factors: \f[F(\zeta) = F_0 + F_1\zeta\f] \f[S(\zeta) = S_0 +
178 : * S_1\zeta\f]
179 : * Where \f{align*}F_0 &= \frac{1}{2} \big\{ (1-s_{outer})R_{outer} +
180 : * (1-s_{inner})R_{inner}\big\}\\
181 : * F_1 &= \partial_{\zeta} F = \frac{1}{2} \big\{ (1-s_{outer})R_{outer} -
182 : * (1-s_{inner})R_{inner}\big\}\\
183 : * S_0 &= \frac{1}{2} \big\{ s_{outer}R_{outer} + s_{inner}R_{inner}\big\}\\
184 : * S_1 &= \partial_{\zeta} S = \frac{1}{2} \big\{ s_{outer}R_{outer} -
185 : * s_{inner}R_{inner}\big\}\f}
186 : *
187 : * The map can then be rewritten as:
188 : * \f[\vec{x}(\xi,\eta,\zeta) = \left\{\frac{F(\zeta)}{\sqrt 3} +
189 : * \frac{S(\zeta)}{\rho}\right\}\begin{bmatrix}
190 : * \Xi\\
191 : * \mathrm{H}\\
192 : * 1\\
193 : * \end{bmatrix}\f]
194 : *
195 : * We provide some common derivatives:
196 : * \f[\partial_{\xi}z = \frac{-S(\zeta)\Xi\Xi'}{\rho^3}\f]
197 : * \f[\partial_{\eta}z = \frac{-S(\zeta)\mathrm{H}\mathrm{H}'}{\rho^3}\f]
198 : * \f[\partial_{\zeta}z = \frac{F'}{\sqrt 3} + \frac{S'}{\rho}\f]
199 : * The Jacobian then is: \f[J =
200 : * \begin{bmatrix}
201 : * \Xi'z + \Xi\partial_{\xi}z & \Xi\partial_{\eta}z & \Xi\partial_{\zeta}z \\
202 : * \mathrm{H}\partial_{\xi}z & \mathrm{H}'z +
203 : * \mathrm{H}\partial_{\eta}z & \mathrm{H}\partial_{\zeta}z\\
204 : * \partial_{\xi}z&\partial_{\eta}z &\partial_{\zeta}z \\
205 : * \end{bmatrix}
206 : * \f]
207 : *
208 : * A common factor that shows up in the inverse jacobian is:
209 : * \f[ T:= \frac{S(\zeta)}{(\partial_{\zeta}z)\rho^3}\f]
210 : *
211 : * The inverse Jacobian then is: \f[J^{-1} =
212 : * \frac{1}{z}\begin{bmatrix}
213 : * \Xi'^{-1} & 0 & -\Xi\Xi'^{-1}\\
214 : * 0 & \mathrm{H}'^{-1} & -\mathrm{H}\mathrm{H}'^{-1}\\
215 : * T\Xi &
216 : * T\mathrm{H} &
217 : * T + F(\partial_{\zeta}z)^{-1}/\sqrt 3\\
218 : * \end{bmatrix}
219 : * \f]
220 : *
221 : * ### Changing the radial distribution of the gridpoints
222 : * By default, Wedge linearly distributes its gridpoints in the radial
223 : * direction. An exponential distribution of gridpoints can be obtained by
224 : * linearly interpolating in the logarithm of the radius, in order to obtain
225 : * a relatively higher resolution at smaller radii. Since this is a radial
226 : * rescaling of Wedge, this option is only supported for fully spherical
227 : * wedges with `sphericity_inner` = `sphericity_outer` = 1.
228 : *
229 : * The linear interpolation done is:
230 : * \f[
231 : * \log r = \frac{1-\zeta}{2}\log R_{inner} +
232 : * \frac{1+\zeta}{2}\log R_{outer}
233 : * \f]
234 : *
235 : * The map then is:
236 : * \f[\vec{x}(\xi,\eta,\zeta) =
237 : * \frac{\sqrt{R_{inner}^{1-\zeta}R_{outer}^{1+\zeta}}}{\rho}\begin{bmatrix}
238 : * \Xi\\
239 : * \mathrm{H}\\
240 : * 1\\
241 : * \end{bmatrix}\f]
242 : *
243 : * The jacobian simplifies similarly.
244 : *
245 : * Alternatively, an inverse radial distribution can be chosen where the linear
246 : * interpolation is:
247 : *
248 : * \f[
249 : * \frac{1}{r} = \frac{R_\mathrm{inner} + R_\mathrm{outer}}{2 R_\mathrm{inner}
250 : * R_\mathrm{outer}} + \frac{R_\mathrm{inner} - R_\mathrm{outer}}{2
251 : * R_\mathrm{inner} R_\mathrm{outer}} \zeta
252 : * \f]
253 : */
254 : template <size_t Dim>
255 1 : class Wedge {
256 : public:
257 0 : static constexpr size_t dim = Dim;
258 0 : enum class WedgeHalves {
259 : /// Use the entire wedge
260 : Both,
261 : /// Use only the upper logical half
262 : UpperOnly,
263 : /// Use only the lower logical half
264 : LowerOnly
265 : };
266 :
267 : /*!
268 : * Constructs a 3D wedge.
269 : * \param radius_inner Distance from the origin to one of the
270 : * corners which lie on the inner surface.
271 : * \param radius_outer Distance from the origin to one of the
272 : * corners which lie on the outer surface.
273 : * \param orientation_of_wedge The orientation of the desired wedge relative
274 : * to the orientation of the default wedge which is a wedge that has its
275 : * curved surfaces pierced by the upper-z axis. The logical xi and eta
276 : * coordinates point in the cartesian x and y directions, respectively.
277 : * \param sphericity_inner Value between 0 and 1 which determines
278 : * whether the inner surface is flat (value of 0), spherical (value of 1) or
279 : * somewhere in between
280 : * \param sphericity_outer Value between 0 and 1 which determines
281 : * whether the outer surface is flat (value of 0), spherical (value of 1) or
282 : * somewhere in between
283 : * \param with_equiangular_map Determines whether to apply a tangent function
284 : * mapping to the logical coordinates (for `true`) or not (for `false`).
285 : * \param halves_to_use Determines whether to construct a full wedge or only
286 : * half a wedge. If constructing only half a wedge, the resulting shape has a
287 : * face normal to the x direction (assuming default OrientationMap). If
288 : * constructing half a wedge, an intermediate affine map is applied to the
289 : * logical xi coordinate such that the interval [-1,1] is mapped to the
290 : * corresponding logical half of the wedge. For example, if `UpperOnly` is
291 : * specified, [-1,1] is mapped to [0,1], and if `LowerOnly` is specified,
292 : * [-1,1] is mapped to [-1,0]. The case of `Both` means a full wedge, with no
293 : * intermediate map applied. In all cases, the logical points returned by the
294 : * inverse map will lie in the range [-1,1] in each dimension. Half wedges are
295 : * currently only useful in constructing domains for binary systems.
296 : * \param radial_distribution Determines how to distribute grid points along
297 : * the radial direction. For wedges that are not exactly spherical, only
298 : * `Distribution::Linear` is currently supported.
299 : * \param opening_angles Determines the angular size of the wedge. The
300 : * default value is pi/2, which corresponds to a wedge size of pi/2. For this
301 : * setting, four Wedges can be put together to cover 2pi in angle along a
302 : * great circle. This option is meant to be used with the equiangular map
303 : * option turned on.
304 : * \param with_adapted_equiangular_map Determines whether to adapt the
305 : * point distribution in the wedge to match its physical angular size. When
306 : * `true`, angular distances are proportional to logical distances. Note
307 : * that it is not possible to use adapted maps in every Wedge of a Sphere
308 : * unless each Wedge has the same size along both angular directions.
309 : */
310 1 : Wedge(double radius_inner, double radius_outer, double sphericity_inner,
311 : double sphericity_outer, OrientationMap<Dim> orientation_of_wedge,
312 : bool with_equiangular_map,
313 : WedgeHalves halves_to_use = WedgeHalves::Both,
314 : Distribution radial_distribution = Distribution::Linear,
315 : const std::array<double, Dim - 1>& opening_angles =
316 : make_array<Dim - 1>(M_PI_2),
317 : bool with_adapted_equiangular_map = true);
318 :
319 0 : Wedge() = default;
320 0 : ~Wedge() = default;
321 0 : Wedge(Wedge&&) = default;
322 0 : Wedge(const Wedge&) = default;
323 0 : Wedge& operator=(const Wedge&) = default;
324 0 : Wedge& operator=(Wedge&&) = default;
325 :
326 : template <typename T>
327 0 : std::array<tt::remove_cvref_wrap_t<T>, Dim> operator()(
328 : const std::array<T, Dim>& source_coords) const;
329 :
330 : /// For a \f$+z\f$-oriented `Wedge`, returns invalid if \f$z<=0\f$
331 : /// or if \f$(x,y,z)\f$ is on or outside the cone defined
332 : /// by \f$(x^2/z^2 + y^2/z^2+1)^{1/2} = -S/F\f$, where
333 : /// \f$S = \frac{1}{2}(s_1 r_1 - s_0 r_0)\f$ and
334 : /// \f$F = \frac{1}{2\sqrt{3}}((1-s_1) r_1 - (1-s_0) r_0)\f$.
335 : /// Here \f$s_0,s_1\f$ and \f$r_0,r_1\f$ are the specified sphericities
336 : /// and radii of the inner and outer \f$z\f$ surfaces. The map is singular on
337 : /// the cone and on the xy plane.
338 : /// The inverse function is only callable with doubles because the inverse
339 : /// might fail if called for a point out of range, and it is unclear
340 : /// what should happen if the inverse were to succeed for some points in a
341 : /// DataVector but fail for other points.
342 1 : std::optional<std::array<double, Dim>> inverse(
343 : const std::array<double, Dim>& target_coords) const;
344 :
345 : template <typename T>
346 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> jacobian(
347 : const std::array<T, Dim>& source_coords) const;
348 :
349 : template <typename T>
350 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, Dim, Frame::NoFrame> inv_jacobian(
351 : const std::array<T, Dim>& source_coords) const;
352 :
353 : // NOLINTNEXTLINE(google-runtime-references)
354 0 : void pup(PUP::er& p);
355 :
356 0 : static constexpr bool is_identity() { return false; }
357 :
358 : private:
359 : // maps between 2D and 3D choices for coordinate axis orientations
360 0 : static constexpr size_t radial_coord =
361 : detail::WedgeCoordOrientation<Dim>::radial_coord;
362 0 : static constexpr size_t polar_coord =
363 : detail::WedgeCoordOrientation<Dim>::polar_coord;
364 0 : static constexpr size_t azimuth_coord =
365 : detail::WedgeCoordOrientation<Dim>::azimuth_coord;
366 :
367 : // factors out calculation of z needed for mapping and jacobian
368 : template <typename T>
369 0 : tt::remove_cvref_wrap_t<T> default_physical_z(const T& zeta,
370 : const T& one_over_rho) const;
371 :
372 : template <size_t LocalDim>
373 : // NOLINTNEXTLINE(readability-redundant-declaration)
374 0 : friend bool operator==(const Wedge<LocalDim>& lhs,
375 : const Wedge<LocalDim>& rhs);
376 :
377 0 : double radius_inner_{std::numeric_limits<double>::signaling_NaN()};
378 0 : double radius_outer_{std::numeric_limits<double>::signaling_NaN()};
379 0 : double sphericity_inner_{std::numeric_limits<double>::signaling_NaN()};
380 0 : double sphericity_outer_{std::numeric_limits<double>::signaling_NaN()};
381 0 : OrientationMap<Dim> orientation_of_wedge_{};
382 0 : bool with_equiangular_map_ = false;
383 0 : WedgeHalves halves_to_use_ = WedgeHalves::Both;
384 0 : Distribution radial_distribution_ = Distribution::Linear;
385 0 : double scaled_frustum_zero_{std::numeric_limits<double>::signaling_NaN()};
386 0 : double sphere_zero_{std::numeric_limits<double>::signaling_NaN()};
387 0 : double scaled_frustum_rate_{std::numeric_limits<double>::signaling_NaN()};
388 0 : double sphere_rate_{std::numeric_limits<double>::signaling_NaN()};
389 0 : std::array<double, Dim - 1> opening_angles_{
390 : make_array<Dim - 1>(std::numeric_limits<double>::signaling_NaN())};
391 0 : std::array<double, Dim - 1> opening_angles_distribution_{
392 : make_array<Dim - 1>(std::numeric_limits<double>::signaling_NaN())};
393 : };
394 :
395 : template <size_t Dim>
396 0 : bool operator!=(const Wedge<Dim>& lhs, const Wedge<Dim>& rhs);
397 : } // namespace domain::CoordinateMaps
|