Line data Source code
1 1 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : /// \file 5 : /// Defines the class CylindricalFlatSide. 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/FocallyLiftedFlatSide.hpp" 16 : #include "Domain/CoordinateMaps/FocallyLiftedMap.hpp" 17 : #include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp" 18 : 19 : /// \cond 20 : namespace PUP { 21 : class er; 22 : } // namespace PUP 23 : /// \endcond 24 : 25 : namespace domain::CoordinateMaps { 26 : 27 : /*! 28 : * \ingroup CoordinateMapsGroup 29 : * 30 : * \brief Map from 3D unit right cylindrical shell to a volume that connects 31 : * a portion of an annulus to a portion of a spherical surface. 32 : * 33 : * \image html CylindricalFlatSide.svg "A cylinder maps to the shaded region." 34 : * 35 : * \details Consider a 2D annulus in 3D space that is normal to the 36 : * \f$z\f$ axis and has (3D) center \f$C_1\f$, inner radius 37 : * \f$R_\mathrm{in}\f$ and outer radius \f$R_\mathrm{out}\f$ 38 : * Also consider a sphere with center \f$C_2\f$, and radius \f$R_2\f$. 39 : * Also let there be a projection point \f$P\f$. 40 : * 41 : * CylindricalFlatSide maps a 3D unit right cylindrical shell (with 42 : * coordinates \f$(\bar{x},\bar{y},\bar{z})\f$ such that 43 : * \f$-1\leq\bar{z}\leq 1\f$ and \f$1 \leq \bar{x}^2+\bar{y}^2 \leq 44 : * 4\f$) to the shaded area in the figure above (with coordinates 45 : * \f$(x,y,z)\f$). The "bottom" of the cylinder \f$\bar{z}=-1\f$ is 46 : * mapped to the interior of the annulus with radii 47 : * \f$R_\mathrm{in}\f$ and \f$R_\mathrm{out}\f$. Curves of constant 48 : * \f$(\bar{x},\bar{y})\f$ are mapped to portions of lines that pass 49 : * through \f$P\f$. Along each of these curves, \f$\bar{z}=-1\f$ is 50 : * mapped to a point inside the annulus and \f$\bar{z}=+1\f$ is mapped to a 51 : * point on the sphere. 52 : * 53 : * CylindricalFlatSide is intended to be composed with Wedge2D maps to 54 : * construct a portion of a cylindrical domain for a binary system. 55 : * 56 : * CylindricalFlatSide is described briefly in the Appendix of 57 : * \cite Buchman:2012dw. 58 : * CylindricalFlatSide is used to construct the blocks labeled 'ME 59 : * cylinder' in Figure 20 of that paper. 60 : * 61 : * CylindricalFlatSide is implemented using `FocallyLiftedMap` 62 : * and `FocallyLiftedInnerMaps::FlatSide`; see those classes for 63 : * details. 64 : * 65 : * ### Restrictions on map parameters. 66 : * 67 : * We demand that: 68 : * - The sphere is at a larger value of \f$z\f$ (plus 5 69 : * percent of the sphere radius) than the plane containing the 70 : * annulus. 71 : * - The projection point \f$z_\mathrm{P}\f$ is 72 : * inside the sphere and more than 15 percent away from the boundary 73 : * of the sphere. 74 : * - The center of the annulus is contained in the circle that results from 75 : * projecting the sphere into the \f$xy\f$ plane. 76 : * - The outer radius of the annulus is larger than 5 percent of the distance 77 : * between the center of the annulus and the projection point. 78 : * - The inner radius of the annulus is less than 95 percent of the outer 79 : * radius, larger than 5 percent of the outer radius, and larger than one 80 : * percent of the distance between the center of the annulus and the 81 : * projection point. The last condition means that the angle subtended by 82 : * the inner radius with respect to the projection point is not too small. 83 : * 84 : * It is possible to construct a valid map without these assumptions, 85 : * but some of these assumptions simplify the code and others eliminate 86 : * edge cases where Jacobians become large or small. 87 : * 88 : */ 89 1 : class CylindricalFlatSide { 90 : public: 91 0 : static constexpr size_t dim = 3; 92 0 : CylindricalFlatSide(const std::array<double, 3>& center_one, 93 : const std::array<double, 3>& center_two, 94 : const std::array<double, 3>& proj_center, 95 : const double inner_radius, const double outer_radius, 96 : const double radius_two); 97 : 98 0 : CylindricalFlatSide() = default; 99 0 : ~CylindricalFlatSide() = default; 100 0 : CylindricalFlatSide(CylindricalFlatSide&&) = default; 101 0 : CylindricalFlatSide(const CylindricalFlatSide&) = default; 102 0 : CylindricalFlatSide& operator=(const CylindricalFlatSide&) = default; 103 0 : CylindricalFlatSide& operator=(CylindricalFlatSide&&) = default; 104 : 105 : template <typename T> 106 0 : std::array<tt::remove_cvref_wrap_t<T>, 3> operator()( 107 : const std::array<T, 3>& source_coords) const; 108 : 109 0 : std::optional<std::array<double, 3>> inverse( 110 : const std::array<double, 3>& target_coords) const; 111 : 112 : template <typename T> 113 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian( 114 : const std::array<T, 3>& source_coords) const; 115 : 116 : template <typename T> 117 0 : tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian( 118 : const std::array<T, 3>& source_coords) const; 119 : 120 : // NOLINTNEXTLINE(google-runtime-references) 121 0 : void pup(PUP::er& p); 122 : 123 0 : static bool is_identity() { return false; } 124 : 125 : private: 126 0 : friend bool operator==(const CylindricalFlatSide& lhs, 127 : const CylindricalFlatSide& rhs); 128 0 : FocallyLiftedMap<FocallyLiftedInnerMaps::FlatSide> impl_; 129 : }; 130 0 : bool operator!=(const CylindricalFlatSide& lhs, const CylindricalFlatSide& rhs); 131 : 132 : } // namespace domain::CoordinateMaps