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/FocallyLiftedFlatEndcap.hpp"
13 : #include "Domain/CoordinateMaps/FocallyLiftedMap.hpp"
14 :
15 : /// \cond
16 : namespace PUP {
17 : class er;
18 : } // namespace PUP
19 : /// \endcond
20 :
21 : namespace domain::CoordinateMaps {
22 :
23 : /*!
24 : * \ingroup CoordinateMapsGroup
25 : *
26 : * \brief Map from a 3D unit right cylinder to a volume that connects
27 : * a flat circular disk (lying inside a sphere) to the far wall of
28 : * that sphere.
29 : *
30 : * \details This is the "interior" counterpart to `CylindricalFlatEndcap`.
31 : * The two maps are identical in structure — both use `FocallyLiftedMap`
32 : * with `FocallyLiftedInnerMaps::FlatEndcap` — but differ in which
33 : * intersection of the projecting ray with the sphere is chosen:
34 : *
35 : * - `CylindricalFlatEndcap`: the flat disk lies *outside* the sphere
36 : * (the sphere is between \f$P\f$ and the flat disk), so
37 : * `source_is_between_focus_and_target = false`.
38 : * - `CylindricalFlatEndcapInterior`: the flat disk lies *inside* the
39 : * sphere (the flat disk is between \f$P\f$ and the sphere's far wall),
40 : * so `source_is_between_focus_and_target = true`.
41 : *
42 : * Consider a 2D circle in 3D space normal to the \f$z\f$ axis with
43 : * (3D) center \f$C_1\f$, a sphere with center \f$C_2\f$ and radius
44 : * \f$R_2\f$, and a projection point \f$P\f$.
45 : *
46 : * The parameter \f$z_\mathrm{extent}\f$ specifies the \f$z\f$-coordinate
47 : * (in the map's frame) of the rim circle where the spherical face of the
48 : * block meets the adjacent hollow-cylinder block. This single number
49 : * determines the radius \f$R_1\f$ of the flat disk:
50 : *
51 : * \f{align}{
52 : * t &= \frac{z_\mathrm{extent} - P_z}{C_1^z - P_z}, \\
53 : * r_\mathrm{rim} &= \sqrt{R_2^2 - (z_\mathrm{extent} - C_2^z)^2}, \\
54 : * R_1 &= \frac{r_\mathrm{rim}}{t}.
55 : * \f}
56 : *
57 : * CylindricalFlatEndcapInterior maps a 3D unit right cylinder
58 : * \f$(\bar{x},\bar{y},\bar{z})\f$ with \f$-1\leq\bar{z}\leq 1\f$ and
59 : * \f$\bar{x}^2+\bar{y}^2\leq 1\f$ so that:
60 : * - \f$\bar{z}=+1\f$ maps to the interior of the disk of radius
61 : * \f$R_1\f$ centred at \f$C_1\f$.
62 : * - \f$\bar{z}=-1\f$ maps to the portion of the sphere on the *far* side
63 : * of the flat disk from \f$P\f$.
64 : * - Curves of constant \f$(\bar{x},\bar{y})\f$ are portions of lines
65 : * passing through \f$P\f$.
66 : * - The rim of the disk (\f$\bar{x}^2+\bar{y}^2=1\f$ on \f$\bar{z}=-1\f$)
67 : * maps to the circle at \f$z = z_\mathrm{extent}\f$ on the sphere.
68 : *
69 : * Note that the \f$\bar{z}\f$ orientation is the *opposite* of
70 : * `CylindricalFlatEndcap`: here \f$\bar{z}=+1\f$ is the flat disk and
71 : * \f$\bar{z}=-1\f$ is the sphere. This reversal is necessary to keep the
72 : * Jacobian determinant positive, because the flat disk sits at a larger
73 : * physical \f$z\f$ than the far sphere wall.
74 : *
75 : * CylindricalFlatEndcapInterior is intended for the Pill domain, where
76 : * the filled-cylinder endcap blocks have their flat face (\f$\bar{z}=+1\f$)
77 : * at the end of the inner cubed-cylinder region (inside the outer domain
78 : * sphere) and their spherical face (\f$\bar{z}=-1\f$) on the outer domain
79 : * sphere.
80 : *
81 : * ### Requirements on map parameters
82 : *
83 : * - \f$P\f$ is sufficiently inside the sphere:
84 : * \f$|P - C_2| \leq 0.95\,R_2\f$.
85 : * - The flat disk lies inside the sphere, at least 5 % of \f$R_2\f$ below
86 : * \f$C_2^z\f$ and at most 95 % of \f$R_2\f$ below \f$C_2^z\f$:
87 : * \f[
88 : * C_2^z - 0.95\,R_2 \;\leq\; C_1^z \;\leq\; C_2^z - 0.05\,R_2.
89 : * \f]
90 : * - The flat disk is below the projection point: \f$C_1^z < P^z\f$.
91 : * - \f$z_\mathrm{extent}\f$ lies strictly on the sphere:
92 : * \f$|z_\mathrm{extent} - C_2^z| < R_2\f$.
93 : * - The focal parameter satisfies \f$t > 1\f$ (the sphere is beyond the
94 : * disk from \f$P\f$).
95 : * - The ratio \f$R_1/R_2\f$ is between 1/10 and 10.
96 : * - The entire flat disk rim lies strictly inside the sphere:
97 : * \f[
98 : * \sqrt{\bigl(\sqrt{(C_1^x-C_2^x)^2+(C_1^y-C_2^y)^2}+R_1\bigr)^2
99 : * +(C_1^z-C_2^z)^2} < R_2.
100 : * \f]
101 : * This is required for the inverse to be defined everywhere: if any rim
102 : * point were outside the sphere, the corresponding ray from \f$P\f$ would
103 : * hit the sphere before reaching the disk (\f$t_\mathrm{sphere}<1\f$),
104 : * violating the `source_is_between_focus_and_target` assumption.
105 : */
106 1 : class CylindricalFlatEndcapInterior {
107 : public:
108 0 : static constexpr size_t dim = 3;
109 :
110 : /*!
111 : * \brief Construct the map.
112 : *
113 : * \param center_one Center of the flat disk (\f$C_1\f$).
114 : * \param center_two Center of the outer sphere (\f$C_2\f$).
115 : * \param proj_center Projection point \f$P\f$.
116 : * \param z_sphere_extent z-coordinate of the rim circle where the spherical
117 : * face meets the adjacent hollow-cylinder block. This determines the flat
118 : * disk radius \f$R_1\f$ via the focal projection.
119 : * \param radius_two Radius of the outer sphere \f$R_2\f$.
120 : */
121 1 : CylindricalFlatEndcapInterior(const std::array<double, 3>& center_one,
122 : const std::array<double, 3>& center_two,
123 : const std::array<double, 3>& proj_center,
124 : double z_sphere_extent, double radius_two);
125 :
126 0 : CylindricalFlatEndcapInterior() = default;
127 0 : ~CylindricalFlatEndcapInterior() = default;
128 0 : CylindricalFlatEndcapInterior(CylindricalFlatEndcapInterior&&) = default;
129 0 : CylindricalFlatEndcapInterior(const CylindricalFlatEndcapInterior&) = default;
130 0 : CylindricalFlatEndcapInterior& operator=(
131 : const CylindricalFlatEndcapInterior&) = default;
132 0 : CylindricalFlatEndcapInterior& operator=(CylindricalFlatEndcapInterior&&) =
133 : default;
134 :
135 : template <typename T>
136 0 : std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
137 :
138 0 : std::optional<std::array<double, 3>> inverse(
139 : const std::array<double, 3>& target_coords) const;
140 :
141 : template <typename T>
142 0 : tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
143 : const std::array<T, 3>& source_coords) const;
144 :
145 : template <typename T>
146 0 : tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
147 : const std::array<T, 3>& source_coords) const;
148 :
149 : // NOLINTNEXTLINE(google-runtime-references)
150 0 : void pup(PUP::er& p);
151 :
152 0 : static bool is_identity() { return false; }
153 :
154 0 : static constexpr bool supports_hessian{false};
155 :
156 : private:
157 0 : friend bool operator==(const CylindricalFlatEndcapInterior& lhs,
158 : const CylindricalFlatEndcapInterior& rhs);
159 0 : FocallyLiftedMap<FocallyLiftedInnerMaps::FlatEndcap> impl_;
160 : };
161 :
162 0 : bool operator!=(const CylindricalFlatEndcapInterior& lhs,
163 : const CylindricalFlatEndcapInterior& rhs);
164 :
165 : } // namespace domain::CoordinateMaps
|
