Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines class SegmentId.
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <cstddef>
      10             : #include <iosfwd>
      11             : #include <pup.h>
      12             : 
      13             : #include "Domain/Structure/Side.hpp"
      14             : #include "Utilities/ConstantExpressions.hpp"
      15             : #include "Utilities/ErrorHandling/Assert.hpp"
      16             : 
      17             : /*!
      18             :  * \ingroup ComputationalDomainGroup
      19             :  * \brief A SegmentId labels a segment of the interval \f$[-1,1]\f$ and is used
      20             :  * to identify the bounds of an Element in a Block in each dimension.
      21             :  *
      22             :  * In \f$d\f$ dimensions, \f$d\f$ SegmentId%s are used to identify an Element.
      23             :  * In each dimension, a segment spans the subinterval \f$[-1 + 2 \frac{i}{N}, -1
      24             :  * + 2 \frac{i+1}{N}]\f$ of the logical coordinates of a Block, where \f$i \f$=
      25             :  * `index` and \f$N = 2^L\f$ where \f$L\f$ = `refinement_level`.
      26             :  *
      27             :  *  \image html SegmentId.svg  "SegmentIds"
      28             :  *
      29             :  * In the figure, The `index` of segments increase from the `lower side` to the
      30             :  * `upper side` in each dimension of a Block, while the `refinement level`
      31             :  * increases as the segments are subdivided.  For example, let the segment
      32             :  * labeled `self` be on `refinement level` \f$L\f$, with `index` \f$i\f$.  Its
      33             :  * `parent` segment is on `refinement level` \f$L-1\f$ with  `index`
      34             :  * \f$\frac{i-1}{2}\f$. The `children` of `self` are on `refinement level`
      35             :  * \f$L+1\f$, and have `index` \f$2i\f$ and \f$2i+1\f$ for the lower and upper
      36             :  * child respectively.  Also labeled on the figure are the `sibling` and
      37             :  * `abutting nibling` (child of sibling) of `self`. These relationships between
      38             :  * segments are important for h-refinement, since in each dimension an Element
      39             :  * can be flagged to split into its two `children` segments, or join with its
      40             :  * `sibling` segment to form its `parent` segment.  As refinement levels of
      41             :  * neighboring elements are kept within one, in the direction of its `sibling`,
      42             :  * a segment can only abut its `sibling` or `abutting nibling`, while on the
      43             :  * opposite side, it can abut a segment on its level, the next-lower, or the
      44             :  * next-higher level.
      45             :  */
      46           1 : class SegmentId {
      47             :  public:
      48             :   /// Default constructor needed for Charm++ serialization.
      49           1 :   SegmentId() = default;
      50           0 :   SegmentId(const SegmentId& segment_id) = default;
      51           0 :   SegmentId(SegmentId&& segment_id) = default;
      52           0 :   ~SegmentId() = default;
      53           0 :   SegmentId& operator=(const SegmentId& segment_id) = default;
      54           0 :   SegmentId& operator=(SegmentId&& segment_id) = default;
      55             : 
      56           0 :   SegmentId(size_t refinement_level, size_t index);
      57             : 
      58           0 :   constexpr size_t refinement_level() const { return refinement_level_; }
      59             : 
      60           0 :   constexpr size_t index() const { return index_; }
      61             : 
      62           0 :   SegmentId id_of_parent() const;
      63             : 
      64           0 :   SegmentId id_of_child(Side side) const;
      65             : 
      66             :   /// The other child of the parent of this segment
      67           1 :   SegmentId id_of_sibling() const;
      68             : 
      69             :   /// The child of the sibling of this segment that shares an endpoint with it
      70           1 :   SegmentId id_of_abutting_nibling() const;
      71             : 
      72             :   /// The side on which this segment shares an endpoint with its sibling
      73           1 :   Side side_of_sibling() const;
      74             : 
      75             :   /// The id this segment would have if the coordinate axis were flipped.
      76           1 :   SegmentId id_if_flipped() const;
      77             : 
      78             :   /// The block logical coordinate of the endpoint of the segment on the given
      79             :   /// Side.
      80           1 :   double endpoint(Side side) const;
      81             : 
      82             :   /// The block logical coordinate of the midpoint of the segment
      83           1 :   double midpoint() const {
      84             :     return -1.0 + (1.0 + 2.0 * index_) / two_to_the(refinement_level_);
      85             :   }
      86             : 
      87             :   /// Does the segment overlap with another?
      88           1 :   bool overlaps(const SegmentId& other) const;
      89             : 
      90             :   // NOLINTNEXTLINE
      91           0 :   void pup(PUP::er& p);
      92             : 
      93             :  private:
      94           0 :   static constexpr size_t max_refinement_level = 16;
      95           0 :   size_t refinement_level_;
      96           0 :   size_t index_;
      97             : };
      98             : 
      99             : /// Output operator for SegmentId.
     100           1 : std::ostream& operator<<(std::ostream& os, const SegmentId& id);
     101             : 
     102             : /// Equivalence operator for SegmentId.
     103             : bool operator==(const SegmentId& lhs, const SegmentId& rhs);
     104             : 
     105             : /// Inequivalence operator for SegmentId.
     106             : bool operator!=(const SegmentId& lhs, const SegmentId& rhs);
     107             : 
     108             : //##############################################################################
     109             : // INLINE DEFINITIONS
     110             : //##############################################################################
     111             : 
     112             : inline SegmentId SegmentId::id_of_parent() const {
     113             :   ASSERT(0 != refinement_level_,
     114             :          "Cannot call id_of_parent() on root refinement level!");
     115             :   // The parent has half as many segments as the child.
     116             :   return {refinement_level() - 1, index() / 2};
     117             : }
     118             : 
     119             : inline SegmentId SegmentId::id_of_child(Side side) const {
     120             :   // We cannot ASSERT on the maximum level because it's only known at runtime
     121             :   // and only known elsewhere in the code, not by SegmentId. I.e. SegmentId is
     122             :   // too low-level to know about this.
     123             :   // The child has twice as many segments as the parent, so for a particular
     124             :   // parent segment, there is both an upper and lower child segment.
     125             :   if (Side::Lower == side) {
     126             :     return {refinement_level() + 1, index() * 2};
     127             :   }
     128             :   return {refinement_level() + 1, 1 + index() * 2};
     129             : }
     130             : 
     131             : inline SegmentId SegmentId::id_of_sibling() const {
     132             :   ASSERT(0 != refinement_level(),
     133             :          "The segment on the root refinement level has no sibling");
     134             :   return {refinement_level(), (0 == index() % 2 ? index() + 1 : index() - 1)};
     135             : }
     136             : 
     137             : inline SegmentId SegmentId::id_of_abutting_nibling() const {
     138             :   ASSERT(0 != refinement_level(),
     139             :          "The segment on the root refinement level has no abutting nibling");
     140             :   return {refinement_level() + 1,
     141             :           (0 == index() % 2 ? 2 * index() + 2 : 2 * index() - 1)};
     142             : }
     143             : 
     144             : inline Side SegmentId::side_of_sibling() const {
     145             :   ASSERT(0 != refinement_level(),
     146             :          "The segment on the root refinement level has no sibling");
     147             :   return 0 == index() % 2 ? Side::Upper : Side::Lower;
     148             : }
     149             : 
     150             : inline SegmentId SegmentId::id_if_flipped() const {
     151             :   return {refinement_level(), two_to_the(refinement_level()) - 1 - index()};
     152             : }
     153             : 
     154             : inline double SegmentId::endpoint(Side side) const {
     155             :   if (Side::Lower == side) {
     156             :     return -1.0 + (2.0 * index()) / two_to_the(refinement_level());
     157             :   }
     158             :   return -1.0 + (2.0 * index() + 2.0) / two_to_the(refinement_level());
     159             : }
     160             : 
     161             : inline bool SegmentId::overlaps(const SegmentId& other) const {
     162             :   const size_t this_denom = two_to_the(refinement_level());
     163             :   const size_t other_denom = two_to_the(other.refinement_level());
     164             :   return index() * other_denom < (other.index() + 1) * this_denom and
     165             :          other.index() * this_denom < (index() + 1) * other_denom;
     166             : }
     167             : 
     168             : // These are defined so that a SegmentId can be used as part of a key of an
     169             : // unordered_set or unordered_map.
     170             : 
     171             : // hash_value is called by boost::hash and related functions.
     172           0 : size_t hash_value(const SegmentId& s);
     173             : 
     174             : namespace std {
     175             : template <>
     176             : struct hash<SegmentId> {
     177             :   size_t operator()(const SegmentId& segment_id) const;
     178             : };
     179             : }  // namespace std
     180             : 
     181           1 : inline bool operator==(const SegmentId& lhs, const SegmentId& rhs) {
     182             :   return (lhs.refinement_level() == rhs.refinement_level() and
     183             :           lhs.index() == rhs.index());
     184             : }
     185             : 
     186           1 : inline bool operator!=(const SegmentId& lhs, const SegmentId& rhs) {
     187             :   return not(lhs == rhs);
     188             : }