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 : }