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 <limits>
       9             : #include <ostream>
      10             : 
      11             : #include "Domain/Structure/Side.hpp"
      12             : #include "Utilities/ConstantExpressions.hpp"
      13             : #include "Utilities/Gsl.hpp"
      14             : #include "Utilities/Requires.hpp"
      15             : #include "Utilities/StdHelpers.hpp"
      16             : 
      17             : /*!
      18             :  * \brief An element of dimension `ElementDim` on the boundary of a hypercube of
      19             :  * dimension `HypercubeDim`
      20             :  *
      21             :  * A hypercube of dimension \f$n\f$ (`HypercubeDim`) is composed of
      22             :  * \f$2^{n-k}\binom{n}{k}\f$ elements of dimension \f$k \leq n\f$
      23             :  * (`ElementDim`). For example, a 3D cube has 8 vertices (\f$k=0\f$), 12 edges
      24             :  * (\f$k=1\f$), 6 faces (\f$k=2\f$) and 1 cell (\f$k=3\f$). Each element is
      25             :  * identified by the \f$k\f$ dimensions it shares with the parent hypercube
      26             :  * and \f$n - k\f$ indices that specify whether it is located on the lower or
      27             :  * upper side of the parent hypercube's remaining dimensions.
      28             :  */
      29             : template <size_t ElementDim, size_t HypercubeDim>
      30           1 : struct HypercubeElement {
      31             :   static_assert(ElementDim <= HypercubeDim);
      32             : 
      33           0 :   HypercubeElement(std::array<size_t, ElementDim> dimensions_in_parent,
      34             :                    std::array<Side, HypercubeDim - ElementDim> index);
      35             : 
      36             :   template <size_t NumIndices = HypercubeDim - ElementDim,
      37             :             Requires<NumIndices == 0> = nullptr>
      38           0 :   HypercubeElement() {
      39             :     for (size_t d = 0; d < ElementDim; ++d) {
      40             :       gsl::at(dimensions_in_parent_, d) = d;
      41             :     }
      42             :   }
      43             : 
      44             :   template <size_t LocalElementDim = ElementDim,
      45             :             Requires<LocalElementDim == 0> = nullptr>
      46           0 :   explicit HypercubeElement(std::array<Side, HypercubeDim> index);
      47             : 
      48             :   template <typename... Indices, size_t LocalElementDim = ElementDim,
      49             :             size_t LocalHypercubeDim = HypercubeDim,
      50             :             Requires<(LocalElementDim == 0 and LocalHypercubeDim > 0 and
      51             :                       sizeof...(Indices) == LocalHypercubeDim)> = nullptr>
      52           0 :   explicit HypercubeElement(Indices... indices)
      53             :       : index_{{static_cast<Side>(indices)...}} {}
      54             : 
      55             :   template <size_t LocalElementDim = ElementDim,
      56             :             Requires<LocalElementDim == 1> = nullptr>
      57           0 :   HypercubeElement(size_t dim_in_parent,
      58             :                    std::array<Side, HypercubeDim - 1> index);
      59             : 
      60             :   /// @{
      61             :   /// The parent hypercube's dimensions that this element shares
      62           1 :   const std::array<size_t, ElementDim>& dimensions_in_parent() const;
      63             : 
      64             :   template <size_t LocalElementDim = ElementDim,
      65             :             Requires<LocalElementDim == 1> = nullptr>
      66           1 :   size_t dimension_in_parent() const;
      67             :   /// @}
      68             : 
      69             :   /// @{
      70             :   /// Whether this element is located on the lower or upper side in those
      71             :   /// dimensions that it does not share with its parent hypercube
      72           1 :   const std::array<Side, HypercubeDim - ElementDim>& index() const;
      73             : 
      74           1 :   const Side& side_in_parent_dimension(size_t d) const;
      75             : 
      76             :   template <size_t NumIndices = HypercubeDim - ElementDim,
      77             :             Requires<NumIndices == 1> = nullptr>
      78           1 :   const Side& side() const {
      79             :     return index_[0];
      80             :   }
      81             :   /// @}
      82             : 
      83           0 :   bool operator==(const HypercubeElement& rhs) const {
      84             :     return dimensions_in_parent_ == rhs.dimensions_in_parent_ and
      85             :            index_ == rhs.index_;
      86             :   }
      87             : 
      88           0 :   bool operator!=(const HypercubeElement& rhs) const {
      89             :     return not(*this == rhs);
      90             :   }
      91             : 
      92             :  private:
      93           0 :   std::array<size_t, ElementDim> dimensions_in_parent_{};
      94           0 :   std::array<Side, HypercubeDim - ElementDim> index_{};
      95             : };
      96             : 
      97             : template <size_t ElementDim, size_t HypercubeDim>
      98           0 : std::ostream& operator<<(
      99             :     std::ostream& os,
     100             :     const HypercubeElement<ElementDim, HypercubeDim>& element);
     101             : 
     102             : /// A vertex in a `Dim`-dimensional hypercube
     103             : template <size_t Dim>
     104           1 : using Vertex = HypercubeElement<0, Dim>;
     105             : 
     106             : /// An edge in a `Dim`-dimensional hypercube
     107             : template <size_t Dim>
     108           1 : using Edge = HypercubeElement<1, Dim>;
     109             : 
     110             : /// A face in a `Dim`-dimensional hypercube
     111             : template <size_t Dim>
     112           1 : using Face = HypercubeElement<2, Dim>;
     113             : 
     114             : /// A cell in a `Dim`-dimensional hypercube
     115             : template <size_t Dim>
     116           1 : using Cell = HypercubeElement<3, Dim>;
     117             : 
     118             : /*!
     119             :  * \brief Iterator over all `ElementDim`-dimensional elements on the boundary of
     120             :  * a `HypercubeDim`-dimensional hypercube.
     121             :  *
     122             :  * \see `HypercubeElement`
     123             :  */
     124             : template <size_t ElementDim, size_t HypercubeDim>
     125           1 : struct HypercubeElementsIterator {
     126             :   static_assert(
     127             :       ElementDim <= HypercubeDim,
     128             :       "Hypercube element dimension must not exceed hypercube dimension.");
     129             : 
     130             :  public:
     131           0 :   static constexpr size_t num_indices = HypercubeDim - ElementDim;
     132             : 
     133             :   /// The number of `ElementDim`-dimensional elements on the boundary of a
     134             :   /// `HypercubeDim`-dimensional hypercube.
     135           1 :   static constexpr size_t size() {
     136             :     return two_to_the(num_indices) * factorial(HypercubeDim) /
     137             :            factorial(ElementDim) / factorial(num_indices);
     138             :   }
     139             : 
     140           0 :   HypercubeElementsIterator();
     141             : 
     142           0 :   static HypercubeElementsIterator begin();
     143             : 
     144           0 :   static HypercubeElementsIterator end();
     145             : 
     146           0 :   HypercubeElementsIterator& operator++();
     147             : 
     148           0 :   HypercubeElementsIterator operator++(int);
     149             : 
     150           0 :   HypercubeElement<ElementDim, HypercubeDim> operator*() const;
     151             : 
     152             :  private:
     153             :   template <size_t LocalElementDim = ElementDim,
     154             :             Requires<(LocalElementDim > 0)> = nullptr>
     155           0 :   void increment_dimension_in_parent(size_t d);
     156             : 
     157             :   template <size_t LocalElementDim, size_t LocalHypercubeDim>
     158             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     159           0 :   friend bool operator==(
     160             :       const HypercubeElementsIterator<LocalElementDim, LocalHypercubeDim>& lhs,
     161             :       const HypercubeElementsIterator<LocalElementDim, LocalHypercubeDim>& rhs);
     162             : 
     163           0 :   std::array<size_t, ElementDim> dimensions_in_parent_{};
     164           0 :   size_t index_ = std::numeric_limits<size_t>::max();
     165             : };
     166             : 
     167             : template <size_t ElementDim, size_t HypercubeDim>
     168           0 : bool operator!=(const HypercubeElementsIterator<ElementDim, HypercubeDim>& lhs,
     169             :                 const HypercubeElementsIterator<ElementDim, HypercubeDim>& rhs);
     170             : 
     171             : /// Iterate over all vertices in a `Dim`-dimensional hypercube
     172             : template <size_t Dim>
     173           1 : using VertexIterator = HypercubeElementsIterator<0, Dim>;
     174             : 
     175             : /// Iterate over all edges in a `Dim`-dimensional hypercube
     176             : template <size_t Dim>
     177           1 : using EdgeIterator = HypercubeElementsIterator<1, Dim>;
     178             : 
     179             : /// Iterate over all faces in a `Dim`-dimensional hypercube
     180             : template <size_t Dim>
     181           1 : using FaceIterator = HypercubeElementsIterator<2, Dim>;
     182             : 
     183             : /// Iterate over all cells in a `Dim`-dimensional hypercube
     184             : template <size_t Dim>
     185           1 : using CellIterator = HypercubeElementsIterator<3, Dim>;