Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines the class CylindricalFlatEndcap.
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/FocallyLiftedFlatEndcap.hpp"
16 : #include "Domain/CoordinateMaps/FocallyLiftedMap.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 3D unit right cylinder to a volume that connects
30 : * a portion of a circle to a portion of a spherical surface.
31 : *
32 : * \image html CylindricalFlatEndcap.svg "A cylinder maps to the shaded region."
33 : *
34 : * \details Consider a 2D circle in 3D space that is normal to the
35 : * \f$z\f$ axis and has (3D) center \f$C_1\f$ and radius \f$R_1\f$.
36 : * Also consider a sphere with center \f$C_2\f$, and radius \f$R_2\f$.
37 : * Also let there be a projection point \f$P\f$.
38 : *
39 : * CylindricalFlatEndcap maps a 3D unit right cylinder (with coordinates
40 : * \f$(\bar{x},\bar{y},\bar{z})\f$ such that \f$-1\leq\bar{z}\leq 1\f$
41 : * and \f$\bar{x}^2+\bar{y}^2 \leq 1\f$) to the shaded area
42 : * in the figure above (with coordinates \f$(x,y,z)\f$). The "bottom"
43 : * of the cylinder \f$\bar{z}=-1\f$ is mapped to the interior of the
44 : * circle of radius \f$R_1\f$. Curves of constant
45 : * \f$(\bar{x},\bar{y})\f$ are mapped to portions of lines that pass
46 : * through \f$P\f$. Along each of these curves, \f$\bar{z}=-1\f$ is
47 : * mapped to a point on the circle and \f$\bar{z}=+1\f$ is mapped to a
48 : * point on the sphere.
49 : *
50 : * CylindricalFlatEndcap is intended to be composed with Wedge2D maps to
51 : * construct a portion of a cylindrical domain for a binary system.
52 : *
53 : * CylindricalFlatEndcap is described briefly in the Appendix of
54 : * \cite Buchman:2012dw.
55 : * CylindricalFlatEndcap is used to construct the blocks labeled 'MA
56 : * wedge' and 'MB wedge' in Figure 20 of that paper.
57 : *
58 : * CylindricalFlatEndcap is implemented using `FocallyLiftedMap`
59 : * and `FocallyLiftedInnerMaps::FlatEndcap`; see those classes for
60 : * details.
61 : *
62 : * ### Restrictions on map parameters.
63 : *
64 : * The following restrictions are made so that the map is not singular
65 : * or close to singular. It is possible to construct a valid map
66 : * without these assumptions, but the assumptions simplify the code and
67 : * avoid problematic edge cases, and the expected use cases obey
68 : * these restrictions.
69 : *
70 : * We demand that
71 : * - The plane containing the circle is below (i.e. at a smaller value
72 : * of \f$z\f$ than) the sphere, by an amount at least 5% of the sphere
73 : * radius \f$R_2\f$ but not more than 5 times the sphere radius \f$R_2\f$.
74 : * - \f$P\f$ is inside the sphere but not too close to its surface;
75 : * specifically, we demand that \f$|P-C_2|\leq 0.95 R_2\f$.
76 : * - The ratio \f$R_1/R_2\f$ is between 10 and 1/10, inclusive.
77 : * - The x and y components of \f$C_2-C_1\f$ both have magnitudes
78 : * smaller than or equal to \f$R_1+R_2\f$.
79 : *
80 : */
81 1 : class CylindricalFlatEndcap {
82 : public:
83 0 : static constexpr size_t dim = 3;
84 0 : CylindricalFlatEndcap(const std::array<double, 3>& center_one,
85 : const std::array<double, 3>& center_two,
86 : const std::array<double, 3>& proj_center,
87 : double radius_one, double radius_two);
88 :
89 0 : CylindricalFlatEndcap() = default;
90 0 : ~CylindricalFlatEndcap() = default;
91 0 : CylindricalFlatEndcap(CylindricalFlatEndcap&&) = default;
92 0 : CylindricalFlatEndcap(const CylindricalFlatEndcap&) = default;
93 0 : CylindricalFlatEndcap& operator=(const CylindricalFlatEndcap&) = default;
94 0 : CylindricalFlatEndcap& operator=(CylindricalFlatEndcap&&) = default;
95 :
96 : template <typename T>
97 0 : std::array<T, 3> operator()(const std::array<T, 3>& source_coords) const;
98 :
99 0 : std::optional<std::array<double, 3>> inverse(
100 : const std::array<double, 3>& target_coords) const;
101 :
102 : template <typename T>
103 0 : tnsr::Ij<T, 3, Frame::NoFrame> jacobian(
104 : const std::array<T, 3>& source_coords) const;
105 :
106 : template <typename T>
107 0 : tnsr::Ij<T, 3, Frame::NoFrame> inv_jacobian(
108 : const std::array<T, 3>& source_coords) const;
109 :
110 : // NOLINTNEXTLINE(google-runtime-references)
111 0 : void pup(PUP::er& p);
112 :
113 0 : static bool is_identity() { return false; }
114 :
115 0 : static constexpr bool supports_hessian{false};
116 :
117 : private:
118 0 : friend bool operator==(const CylindricalFlatEndcap& lhs,
119 : const CylindricalFlatEndcap& rhs);
120 0 : FocallyLiftedMap<FocallyLiftedInnerMaps::FlatEndcap> impl_;
121 : };
122 0 : bool operator!=(const CylindricalFlatEndcap& lhs,
123 : const CylindricalFlatEndcap& rhs);
124 :
125 : } // namespace domain::CoordinateMaps
|
