Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines the class template Mesh.
6 :
7 : #pragma once
8 :
9 : #include <array>
10 : #include <cstddef>
11 :
12 : #include "DataStructures/Index.hpp"
13 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
14 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
15 : #include "Options/String.hpp"
16 : #include "Utilities/ErrorHandling/Assert.hpp"
17 : #include "Utilities/Gsl.hpp"
18 : #include "Utilities/Requires.hpp"
19 : #include "Utilities/TypeTraits.hpp"
20 : #include "Utilities/TypeTraits/IsInteger.hpp"
21 :
22 : /// \cond
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : /*!
29 : * \ingroup DataStructuresGroup
30 : * \brief Holds the number of grid points, basis, and quadrature in each
31 : * direction of the computational grid.
32 : *
33 : * \details A mesh encapsulates all information necessary to construct the
34 : * placement of grid points in the computational domain. It does so through a
35 : * choice of basis functions, quadrature and number of points \f$N\f$ in each
36 : * dimension. The grid points are the associated collocation points and can be
37 : * obtained from Spectral::collocation_points(const Mesh<1>&):
38 : *
39 : * \snippet Test_Spectral.cpp get_points_for_mesh
40 : *
41 : * A simulated physical field can be represented by a DataVector of length
42 : * number_of_grid_points() that holds the field value on each point of
43 : * the computational grid. These values are identical to the field's nodal
44 : * expansion coefficients. They approximate the field by a polynomial of degree
45 : * \f$p=N-1\f$ through a linear combination of Lagrange polynomials.
46 : *
47 : * \note Because we use a compact bit representation for the mesh, the number
48 : * of grid points/extents must be fewer than 256 per dimension.
49 : *
50 : * \tparam Dim the number of dimensions of the computational grid.
51 : */
52 : template <size_t Dim>
53 1 : class Mesh {
54 : public:
55 0 : static constexpr size_t dim = Dim;
56 :
57 0 : struct Extents {
58 0 : using type = size_t;
59 0 : static constexpr Options::String help = {
60 : "The number of collocation points per dimension"};
61 : };
62 :
63 0 : struct Basis {
64 0 : using type = Spectral::Basis;
65 0 : static constexpr Options::String help = {
66 : "The choice of spectral basis to compute the collocation points"};
67 : };
68 :
69 0 : struct Quadrature {
70 0 : using type = Spectral::Quadrature;
71 0 : static constexpr Options::String help = {
72 : "The choice of quadrature to compute the collocation points"};
73 : };
74 :
75 0 : using options = tmpl::list<Extents, Basis, Quadrature>;
76 :
77 0 : static constexpr Options::String help =
78 : "Holds the number of grid points, basis, and quadrature in each "
79 : "direction of the computational grid. "
80 : "A mesh encapsulates all information necessary to construct the "
81 : "placement of grid points in the computational domain. It does so "
82 : "through a choice of basis functions, quadrature and number of points "
83 : "in each dimension.";
84 :
85 0 : Mesh() = default;
86 :
87 : /*!
88 : * \brief Construct a computational grid with the same number of collocation
89 : * points in each dimension.
90 : *
91 : * \param isotropic_extents The number of collocation points in each
92 : * dimension.
93 : * \param basis The choice of spectral basis to compute the
94 : * collocation points
95 : * \param quadrature The choice of quadrature to compute
96 : * the collocation points
97 : *
98 : * \note Because a `Mesh<0>` extends over no dimensions, it has 1 grid point
99 : * independent of the value of `isotropic_extents`, and the `basis` and
100 : * `quadrature` are ignored.
101 : */
102 1 : Mesh(size_t isotropic_extents, Spectral::Basis basis,
103 : Spectral::Quadrature quadrature);
104 :
105 : /*!
106 : * \brief Construct a computational grid where each dimension can have a
107 : * different number of collocation points.
108 : *
109 : * \param extents The number of collocation points per dimension
110 : * \param basis The choice of spectral basis to compute the
111 : * collocation points
112 : * \param quadrature The choice of quadrature to compute
113 : * the collocation points
114 : *
115 : * \note A `Mesh<0>` extends over no dimensions so the `basis` and
116 : * `quadrature` are ignored.
117 : */
118 1 : Mesh(const std::array<size_t, Dim>& extents, Spectral::Basis basis,
119 : Spectral::Quadrature quadrature);
120 :
121 : /*!
122 : * \brief Construct a computational grid where each dimension can have both a
123 : * different number and placement of collocation points.
124 : *
125 : * \param extents The number of collocation points per dimension
126 : * \param bases The choice of spectral bases to compute the
127 : * collocation points per dimension
128 : * \param quadratures The choice of quadratures to compute
129 : * the collocation points per dimension
130 : */
131 1 : Mesh(const std::array<size_t, Dim>& extents,
132 : const std::array<Spectral::Basis, Dim>& bases,
133 : const std::array<Spectral::Quadrature, Dim>& quadratures);
134 :
135 : /*!
136 : * \brief The number of grid points in each dimension of the grid.
137 : */
138 1 : Index<Dim> extents() const;
139 :
140 : /*!
141 : * \brief The number of grid points in dimension \p d of the grid
142 : * (zero-indexed).
143 : */
144 1 : size_t extents(size_t d) const;
145 :
146 : /*!
147 : * \brief The total number of grid points in all dimensions.
148 : *
149 : * \details `DataVector`s that represent field values on the grid have this
150 : * many entries.
151 : *
152 : * \note A zero-dimensional mesh has one grid point, since it is the slice
153 : * through a one-dimensional mesh (a line).
154 : */
155 1 : size_t number_of_grid_points() const;
156 :
157 : /*!
158 : * \brief Returns the 1-dimensional index corresponding to the `Dim`
159 : * dimensional `index`.
160 : *
161 : * The first dimension varies fastest.
162 : *
163 : * \see collapsed_index()
164 : */
165 1 : size_t storage_index(const Index<Dim>& index) const;
166 :
167 : /*!
168 : * \brief The basis chosen in each dimension of the grid.
169 : */
170 1 : std::array<Spectral::Basis, Dim> basis() const;
171 :
172 : /*!
173 : * \brief The basis chosen in dimension \p d of the grid (zero-indexed).
174 : */
175 1 : Spectral::Basis basis(size_t d) const;
176 :
177 : /*!
178 : * \brief The quadrature chosen in each dimension of the grid.
179 : */
180 1 : std::array<Spectral::Quadrature, Dim> quadrature() const;
181 :
182 : /*!
183 : * \brief The quadrature chosen in dimension \p d of the grid (zero-indexed).
184 : */
185 1 : Spectral::Quadrature quadrature(size_t d) const;
186 :
187 : /*!
188 : * \brief Returns a Mesh with dimension \p d removed (zero-indexed).
189 : *
190 : * \see slice_through()
191 : */
192 : // clang-tidy: incorrectly reported redundancy in template expression
193 : template <size_t N = Dim, Requires<(N > 0 and N == Dim)> = nullptr> // NOLINT
194 1 : Mesh<Dim - 1> slice_away(size_t d) const;
195 :
196 : /*!
197 : * \brief Returns the interface Mesh for a DG boundary in direction \p d.
198 : *
199 : * \details Unlike `slice_away`, this function converts ZernikeB3 (ZernikeB2)
200 : * angular bases to SphericalHarmonic (Fourier) when slicing away the radial
201 : * dimension of a B3 ball (B2 disk), reflecting the fact that the outer
202 : * boundary of a B3 ball (B2 disk) is geometrically an S2 (S1) and shares its
203 : * collocation points with SphericalHarmonic (Fourier).
204 : *
205 : * \see slice_away()
206 : */
207 : // clang-tidy: incorrectly reported redundancy in template expression
208 : template <size_t N = Dim, Requires<(N > 0 and N == Dim)> = nullptr> // NOLINT
209 1 : Mesh<Dim - 1> on_interface(size_t d) const;
210 :
211 : /*!
212 : * \brief Returns a Mesh with the dimensions \p d, ... present (zero-indexed).
213 : *
214 : * \details Generally you use this method to obtain a lower-dimensional Mesh
215 : * by slicing through a subset of the dimensions. However, you can also
216 : * reorder dimensions using this method by slicing through the dimensions in
217 : * an order you choose.
218 : *
219 : * \see slice_away()
220 : */
221 : template <typename... D, Requires<(sizeof...(D) <= Dim)> = nullptr>
222 1 : Mesh<sizeof...(D)> slice_through(D... d) const {
223 : static_assert(std::conjunction_v<tt::is_integer<D>...>,
224 : "The dimensions must be integers.");
225 : const std::array<size_t, sizeof...(D)> dims{{static_cast<size_t>(d)...}};
226 : return slice_through(dims);
227 : }
228 :
229 : /*!
230 : * \brief Returns a Mesh with the dimensions \p dims present (zero-indexed).
231 : *
232 : * \see slice_through() The templated overload of this function
233 : */
234 : template <size_t SliceDim, Requires<(SliceDim <= Dim)> = nullptr>
235 1 : Mesh<SliceDim> slice_through(const std::array<size_t, SliceDim>& dims) const;
236 :
237 : /*!
238 : * \brief Returns the Meshes representing 1D slices of this Mesh.
239 : *
240 : * The is the same as the array filled with `slice_through(d)` for
241 : * `d` from `0` to `Dim - 1` except in dimension 0 where
242 : * `slice_through(d)` is not defined.
243 : */
244 1 : std::array<Mesh<1>, Dim> slices() const;
245 :
246 : // clang-tidy: runtime-references
247 : // NOLINTNEXTLINE(google-runtime-references)
248 0 : void pup(PUP::er& p);
249 :
250 : private:
251 : template <size_t LocalDim>
252 0 : friend class Mesh;
253 :
254 : template <size_t LocalDim>
255 0 : friend bool operator==(const Mesh<LocalDim>& lhs, const Mesh<LocalDim>& rhs);
256 :
257 : // - 8 bits per extent = 3 * 8 = 24 bits = 3 bytes
258 : // This limits us to 255 extent per direction. Since FD is 2N-1 DG, this
259 : // means at most 128 DG points, which is should be fine.
260 : // - We encode the basis & quadrature as two sets of 4 bits in a uint8_t.
261 0 : std::array<uint8_t, (Dim > 0 ? Dim : 1)> extents_{};
262 0 : std::array<uint8_t, (Dim > 0 ? Dim : 1)> quadrature_and_basis_{};
263 : };
264 :
265 : /*!
266 : * \ingroup DataStructuresGroup
267 : * \brief Returns `true` if the mesh is isotropic, `false` otherwise.
268 : *
269 : * If `Dim` is zero, then `true` is always returned.
270 : */
271 : template <size_t Dim>
272 1 : bool is_isotropic(const Mesh<Dim>& mesh);
273 :
274 : /// \cond HIDDEN_SYMBOLS
275 : template <size_t Dim>
276 : bool operator==(const Mesh<Dim>& lhs, const Mesh<Dim>& rhs);
277 :
278 : template <size_t Dim>
279 : bool operator!=(const Mesh<Dim>& lhs, const Mesh<Dim>& rhs);
280 :
281 : template <size_t Dim>
282 : std::ostream& operator<<(std::ostream& os, const Mesh<Dim>& mesh);
283 : /// \endcond
|
