Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines the class CylindricalSide.
6 :
7 : #pragma once
8 :
9 : #include <array>
10 : #include <cstddef>
11 : #include <limits>
12 : #include <optional>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Domain/CoordinateMaps/FocallyLiftedMap.hpp"
16 : #include "Domain/CoordinateMaps/FocallyLiftedSide.hpp"
17 :
18 : /// \cond
19 : namespace PUP {
20 : class er;
21 : } // namespace PUP
22 : /// \endcond
23 :
24 : namespace domain::CoordinateMaps {
25 :
26 : /*!
27 : * \ingroup CoordinateMapsGroup
28 : *
29 : * \brief Map from a 3D unit right cylindrical shell to a volume that connects
30 : * portions of two spherical surfaces.
31 : *
32 : * \image html CylindricalSide.svg "2D slice showing mapped (shaded) region."
33 : *
34 : * \details Consider two spheres with centers \f$C_1\f$ and \f$C_2\f$,
35 : * and radii \f$R_1\f$ and \f$R_2\f$. Let sphere 1 be intersected by two
36 : * planes normal to the \f$z\f$ axis and located at \f$z = z_\mathrm{L}\f$
37 : * and \f$z = z_\mathrm{U}\f$, with \f$z_\mathrm{L} < z_\mathrm{U}\f$.
38 : * Also let there be a projection point \f$P\f$.
39 : *
40 : * Note that Sphere 1 and Sphere 2 are not equivalent, because the
41 : * mapped portion of Sphere 1 is bounded by planes of constant
42 : * \f$z\f$ but the mapped portion of Sphere 2 is not (except for
43 : * special choices of \f$C_1\f$, \f$C_2\f$, and \f$P\f$).
44 : *
45 : * CylindricalSide maps a 3D unit right cylindrical shell (with
46 : * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that
47 : * \f$-1\leq\bar{z}\leq 1\f$ and \f$1 \leq \bar{x}^2+\bar{y}^2 \leq
48 : * 4\f$) to the shaded area in each panel of the figure above (with
49 : * coordinates \f$(x,y,z)\f$). The figure shows two different allowed
50 : * possibilities: \f$R_1 > R_2\f$ and \f$R_2 > R_1\f$. Note that the
51 : * two portions of the shaded region in each panel of the figure
52 : * represent different portions of the same block; each panel of the
53 : * figure is to be understood as rotated around the \f$z\f$ axis. The
54 : * inner boundary of the cylindrical shell \f$\bar{x}^2+\bar{y}^2=1\f$
55 : * is mapped to the portion of sphere 1 that has \f$z_\mathrm{L} \leq
56 : * z \leq z_\mathrm{U}\f$. Curves of constant \f$(\bar{z})\f$ along
57 : * the vector \f$(\bar{x},\bar{y})\f$ are mapped to portions of lines
58 : * that pass through \f$P\f$. Along each of these curves,
59 : * \f$\bar{x}^2+\bar{y}^2=1\f$ is mapped to a point on sphere 1 and
60 : * \f$\bar{x}^2+\bar{y}^2=4\f$ is mapped to a point on sphere 2.
61 : *
62 : * CylindricalSide is described briefly in the Appendix of
63 : * \cite Buchman:2012dw. CylindricalSide is used to construct the blocks
64 : * labeled 'CA cylinder', 'EA cylinder', 'CB cylinder', 'EE cylinder',
65 : * and 'EB cylinder' in Figure 20 of that paper. Note that 'CA
66 : * cylinder', 'CB cylinder', and 'EE cylinder' have Sphere 1 contained
67 : * in Sphere2, and 'EA cylinder' and 'EB cylinder' have Sphere 2
68 : * contained in Sphere 1.
69 : *
70 : * CylindricalSide is implemented using `FocallyLiftedMap`
71 : * and `FocallyLiftedInnerMaps::Side`; see those classes for
72 : * details.
73 : *
74 : * ### Restrictions on map parameters.
75 : *
76 : * We demand that:
77 : * - Either Sphere 1 is fully contained inside Sphere 2, or
78 : * Sphere 2 is fully contained inside Sphere 1.
79 : * - \f$P\f$ is contained inside the smaller sphere, and
80 : * between (or on) the planes defined by \f$z_\mathrm{L}\f$ and
81 : * \f$z_\mathrm{U}\f$.
82 : * - If sphere 1 is contained in sphere 2:
83 : * - \f$C_1^z - 0.95 R_1 \leq z_\mathrm{L}\f$
84 : * - \f$z_\mathrm{U} \leq C_1^z + 0.95 R_1\f$
85 : * - If sphere 2 is contained in sphere 1:
86 : * - \f$C_1^z - 0.95 R_1 \leq z_\mathrm{L} \leq C_1^z - 0.2 R_1\f$
87 : * - \f$C_1^z + 0.2 R_1 \leq z_\mathrm{U} \leq C_1^z + 0.95 R_1\f$
88 : *
89 : */
90 1 : class CylindricalSide {
91 : public:
92 0 : static constexpr size_t dim = 3;
93 0 : CylindricalSide(const std::array<double, 3>& center_one,
94 : const std::array<double, 3>& center_two,
95 : const std::array<double, 3>& proj_center, double radius_one,
96 : double radius_two, double z_lower, double z_upper);
97 :
98 0 : CylindricalSide() = default;
99 0 : ~CylindricalSide() = default;
100 0 : CylindricalSide(CylindricalSide&&) = default;
101 0 : CylindricalSide(const CylindricalSide&) = default;
102 0 : CylindricalSide& operator=(const CylindricalSide&) = default;
103 0 : CylindricalSide& operator=(CylindricalSide&&) = default;
104 :
105 : template <typename T>
106 0 : std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
107 :
108 0 : std::optional<std::array<double, 3>> inverse(
109 : const std::array<double, 3>& target_coords) const;
110 :
111 : template <typename T>
112 0 : tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
113 : const std::array<T, 3>& source_coords) const;
114 :
115 : template <typename T>
116 0 : tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
117 : const std::array<T, 3>& source_coords) const;
118 :
119 : // NOLINTNEXTLINE(google-runtime-references)
120 0 : void pup(PUP::er& p);
121 :
122 0 : static bool is_identity() { return false; }
123 :
124 0 : static constexpr bool supports_hessian{false};
125 :
126 : private:
127 0 : friend bool operator==(const CylindricalSide& lhs,
128 : const CylindricalSide& rhs);
129 0 : FocallyLiftedMap<FocallyLiftedInnerMaps::Side> impl_;
130 : };
131 0 : bool operator!=(const CylindricalSide& lhs, const CylindricalSide& rhs);
132 :
133 : } // namespace domain::CoordinateMaps
|
