Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines class ElementId.
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <array>
      10             : #include <cstddef>
      11             : #include <functional>
      12             : #include <iosfwd>
      13             : #include <limits>
      14             : #include <optional>
      15             : 
      16             : #include "Domain/Structure/SegmentId.hpp"
      17             : #include "Domain/Structure/Side.hpp"
      18             : #include "Utilities/Algorithm.hpp"
      19             : #include "Utilities/ErrorHandling/Assert.hpp"
      20             : #include "Utilities/ErrorHandling/Error.hpp"
      21             : #include "Utilities/MakeArray.hpp"
      22             : 
      23             : /// \cond
      24             : namespace PUP {
      25             : class er;
      26             : }  // namespace PUP
      27             : /// \endcond
      28             : 
      29             : /*!
      30             :  * \ingroup ComputationalDomainGroup
      31             :  * \brief An ElementId uniquely labels an Element.
      32             :  *
      33             :  * It is constructed from the BlockId of the Block to which the Element belongs
      34             :  * and the VolumeDim SegmentIds that label the segments of the Block that the
      35             :  * Element covers. An optional `grid_index` identifies elements with the same
      36             :  * BlockId and SegmentIds across multiple grids.
      37             :  *
      38             :  * \details
      39             :  * The `ElementId` serves as an index that is compatible with Charm++ and
      40             :  * therefore must adhere to the restrictions imposed by Charm++. These are:
      41             :  * - `ElementId` must satisfy `std::is_pod`
      42             :  * - `ElementId` must not be larger than the size of three `int`s, i.e.
      43             :  *   `sizeof(ElementId) <= 3 * sizeof(int)`
      44             :  *
      45             :  * The latter restriction can be relaxed with a special compilation flag to
      46             :  * Charm++, but we have not yet needed more elements than can be accounted for
      47             :  * by densely packing bits together. `SegmentId` is responsible for handling the
      48             :  * low-level bit manipulations to create an index that satisfies the size
      49             :  * constraints.
      50             :  */
      51             : template <size_t VolumeDim>
      52           1 : class ElementId {
      53             :  public:
      54           0 :   static constexpr size_t block_id_bits = 24;
      55           0 :   static constexpr size_t grid_index_bits = 8;
      56           0 :   static constexpr size_t refinement_bits = 4;
      57           0 :   static constexpr size_t max_refinement_level = 16;
      58             :   // We need some padding to ensure bit fields align with type boundaries,
      59             :   // otherwise the size of `ElementId` is too large.
      60           0 :   static constexpr size_t padding = 4;
      61             :   static_assert(block_id_bits + 3 * (refinement_bits + max_refinement_level) +
      62             :                         grid_index_bits + padding ==
      63             :                     3 * 8 * sizeof(int),
      64             :                 "Bit representation requires padding or is too large");
      65             :   static_assert(two_to_the(refinement_bits) >= max_refinement_level,
      66             :                 "Not enough bits to represent all refinement levels");
      67             : 
      68           0 :   static constexpr size_t volume_dim = VolumeDim;
      69             : 
      70             :   /// Default constructor needed for Charm++ serialization.
      71           1 :   ElementId() = default;
      72           0 :   ElementId(const ElementId&) = default;
      73           0 :   ElementId& operator=(const ElementId&) = default;
      74           0 :   ElementId(ElementId&&) = default;
      75           0 :   ElementId& operator=(ElementId&&) = default;
      76           0 :   ~ElementId() = default;
      77             : 
      78             :   /// Create the ElementId of the root Element of a Block.
      79           1 :   explicit ElementId(size_t block_id, size_t grid_index = 0);
      80             : 
      81             :   /// Create an arbitrary ElementId.
      82           1 :   ElementId(size_t block_id, std::array<SegmentId, VolumeDim> segment_ids,
      83             :             size_t grid_index = 0);
      84             : 
      85             :   /// Create an ElementId from its string representation (see `operator<<`).
      86           1 :   ElementId(const std::string& grid_name);
      87             : 
      88           0 :   ElementId<VolumeDim> id_of_child(size_t dim, Side side) const;
      89             : 
      90           0 :   ElementId<VolumeDim> id_of_parent(size_t dim) const;
      91             : 
      92           0 :   size_t block_id() const { return block_id_; }
      93             : 
      94           0 :   size_t grid_index() const { return grid_index_; }
      95             : 
      96           0 :   std::array<size_t, VolumeDim> refinement_levels() const {
      97             :     if constexpr (VolumeDim == 1) {
      98             :       return {{refinement_level_xi_}};
      99             :     } else if constexpr (VolumeDim == 2) {
     100             :       return {{refinement_level_xi_, refinement_level_eta_}};
     101             :     } else if constexpr (VolumeDim == 3) {
     102             :       return {{refinement_level_xi_, refinement_level_eta_,
     103             :                refinement_level_zeta_}};
     104             :     }
     105             :   }
     106             : 
     107           0 :   std::array<SegmentId, VolumeDim> segment_ids() const {
     108             :     if constexpr (VolumeDim == 1) {
     109             :       return {{SegmentId{refinement_level_xi_, index_xi_}}};
     110             :     } else if constexpr (VolumeDim == 2) {
     111             :       return {{SegmentId{refinement_level_xi_, index_xi_},
     112             :                SegmentId{refinement_level_eta_, index_eta_}}};
     113             :     } else if constexpr (VolumeDim == 3) {
     114             :       return {{SegmentId{refinement_level_xi_, index_xi_},
     115             :                SegmentId{refinement_level_eta_, index_eta_},
     116             :                SegmentId{refinement_level_zeta_, index_zeta_}}};
     117             :     }
     118             :   }
     119             : 
     120           0 :   SegmentId segment_id(const size_t dim) const {
     121             :     ASSERT(dim < VolumeDim, "Dimension must be smaller than "
     122             :                                 << VolumeDim << ", but is: " << dim);
     123             :     switch (dim) {
     124             :       case 0:
     125             :         return {refinement_level_xi_, index_xi_};
     126             :       case 1:
     127             :         return {refinement_level_eta_, index_eta_};
     128             :       case 2:
     129             :         return {refinement_level_zeta_, index_zeta_};
     130             :       default:
     131             :         ERROR("Invalid dimension: " << dim);
     132             :     }
     133             :   }
     134             : 
     135             :   /// Returns an ElementId meant for identifying data on external boundaries,
     136             :   /// which should never correspond to the Id of an actual element.
     137           1 :   static ElementId<VolumeDim> external_boundary_id();
     138             : 
     139             :  private:
     140           0 :   uint32_t block_id_ : block_id_bits;
     141           0 :   uint32_t grid_index_ : grid_index_bits;  // end first 32 bits
     142           0 :   uint32_t index_xi_ : max_refinement_level;
     143           0 :   uint32_t index_eta_ : max_refinement_level;
     144           0 :   uint32_t index_zeta_ : max_refinement_level;
     145           0 :   uint32_t empty_ : padding;
     146           0 :   uint32_t refinement_level_xi_ : refinement_bits;  // end second 32 bits
     147           0 :   uint32_t refinement_level_eta_ : refinement_bits;
     148           0 :   uint32_t refinement_level_zeta_ : refinement_bits;  // end third 32 bits
     149             : };
     150             : 
     151             : /// \cond
     152             : // macro that generate the pup operator for SegmentId
     153             : PUPbytes(ElementId<1>)      // NOLINT
     154             :     PUPbytes(ElementId<2>)  // NOLINT
     155             :     PUPbytes(ElementId<3>)  // NOLINT
     156             :     /// \endcond
     157             : 
     158             :     /// Output operator for ElementId.
     159             :     template <size_t VolumeDim>
     160           1 :     std::ostream& operator<<(std::ostream& os, const ElementId<VolumeDim>& id);
     161             : 
     162             : /// Equivalence operator for ElementId.
     163             : template <size_t VolumeDim>
     164             : bool operator==(const ElementId<VolumeDim>& lhs,
     165             :                 const ElementId<VolumeDim>& rhs);
     166             : 
     167             : /// Inequivalence operator for ElementId.
     168             : template <size_t VolumeDim>
     169             : bool operator!=(const ElementId<VolumeDim>& lhs,
     170             :                 const ElementId<VolumeDim>& rhs);
     171             : 
     172             : /// Define an ordering of element IDs first by grid index, then by block ID,
     173             : /// then by segment ID in each dimension in turn. In each dimension, segment IDs
     174             : /// are ordered first by refinement level (which will typically be the same when
     175             : /// comparing two element IDs), and second by index. There's no particular
     176             : /// reason for this choice of ordering. For applications such as distributing
     177             : /// elements among cores, orderings such as defined by
     178             : /// `domain::BlockZCurveProcDistribution` may be more appropriate.
     179             : template <size_t VolumeDim>
     180           1 : bool operator<(const ElementId<VolumeDim>& lhs,
     181             :                const ElementId<VolumeDim>& rhs);
     182             : 
     183             : template <size_t VolumeDim>
     184           0 : bool operator>(const ElementId<VolumeDim>& lhs,
     185             :                const ElementId<VolumeDim>& rhs) {
     186             :   return rhs < lhs;
     187             : }
     188             : template <size_t VolumeDim>
     189           0 : bool operator<=(const ElementId<VolumeDim>& lhs,
     190             :                 const ElementId<VolumeDim>& rhs) {
     191             :   return !(lhs > rhs);
     192             : }
     193             : template <size_t VolumeDim>
     194           0 : bool operator>=(const ElementId<VolumeDim>& lhs,
     195             :                 const ElementId<VolumeDim>& rhs) {
     196             :   return !(lhs < rhs);
     197             : }
     198             : 
     199             : /// @{
     200             : /// \brief Returns a bool if the element is the zeroth element in the domain.
     201             : ///
     202             : /// \details An element is considered to be the zeroth element if its ElementId
     203             : /// `id` has
     204             : /// 1. id.block_id() == 0
     205             : /// 2. All id.segment_ids() have SegmentId.index() == 0
     206             : /// 3. If the argument `grid_index` is specified, id.grid_index() == grid_index.
     207             : ///
     208             : /// This means that the only element in a domain that this function will return
     209             : /// `true` for is the element in the lower corner of Block0 of that domain. The
     210             : /// `grid_index` will determine which domain is used for the comparison. During
     211             : /// evolutions, only one domain will be active at a time so it doesn't make
     212             : /// sense to compare the `grid_index`. However, during an elliptic solve
     213             : /// when there are multiple grids, this `grid_index` is useful for specifying
     214             : /// only one element over all domains.
     215             : ///
     216             : /// This function is useful if you need a unique element in the domain because
     217             : /// only one element in the whole domain can be the zeroth element.
     218             : ///
     219             : /// \parblock
     220             : /// \warning If you have multiple grids and you don't specify the `grid_index`
     221             : /// argument, this function will return `true` for one element in every grid
     222             : /// and thus can't be used to determine a unique element in a simulation; only a
     223             : /// unique element in each grid.
     224             : /// \endparblock
     225             : /// \parblock
     226             : /// \warning If the domain is re-gridded, a different ElementId may represent
     227             : /// the zeroth element.
     228             : /// \endparblock
     229             : template <size_t Dim>
     230           1 : bool is_zeroth_element(const ElementId<Dim>& id,
     231             :                        const std::optional<size_t>& grid_index);
     232             : 
     233             : // This overload is added (instead of adding a default value for grid_index)
     234             : // in order to avoid adding DomainStructures as a dependency of Parallel
     235             : // by using a forward declaration in Parallel/DistributedObject.hpp
     236             : template <size_t Dim>
     237           1 : bool is_zeroth_element(const ElementId<Dim>& id);
     238             : /// @}
     239             : // ######################################################################
     240             : // INLINE DEFINITIONS
     241             : // ######################################################################
     242             : 
     243             : template <size_t VolumeDim>
     244           0 : size_t hash_value(const ElementId<VolumeDim>& id);
     245             : 
     246             : // clang-tidy: do not modify namespace std
     247             : namespace std {  // NOLINT
     248             : template <size_t VolumeDim>
     249             : struct hash<ElementId<VolumeDim>> {
     250             :   size_t operator()(const ElementId<VolumeDim>& id) const;
     251             : };
     252             : }  // namespace std
     253             : 
     254             : template <size_t VolumeDim>
     255           1 : inline bool operator==(const ElementId<VolumeDim>& lhs,
     256             :                        const ElementId<VolumeDim>& rhs) {
     257             :   return lhs.block_id() == rhs.block_id() and
     258             :          lhs.segment_ids() == rhs.segment_ids() and
     259             :          lhs.grid_index() == rhs.grid_index();
     260             : }
     261             : 
     262             : template <size_t VolumeDim>
     263           1 : inline bool operator!=(const ElementId<VolumeDim>& lhs,
     264             :                        const ElementId<VolumeDim>& rhs) {
     265             :   return not(lhs == rhs);
     266             : }