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"  // IWYU pragma: keep
      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             :  * \tparam Dim the number of dimensions of the computational grid.
      48             :  */
      49             : template <size_t Dim>
      50           1 : class Mesh {
      51             :  public:
      52           0 :   static constexpr size_t dim = Dim;
      53             : 
      54           0 :   struct Extents {
      55           0 :     using type = size_t;
      56           0 :     static constexpr Options::String help = {
      57             :         "The number of collocation points per dimension"};
      58             :   };
      59             : 
      60           0 :   struct Basis {
      61           0 :     using type = Spectral::Basis;
      62           0 :     static constexpr Options::String help = {
      63             :         "The choice of spectral basis to compute the collocation points"};
      64             :   };
      65             : 
      66           0 :   struct Quadrature {
      67           0 :     using type = Spectral::Quadrature;
      68           0 :     static constexpr Options::String help = {
      69             :         "The choice of quadrature to compute the collocation points"};
      70             :   };
      71             : 
      72           0 :   using options = tmpl::list<Extents, Basis, Quadrature>;
      73             : 
      74           0 :   static constexpr Options::String help =
      75             :       "Holds the number of grid points, basis, and quadrature in each "
      76             :       "direction of the computational grid. "
      77             :       "A mesh encapsulates all information necessary to construct the "
      78             :       "placement of grid points in the computational domain. It does so "
      79             :       "through a choice of basis functions, quadrature and number of points "
      80             :       "in each dimension.";
      81             : 
      82           0 :   Mesh() = default;
      83             : 
      84             :   /*!
      85             :    * \brief Construct a computational grid with the same number of collocation
      86             :    * points in each dimension.
      87             :    *
      88             :    * \param isotropic_extents The number of collocation points in each
      89             :    * dimension.
      90             :    * \param basis The choice of spectral basis to compute the
      91             :    * collocation points
      92             :    * \param quadrature The choice of quadrature to compute
      93             :    * the collocation points
      94             :    *
      95             :    * \note Because a `Mesh<0>` extends over no dimensions, it has 1 grid point
      96             :    * independent of the value of `isotropic_extents`.
      97             :    */
      98           1 :   Mesh(const size_t isotropic_extents, const Spectral::Basis basis,
      99             :        const Spectral::Quadrature quadrature)
     100             :       : extents_(isotropic_extents) {
     101             :     ASSERT(basis != Spectral::Basis::SphericalHarmonic,
     102             :            "SphericalHarmonic is not a valid basis for the Mesh");
     103             :     bases_.fill(basis);
     104             :     quadratures_.fill(quadrature);
     105             :   }
     106             : 
     107             :   /*!
     108             :    * \brief Construct a computational grid where each dimension can have a
     109             :    * different number of collocation points.
     110             :    *
     111             :    * \param extents The number of collocation points per dimension
     112             :    * \param basis The choice of spectral basis to compute the
     113             :    * collocation points
     114             :    * \param quadrature The choice of quadrature to compute
     115             :    * the collocation points
     116             :    */
     117           1 :   Mesh(std::array<size_t, Dim> extents, const Spectral::Basis basis,
     118             :        const Spectral::Quadrature quadrature)
     119             :       : extents_(std::move(extents)) {
     120             :     ASSERT(basis != Spectral::Basis::SphericalHarmonic,
     121             :            "SphericalHarmonic is not a valid basis for the Mesh");
     122             :     bases_.fill(basis);
     123             :     quadratures_.fill(quadrature);
     124             :   }
     125             : 
     126             :   /*!
     127             :    * \brief Construct a computational grid where each dimension can have both a
     128             :    * different number and placement of collocation points.
     129             :    *
     130             :    * \param extents The number of collocation points per dimension
     131             :    * \param bases The choice of spectral bases to compute the
     132             :    * collocation points per dimension
     133             :    * \param quadratures The choice of quadratures to compute
     134             :    * the collocation points per dimension
     135             :    */
     136           1 :   Mesh(std::array<size_t, Dim> extents, std::array<Spectral::Basis, Dim> bases,
     137             :        std::array<Spectral::Quadrature, Dim> quadratures)
     138             :       : extents_(std::move(extents)), quadratures_(std::move(quadratures)) {
     139             :     for (auto it = bases.begin(); it != bases.end(); it++) {
     140             :       ASSERT(*it != Spectral::Basis::SphericalHarmonic,
     141             :              "SphericalHarmonic is not a valid basis for the Mesh");
     142             :     }
     143             :     bases_ = std::move(bases);
     144             :   }
     145             : 
     146             :   /*!
     147             :    * \brief The number of grid points in each dimension of the grid.
     148             :    */
     149           1 :   const Index<Dim>& extents() const { return extents_; }
     150             : 
     151             :   /*!
     152             :    * \brief The number of grid points in dimension \p d of the grid
     153             :    * (zero-indexed).
     154             :    */
     155           1 :   size_t extents(const size_t d) const { return extents_[d]; }
     156             : 
     157             :   /*!
     158             :    * \brief The total number of grid points in all dimensions.
     159             :    *
     160             :    * \details `DataVector`s that represent field values on the grid have this
     161             :    * many entries.
     162             :    *
     163             :    * \note A zero-dimensional mesh has one grid point, since it is the slice
     164             :    * through a one-dimensional mesh (a line).
     165             :    */
     166           1 :   size_t number_of_grid_points() const { return extents_.product(); }
     167             : 
     168             :   /*!
     169             :    * \brief Returns the 1-dimensional index corresponding to the `Dim`
     170             :    * dimensional `index`.
     171             :    *
     172             :    * The first dimension varies fastest.
     173             :    *
     174             :    * \see collapsed_index()
     175             :    */
     176           1 :   size_t storage_index(const Index<Dim>& index) const {
     177             :     return collapsed_index(index, extents_);
     178             :   }
     179             : 
     180             :   /*!
     181             :    * \brief The basis chosen in each dimension of the grid.
     182             :    */
     183           1 :   const std::array<Spectral::Basis, Dim>& basis() const { return bases_; }
     184             : 
     185             :   /*!
     186             :    * \brief The basis chosen in dimension \p d of the grid (zero-indexed).
     187             :    */
     188           1 :   Spectral::Basis basis(const size_t d) const { return gsl::at(bases_, d); }
     189             : 
     190             :   /*!
     191             :    * \brief The quadrature chosen in each dimension of the grid.
     192             :    */
     193           1 :   const std::array<Spectral::Quadrature, Dim>& quadrature() const {
     194             :     return quadratures_;
     195             :   }
     196             : 
     197             :   /*!
     198             :    * \brief The quadrature chosen in dimension \p d of the grid (zero-indexed).
     199             :    */
     200           1 :   Spectral::Quadrature quadrature(const size_t d) const {
     201             :     return gsl::at(quadratures_, d);
     202             :   }
     203             : 
     204             :   /*!
     205             :    * \brief Returns a Mesh with dimension \p d removed (zero-indexed).
     206             :    *
     207             :    * \see slice_through()
     208             :    */
     209             :   // clang-tidy: incorrectly reported redundancy in template expression
     210             :   template <size_t N = Dim, Requires<(N > 0 and N == Dim)> = nullptr>  // NOLINT
     211           1 :   Mesh<Dim - 1> slice_away(size_t d) const;
     212             : 
     213             :   /*!
     214             :    * \brief Returns a Mesh with the dimensions \p d, ... present (zero-indexed).
     215             :    *
     216             :    * \details Generally you use this method to obtain a lower-dimensional Mesh
     217             :    * by slicing through a subset of the dimensions. However, you can also
     218             :    * reorder dimensions using this method by slicing through the dimensions in
     219             :    * an order you choose.
     220             :    *
     221             :    * \see slice_away()
     222             :    */
     223             :   template <typename... D, Requires<(sizeof...(D) <= Dim)> = nullptr>
     224           1 :   Mesh<sizeof...(D)> slice_through(D... d) const {
     225             :     static_assert(std::conjunction_v<tt::is_integer<D>...>,
     226             :                   "The dimensions must be integers.");
     227             :     const std::array<size_t, sizeof...(D)> dims{{static_cast<size_t>(d)...}};
     228             :     return slice_through(dims);
     229             :   }
     230             : 
     231             :   /*!
     232             :    * \brief Returns a Mesh with the dimensions \p dims present (zero-indexed).
     233             :    *
     234             :    * \see slice_through() The templated overload of this function
     235             :    */
     236             :   template <size_t SliceDim, Requires<(SliceDim <= Dim)> = nullptr>
     237           1 :   Mesh<SliceDim> slice_through(const std::array<size_t, SliceDim>& dims) const;
     238             : 
     239             :   /*!
     240             :    * \brief Returns the Meshes representing 1D slices of this Mesh.
     241             :    *
     242             :    * The is the same as the array filled with `slice_through(d)` for
     243             :    * `d` from `0` to `Dim - 1` except in dimension 0 where
     244             :    * `slice_through(d)` is not defined.
     245             :    */
     246           1 :   std::array<Mesh<1>, Dim> slices() const;
     247             : 
     248             :   // clang-tidy: runtime-references
     249             :   // NOLINTNEXTLINE(google-runtime-references)
     250           0 :   void pup(PUP::er& p);
     251             : 
     252             :  private:
     253           0 :   Index<Dim> extents_{};
     254           0 :   std::array<Spectral::Basis, Dim> bases_{};
     255           0 :   std::array<Spectral::Quadrature, Dim> quadratures_{};
     256             : };
     257             : 
     258             : /*!
     259             :  * \ingroup DataStructuresGroup
     260             :  * \brief Returns `true` if the mesh is isotropic, `false` otherwise.
     261             :  *
     262             :  * If `Dim` is zero, then `true` is always returned.
     263             :  */
     264             : template <size_t Dim>
     265           1 : bool is_isotropic(const Mesh<Dim>& mesh);
     266             : 
     267             : /// \cond HIDDEN_SYMBOLS
     268             : template <size_t Dim>
     269             : bool operator==(const Mesh<Dim>& lhs, const Mesh<Dim>& rhs);
     270             : 
     271             : template <size_t Dim>
     272             : bool operator!=(const Mesh<Dim>& lhs, const Mesh<Dim>& rhs);
     273             : 
     274             : template <size_t Dim>
     275             : std::ostream& operator<<(std::ostream& os, const Mesh<Dim>& mesh);
     276             : /// \endcond