Wedge3D.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <boost/optional.hpp>
8 #include <cstddef>
9 #include <limits>
10 
12 #include "Domain/OrientationMap.hpp"
13 #include "Utilities/TypeTraits.hpp"
14 
15 namespace PUP {
16 class er;
17 } // namespace PUP
18 
19 namespace domain {
20 namespace CoordinateMaps {
21 
22 /*!
23  * \ingroup CoordinateMapsGroup
24  *
25  * \brief Three dimensional map from the cube to a wedge.
26  * \image html Shell.png "A shell can be constructed out of six wedges."
27  *
28  * \details The mapping that goes from a reference cube to a three-dimensional
29  * wedge centered on a coordinate axis covering a volume between an inner
30  * surface and outer surface. Each surface can be given a curvature
31  * between flat (a sphericity of 0) or spherical (a sphericity of 1).
32  *
33  * The first two logical coordinates correspond to the two angular coordinates,
34  * and the third to the radial coordinate.
35  *
36  * The Wedge3D map is constructed by linearly interpolating between a bulged
37  * face of radius `radius_of_inner_surface` to a bulged face of
38  * radius `radius_of_outer_surface`, where the radius of each bulged face
39  * is defined to be the radius of the sphere circumscribing the bulge.
40  *
41  * We make a choice here as to whether we wish to use the logical coordinates
42  * parameterizing these surface as they are, in which case we have the
43  * equidistant choice of coordinates, or whether to apply a tangent map to them
44  * which leads us to the equiangular choice of coordinates. In terms of the
45  * logical coordinates, the equiangular coordinates are:
46  *
47  * \f[\textrm{equiangular xi} : \Xi(\xi) = \textrm{tan}(\xi\pi/4)\f]
48  *
49  * \f[\textrm{equiangular eta} : \mathrm{H}(\eta) = \textrm{tan}(\eta\pi/4)\f]
50  *
51  * With derivatives:
52  *
53  * \f[\Xi'(\xi) = \frac{\pi}{4}(1+\Xi^2)\f]
54  *
55  * \f[\mathrm{H}'(\eta) = \frac{\pi}{4}(1+\mathrm{H}^2)\f]
56  *
57  * The equidistant coordinates are:
58  *
59  * \f[ \textrm{equidistant xi} : \Xi = \xi\f]
60  *
61  * \f[ \textrm{equidistant eta} : \mathrm{H} = \eta\f]
62  *
63  * with derivatives:
64  *
65  * <center>\f$\Xi'(\xi) = 1\f$, and \f$\mathrm{H}'(\eta) = 1\f$</center>
66  *
67  * We also define the variable \f$\rho\f$, given by:
68  *
69  * \f[\textrm{rho} : \rho = \sqrt{1+\Xi^2+\mathrm{H}^2}\f]
70  *
71  * ### The Spherical Face Map
72  * The surface map for the spherical face of radius \f$R\f$ lying in the
73  * \f$+z\f$
74  * direction in either choice of coordinates is then given by:
75  *
76  * \f[\vec{\sigma}_{spherical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\f]
77  * Where
78  * \f[
79  * \vec{x}(\xi,\eta) =
80  * \begin{bmatrix}
81  * x(\xi,\eta)\\
82  * y(\xi,\eta)\\
83  * z(\xi,\eta)\\
84  * \end{bmatrix} = \frac{R}{\rho}
85  * \begin{bmatrix}
86  * \Xi\\
87  * \mathrm{H}\\
88  * 1\\
89  * \end{bmatrix}\f]
90  *
91  * ### The Bulged Face Map
92  * The bulged surface is itself constructed by linearly interpolating between
93  * a cubical face and a spherical face. The surface map for the cubical face
94  * of side length \f$2L\f$ lying in the \f$+z\f$ direction is given by:
95  *
96  * \f[\vec{\sigma}_{cubical}: \vec{\xi} \rightarrow \vec{x}(\vec{\xi})\f]
97  * Where
98  * \f[
99  * \vec{x}(\xi,\eta) =
100  * \begin{bmatrix}
101  * x(\xi,\eta)\\
102  * y(\xi,\eta)\\
103  * L\\
104  * \end{bmatrix} = L
105  * \begin{bmatrix}
106  * \Xi\\
107  * \mathrm{H}\\
108  * 1\\
109  * \end{bmatrix}\f]
110  *
111  * To construct the bulged map we interpolate between this cubical face map
112  * and a spherical face map of radius \f$R\f$, with the
113  * interpolation parameter being \f$s\f$. The surface map for the bulged face
114  * lying in the \f$+z\f$ direction is then given by:
115  *
116  * \f[\vec{\sigma}_{bulged}(\xi,\eta) = {(1-s)L + \frac{sR}{\rho}}
117  * \begin{bmatrix}
118  * \Xi\\
119  * \mathrm{H}\\
120  * 1\\
121  * \end{bmatrix}\f]
122  *
123  * We constrain L by demanding that the spherical face circumscribe the cube.
124  * With this condition, we have \f$L = R/\sqrt3\f$.
125  * \note This differs from the choice in SpEC where it is demanded that the
126  * surfaces touch at the center, which leads to \f$L = R\f$.
127  *
128  * ### The Full Volume Map
129  * The final map for the wedge which lies along the \f$+z\f$ is obtained by
130  * interpolating between the two surfaces with the
131  * interpolation parameter being the logical coordinate \f$\zeta\f$. This
132  * results in:
133  *
134  * \f[\vec{x}(\xi,\eta,\zeta) =
135  * \frac{1}{2}\left\{(1-\zeta)\Big[(1-s_{inner})\frac{R_{inner}}{\sqrt 3}
136  * + s_{inner}\frac{R_{inner}}{\rho}\Big] +
137  * (1+\zeta)\Big[(1-s_{outer})\frac{R_{outer}}{\sqrt 3} +s_{outer}
138  * \frac{R_{outer}}{\rho}\Big] \right\}\begin{bmatrix}
139  * \Xi\\
140  * \mathrm{H}\\
141  * 1\\
142  * \end{bmatrix}\f]
143  *
144  * We will define the variables \f$F(\zeta)\f$ and \f$S(\zeta)\f$, the frustum
145  * and sphere factors: \f[F(\zeta) = F_0 + F_1\zeta\f] \f[S(\zeta) = S_0 +
146  * S_1\zeta\f]
147  * Where \f{align*}F_0 &= \frac{1}{2} \big\{ (1-s_{outer})R_{outer} +
148  * (1-s_{inner})R_{inner}\big\}\\
149  * F_1 &= \partial_{\zeta} F = \frac{1}{2} \big\{ (1-s_{outer})R_{outer} -
150  * (1-s_{inner})R_{inner}\big\}\\
151  * S_0 &= \frac{1}{2} \big\{ s_{outer}R_{outer} + s_{inner}R_{inner}\big\}\\
152  * S_1 &= \partial_{\zeta} S = \frac{1}{2} \big\{ s_{outer}R_{outer} -
153  * s_{inner}R_{inner}\big\}\f}
154  *
155  * The map can then be rewritten as:
156  * \f[\vec{x}(\xi,\eta,\zeta) = \left\{\frac{F(\zeta)}{\sqrt 3} +
157  * \frac{S(\zeta)}{\rho}\right\}\begin{bmatrix}
158  * \Xi\\
159  * \mathrm{H}\\
160  * 1\\
161  * \end{bmatrix}\f]
162  *
163  * We provide some common derivatives:
164  * \f[\partial_{\xi}z = \frac{-S(\zeta)\Xi\Xi'}{\rho^3}\f]
165  * \f[\partial_{\eta}z = \frac{-S(\zeta)\mathrm{H}\mathrm{H}'}{\rho^3}\f]
166  * \f[\partial_{\zeta}z = \frac{F'}{\sqrt 3} + \frac{S'}{\rho}\f]
167  * The Jacobian then is: \f[J =
168  * \begin{bmatrix}
169  * \Xi'z + \Xi\partial_{\xi}z & \Xi\partial_{\eta}z & \Xi\partial_{\zeta}z \\
170  * \mathrm{H}\partial_{\xi}z & \mathrm{H}'z +
171  * \mathrm{H}\partial_{\eta}z & \mathrm{H}\partial_{\zeta}z\\
172  * \partial_{\xi}z&\partial_{\eta}z &\partial_{\zeta}z \\
173  * \end{bmatrix}
174  * \f]
175  *
176  * A common factor that shows up in the inverse jacobian is:
177  * \f[ T:= \frac{S(\zeta)}{(\partial_{\zeta}z)\rho^3}\f]
178  *
179  * The inverse Jacobian then is: \f[J^{-1} =
180  * \frac{1}{z}\begin{bmatrix}
181  * \Xi'^{-1} & 0 & -\Xi\Xi'^{-1}\\
182  * 0 & \mathrm{H}'^{-1} & -\mathrm{H}\mathrm{H}'^{-1}\\
183  * T\Xi &
184  * T\mathrm{H} &
185  * T + F(\partial_{\zeta}z)^{-1}/\sqrt 3\\
186  * \end{bmatrix}
187  * \f]
188  *
189  * ### Changing the radial distribution of the gridpoints
190  * By default, Wedge3D linearly distributes its gridpoints in the radial
191  * direction. An exponential distribution of gridpoints can be obtained by
192  * linearly interpolating in the logarithm of the radius, in order to obtain
193  * a relatively higher resolution at smaller radii. Since this is a radial
194  * rescaling of Wedge3D, this option is only supported for fully spherical
195  * wedges with `sphericity_inner` = `sphericity_outer` = 1.
196  *
197  * The linear interpolation done is:
198  * \f[
199  * \log r = \frac{1-\zeta}{2}\log R_{inner} +
200  * \frac{1+\zeta}{2}\log R_{outer}
201  * \f]
202  *
203  * The map then is:
204  * \f[\vec{x}(\xi,\eta,\zeta) =
205  * \frac{\sqrt{R_{inner}^{1-\zeta}R_{outer}^{1+\zeta}}}{\rho}\begin{bmatrix}
206  * \Xi\\
207  * \mathrm{H}\\
208  * 1\\
209  * \end{bmatrix}\f]
210  *
211  * The jacobian simplifies similarly.
212  *
213  */
214 class Wedge3D {
215  public:
216  static constexpr size_t dim = 3;
217  enum class WedgeHalves {
218  /// Use the entire wedge
219  Both,
220  /// Use only the upper logical half
221  UpperOnly,
222  /// Use only the lower logical half
223  LowerOnly
224  };
225 
226  /*!
227  * Constructs a 3D wedge.
228  * \param radius_inner Distance from the origin to one of the
229  * corners which lie on the inner surface.
230  * \param radius_outer Distance from the origin to one of the
231  * corners which lie on the outer surface.
232  * \param orientation_of_wedge The orientation of the desired wedge relative
233  * to the orientation of the default wedge which is a wedge that has its
234  * curved surfaces pierced by the upper-z axis. The logical xi and eta
235  * coordinates point in the cartesian x and y directions, respectively.
236  * \param sphericity_inner Value between 0 and 1 which determines
237  * whether the inner surface is flat (value of 0), spherical (value of 1) or
238  * somewhere in between
239  * \param sphericity_outer Value between 0 and 1 which determines
240  * whether the outer surface is flat (value of 0), spherical (value of 1) or
241  * somewhere in between
242  * \param with_equiangular_map Determines whether to apply a tangent function
243  * mapping to the logical coordinates (for `true`) or not (for `false`).
244  * \param halves_to_use Determines whether to construct a full wedge or only
245  * half a wedge. If constructing only half a wedge, the resulting shape has a
246  * face normal to the x direction (assuming default OrientationMap). If
247  * constructing half a wedge, an intermediate affine map is applied to the
248  * logical xi coordinate such that the interval [-1,1] is mapped to the
249  * corresponding logical half of the wedge. For example, if `UpperOnly` is
250  * specified, [-1,1] is mapped to [0,1], and if `LowerOnly` is specified,
251  * [-1,1] is mapped to [-1,0]. The case of `Both` means a full wedge, with no
252  * intermediate map applied. In all cases, the logical points returned by the
253  * inverse map will lie in the range [-1,1] in each dimension. Half wedges are
254  * currently only useful in constructing domains for binary systems.
255  * \param with_logarithmic_map Determines whether to apply an exponential
256  * function mapping to the "sphere factor", the effect of which is to
257  * distribute the radial gridpoints logarithmically in physical space.
258  */
259  Wedge3D(double radius_inner, double radius_outer,
260  OrientationMap<3> orientation_of_wedge, double sphericity_inner,
261  double sphericity_outer, bool with_equiangular_map,
262  WedgeHalves halves_to_use = WedgeHalves::Both,
263  bool with_logarithmic_map = false) noexcept;
264 
265  Wedge3D() = default;
266  ~Wedge3D() = default;
267  Wedge3D(Wedge3D&&) = default;
268  Wedge3D(const Wedge3D&) = default;
269  Wedge3D& operator=(const Wedge3D&) = default;
270  Wedge3D& operator=(Wedge3D&&) = default;
271 
272  template <typename T>
274  const std::array<T, 3>& source_coords) const noexcept;
275 
276  /// For a \f$+z\f$-oriented `Wedge3D`, returns invalid if \f$z<=0\f$
277  /// or if \f$(x,y,z)\f$ is on or outside the cone defined
278  /// by \f$(x^2/z^2 + y^2/z^2+1)^{1/2} = -S/F\f$, where
279  /// \f$S = \frac{1}{2}(s_1 r_1 - s_0 r_0)\f$ and
280  /// \f$F = \frac{1}{2\sqrt{3}}((1-s_1) r_1 - (1-s_0) r_0)\f$.
281  /// Here \f$s_0,s_1\f$ and \f$r_0,r_1\f$ are the specified sphericities
282  /// and radii of the inner and outer \f$z\f$ surfaces. The map is singular on
283  /// the cone and on the xy plane.
284  boost::optional<std::array<double, 3>> inverse(
285  const std::array<double, 3>& target_coords) const noexcept;
286 
287  template <typename T>
288  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian(
289  const std::array<T, 3>& source_coords) const noexcept;
290 
291  template <typename T>
292  tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian(
293  const std::array<T, 3>& source_coords) const noexcept;
294 
295  // clang-tidy: google runtime references
296  void pup(PUP::er& p) noexcept; // NOLINT
297 
298  bool is_identity() const noexcept { return false; }
299 
300  private:
301  // factors out calculation of z needed for mapping and jacobian
302  template <typename T>
303  tt::remove_cvref_wrap_t<T> default_physical_z(const T& zeta,
304  const T& one_over_rho) const
305  noexcept;
306  friend bool operator==(const Wedge3D& lhs, const Wedge3D& rhs) noexcept;
307 
308  double radius_inner_{std::numeric_limits<double>::signaling_NaN()};
309  double radius_outer_{std::numeric_limits<double>::signaling_NaN()};
310  OrientationMap<3> orientation_of_wedge_{};
311  double sphericity_inner_{std::numeric_limits<double>::signaling_NaN()};
312  double sphericity_outer_{std::numeric_limits<double>::signaling_NaN()};
313  bool with_equiangular_map_ = false;
314  WedgeHalves halves_to_use_ = WedgeHalves::Both;
315  bool with_logarithmic_map_ = false;
316  double scaled_frustum_zero_{std::numeric_limits<double>::signaling_NaN()};
317  double sphere_zero_{std::numeric_limits<double>::signaling_NaN()};
318  double scaled_frustum_rate_{std::numeric_limits<double>::signaling_NaN()};
319  double sphere_rate_{std::numeric_limits<double>::signaling_NaN()};
320 };
321 bool operator!=(const Wedge3D& lhs, const Wedge3D& rhs) noexcept;
322 } // namespace CoordinateMaps
323 } // namespace domain
Definition: Strahlkorper.hpp:14
Three dimensional map from the cube to a wedge.
Definition: Wedge3D.hpp:214
Definition: BlockId.hpp:16
T signaling_NaN(T... args)
WedgeHalves
Definition: Wedge3D.hpp:217
Defines a list of useful type aliases for tensors.
Represents an index that is not in a known frame, e.g. some internal intermediate frame that is irrel...
Definition: IndexType.hpp:48
Defines type traits, some of which are future STL type_traits header.