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 <optional>
9 :
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Domain/CoordinateMaps/Distribution.hpp"
12 : #include "Utilities/Serialization/PupStlCpp17.hpp"
13 :
14 : /// \cond
15 : namespace PUP {
16 : class er;
17 : } // namespace PUP
18 : /// \endcond
19 :
20 : namespace domain::CoordinateMaps {
21 :
22 : /*!
23 : * \ingroup CoordinateMapsGroup
24 : * \brief Maps \f$\xi\f$ in the 1D interval \f$[A, B]\f$ to \f$x\f$ in the
25 : * interval \f$[a, b]\f$ according to a `domain::CoordinateMaps::Distribution`.
26 : *
27 : * \details The mapping takes a `domain::CoordinateMaps::Distribution` and
28 : * distributes the grid points accordingly.
29 : *
30 : * The formula for the mapping is, in case of a `Linear` distribution
31 : * \f{align}
32 : * x &= \frac{b}{B-A} (\xi-A) + \frac{a}{B-A} (B-\xi)\\
33 : * \xi &=\frac{B}{b-a} (x-a) + \frac{A}{b-a} (b-x)
34 : * \f}
35 : *
36 : * For every other distribution we use this linear mapping to map onto an
37 : * interval \f$[-1, 1]\f$, so we define:
38 : * \f{align}
39 : * f(\xi) &:= \frac{A+B-2\xi}{A-B} \in [-1, 1]\\
40 : * g(x) &:= \frac{a+b-2x}{a-b} \in [-1, 1]
41 : * \f}
42 : *
43 : * With this an `Equiangular` distribution is described by
44 : *
45 : * \f{align}
46 : * x &= \frac{a}{2} \left(1-\mathrm{tan}\left(\frac{\pi}{4}f(\xi)\right)\right)
47 : * + \frac{b}{2} \left(1+\mathrm{tan}\left(\frac{\pi}{4}f(\xi)\right)\right)\\
48 : * \xi &=
49 : * \frac{A}{2} \left(1-\frac{4}{\pi}\mathrm{arctan}\left(g(x)\right)\right) +
50 : * \frac{B}{2} \left(1+\frac{4}{\pi}\mathrm{arctan}\left(g(x)\right)\right)
51 : * \f}
52 : *
53 : * \note The equiangular distribution is intended to be used with the `Wedge`
54 : * map when equiangular coordinates are chosen for those maps. For more
55 : * information on this choice of coordinates, see the documentation for `Wedge`.
56 : *
57 : *
58 : * For both the `Logarithmic` and `Inverse` distribution, we first specify a
59 : * position for the singularity \f$c:=\f$`singularity_pos` outside the target
60 : * interval \f$[a, b]\f$.
61 : *
62 : * The `Logarithmic` distribution further requires the introduction of a
63 : * variable \f$\sigma\f$ dependent on whether \f$c\f$ is left (\f$ \sigma =
64 : * 1\f$) or right (\f$ \sigma = -1\f$) of the target interval. With this, the
65 : * `Logarithmic` distribution is described by
66 : *
67 : * \f{align}
68 : * x &= \sigma\, {\rm exp}\left(\frac{\ln(b-c)+\ln(a-c)}{2} + f(\xi)
69 : * \frac{\ln(b-c)-\ln(a-c)}{2}\right) + c\\
70 : * \xi &= \frac{B-A}{2}\frac{2\ln(\sigma [x-c]) - \ln(b-c) - \ln(a-c)}{\ln(b-c)
71 : * - \ln(a-c)} + \frac{B+A}{2} \f}
72 : *
73 : * and the `Inverse` distribution by
74 : *
75 : * \f{align}
76 : * x &= \frac{2(a-c)(b-c)}{a+b-2c+(a-b)f(\xi)} + c\\
77 : * \xi &= -\frac{A-B}{2(a-b)} \left(\frac{2(a-c)(b-c)}{x-c} - (a-c) -
78 : * (b-c)\right) + \frac{A+B}{2}
79 : * \f}
80 : */
81 1 : class Interval {
82 : public:
83 0 : static constexpr size_t dim = 1;
84 :
85 0 : Interval(double A, double B, double a, double b, Distribution distribution,
86 : std::optional<double> singularity_pos = std::nullopt);
87 :
88 0 : Interval() = default;
89 :
90 : template <typename T>
91 0 : std::array<T, 1> operator()(const std::array<T, 1>& source_coords) const;
92 :
93 0 : std::optional<std::array<double, 1>> inverse(
94 : const std::array<double, 1>& target_coords) const;
95 :
96 : template <typename T>
97 0 : tnsr::Ij<T, 1, Frame::NoFrame> jacobian(
98 : const std::array<T, 1>& source_coords) const;
99 :
100 : template <typename T>
101 0 : tnsr::Ij<T, 1, Frame::NoFrame> inv_jacobian(
102 : const std::array<T, 1>& source_coords) const;
103 :
104 : // NOLINTNEXTLINE(google-runtime-references)
105 0 : void pup(PUP::er& p);
106 :
107 0 : bool is_identity() const { return is_identity_; }
108 :
109 0 : static constexpr bool supports_hessian{true};
110 :
111 : private:
112 0 : friend bool operator==(const Interval& lhs, const Interval& rhs);
113 :
114 0 : double A_{std::numeric_limits<double>::signaling_NaN()};
115 0 : double B_{std::numeric_limits<double>::signaling_NaN()};
116 0 : double a_{std::numeric_limits<double>::signaling_NaN()};
117 0 : double b_{std::numeric_limits<double>::signaling_NaN()};
118 0 : Distribution distribution_{Distribution::Linear};
119 0 : std::optional<double> singularity_pos_{std::nullopt};
120 0 : bool is_identity_{false};
121 : };
122 :
123 0 : inline bool operator!=(const CoordinateMaps::Interval& lhs,
124 : const CoordinateMaps::Interval& rhs) {
125 : return not(lhs == rhs);
126 : }
127 :
128 : } // namespace domain::CoordinateMaps
|
