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