Line data Source code
1 : #pragma once
2 :
3 : #include <iostream>
4 : #include <iterator>
5 : #include <string_view>
6 : #include <type_traits>
7 :
8 : #include "tarch/la/Vector.h"
9 :
10 : namespace coupling {
11 : namespace indexing {
12 :
13 : /**
14 : * Expresses type parametrisation of CellIndex specialisations.
15 : *
16 : * vector: implies representation as vector, false implies scalar index. \n
17 : * local: implies indexing restricted to local MD domain. \n
18 : * md2macro: implies indexing restricted to cells that are sent from MD to macro
19 : * solver. \n noGhost: implies ghost layer cells to not be included in indexing.
20 : * \n
21 : *
22 : * @author Felix Maurer
23 : */
24 : enum class IndexTrait { vector, local, md2macro, noGhost };
25 : namespace TraitOperations {
26 :
27 : /**
28 : * Returns true iff the template and runtime argument match. Curried operator==
29 : * for above enum class.
30 : */
31 : template <IndexTrait t1> constexpr bool is_same(const IndexTrait& t2) { return t1 == t2; }
32 :
33 : /**
34 : * Checks if a given parameter pack of at least two IndexTrait entries is in the
35 : * correct order. The correct order is the order in which the Traits occur in
36 : * IndexTrait's declaration, i.e.:\n
37 : *
38 : * *vector < local < md2macro < noGhost*
39 : *
40 : * This order is naturally induced by operator> on the enum class IndexTrait.
41 : *
42 : * @tparam Parameter pack of at least two IndexTraits
43 : * @return true iff the parameter pack is ordered correctly.
44 : */
45 : template <IndexTrait t1, IndexTrait t2, IndexTrait... rest> constexpr bool is_ordered() {
46 : if constexpr (sizeof...(rest) == 0) {
47 : return t1 < t2;
48 : } else {
49 : return t1 < t2 and is_ordered<t2, rest...>();
50 : }
51 : }
52 :
53 : } // namespace TraitOperations
54 :
55 : /**
56 : * Index used to describe spatial location of a CouplingCell.
57 : * Since various different ways of expressing this location are useful for
58 : * different applications, IndexTraits are used to describe the context of this
59 : * index. \n
60 : *
61 : * All commonly used (arithmetic) operations on CouplingCell indices are
62 : * provided as well as seamless conversion between any two ways of expressing
63 : * these indices. (cf. user-defined conversion function below)\n
64 : *
65 : * @tparam dim number of dimensions of the coupled simulation
66 : * @tparam traits... index type parametrisation (expressed via IndexTraits) used
67 : * by this specific index
68 : *
69 : * @author Felix Maurer
70 : */
71 : template <unsigned int dim, IndexTrait... traits> class CellIndex;
72 :
73 : /**
74 : * Base CellIndex specialisation (mainly) used for conversions.
75 : */
76 : template <unsigned int dim> using BaseIndex = CellIndex<dim, IndexTrait::vector>;
77 :
78 : // TODO: refactor as member functions
79 : template <unsigned int dim, IndexTrait... traits> unsigned int convertToScalar(const CellIndex<dim, traits...>&);
80 :
81 : template <unsigned int dim, IndexTrait... traits> tarch::la::Vector<dim, int> convertToVector(const CellIndex<dim, traits...>&);
82 : } // namespace indexing
83 : } // namespace coupling
84 :
85 33488 : template <unsigned int dim, coupling::indexing::IndexTrait... traits> class coupling::indexing::CellIndex {
86 : public:
87 : /**
88 : * Check at compile time wether traits... is in proper order, i.e.
89 : *
90 : * \forall i<j: traits[i] < traits[j]
91 : *
92 : * where '<' is the standard operator< on enum classes.
93 : * This means: vector < local < md2macro < noGhost
94 : *
95 : * Note that this causes duplicate IndexTraits in 'traits' to be not accepted.
96 : */
97 : static constexpr bool checkIndexTraitOrder() {
98 : if constexpr (sizeof...(traits) > 1) {
99 : return coupling::indexing::TraitOperations::is_ordered<traits...>();
100 : } else {
101 : return true;
102 : }
103 : }
104 : static_assert(checkIndexTraitOrder(), "Invalid order in IndexTrait parameter pack! Correct oder: "
105 : "IndexTrait::vector < IndexTrait::local < IndexTrait::md2macro "
106 : "< IndexTrait::noGhost");
107 :
108 : /**
109 : * The type of this CellIndex's underlying index representation.
110 : */
111 : using value_T = std::conditional_t<(coupling::indexing::TraitOperations::is_same<coupling::indexing::IndexTrait::vector>(traits) or ...),
112 : tarch::la::Vector<dim, int>, unsigned int>;
113 :
114 : /**
115 : * Constructors
116 : */
117 716 : CellIndex() = default;
118 2558604 : CellIndex(const CellIndex& ci) = default;
119 2678444 : explicit CellIndex(const value_T& i) : _index(i) {}
120 :
121 : /**
122 : * Conversion function: Convert to CellIndex of same dim but different
123 : * IndexType.
124 : *
125 : * @tparam convert_to_T IndexType parameter of the CellIndex specialisation to
126 : * convert to
127 : * @return CellIndex of different template parametrisation.
128 : */
129 : template <coupling::indexing::IndexTrait... converted_traits> operator CellIndex<dim, converted_traits...>() const;
130 :
131 : /**
132 : * Access to primive value_T of this index.
133 : * Should be used as sparingly as possible since it can lead to bugs
134 : * preventable by using CellIndex instances instead of primitives.
135 : *
136 : * @return unsigned integer/vector representation of this index.
137 : */
138 46390523 : value_T get() const { return _index; } // TODO: why (value_T) cast?
139 :
140 : /**
141 : * Increments the index by one.
142 : * Note that this does NOT increment indices in vector representation in all
143 : * directions.
144 : */
145 1246344 : CellIndex& operator++() {
146 : if constexpr (std::is_same_v<value_T, tarch::la::Vector<dim, int>>) {
147 659432 : ++_index[0];
148 659432 : if (_index[0] == (int)numberCellsInDomain[0]) {
149 77248 : _index[0] = 0;
150 77248 : ++_index[1];
151 77248 : if (_index[1] == (int)numberCellsInDomain[1]) {
152 16056 : _index[1] = 0;
153 : if constexpr (dim == 3)
154 16056 : ++_index[2];
155 : }
156 : }
157 : } else
158 259356 : ++_index;
159 :
160 659432 : return *this;
161 : }
162 : /**
163 : * Decrements the index by one.
164 : * Note that this does NOT decrement indices in vector representation in all
165 : * directions.
166 : */
167 433792 : CellIndex& operator--() {
168 : if constexpr (std::is_same_v<value_T, tarch::la::Vector<dim, int>>) {
169 216896 : --_index[0];
170 216896 : if (_index[0] < 0) {
171 21696 : _index[0] = numberCellsInDomain[0] - 1;
172 21696 : --_index[1];
173 21696 : if (_index[1] < 0) {
174 2432 : _index[1] = numberCellsInDomain[1] - 1;
175 : if constexpr (dim == 3)
176 2432 : --_index[2];
177 : }
178 : }
179 : } else
180 216896 : --_index;
181 :
182 216896 : return *this;
183 : }
184 :
185 : /*
186 : * Any two indices fulfill some relation iff the unsigned integers underlying
187 : * their CellIndex<dim, {}> equivalents fulfill that relation. Note that this
188 : * is different than what the static member 'contains(...)' does.
189 : *
190 : * @param CellIndex index to compare this index to
191 : */
192 2881304 : bool operator==(const CellIndex& i) const { return _index == i.get(); }
193 1164200 : bool operator!=(const CellIndex& i) const { return _index != i.get(); }
194 : bool operator<(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) < convertToScalar<dim, traits...>(i); };
195 : bool operator<=(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) <= convertToScalar<dim, traits...>(i); };
196 : bool operator>(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) > convertToScalar<dim, traits...>(i); };
197 : bool operator>=(const CellIndex& i) const { return convertToScalar<dim, traits...>(*this) >= convertToScalar<dim, traits...>(i); };
198 :
199 : /**
200 : * Initialises all static members dependant only on upperBoundary and
201 : * lowerBoundary
202 : */
203 12224 : static void setDomainParameters() {
204 110016 : auto cellnum = upperBoundary.get() - lowerBoundary.get() + tarch::la::Vector<dim, int>{1};
205 48896 : for (unsigned int d = 0; d < dim; d++)
206 37044 : cellnum[d] = std::max(0, cellnum[d]);
207 48896 : numberCellsInDomain = tarch::la::Vector<dim, unsigned int>{cellnum};
208 :
209 12224 : linearNumberCellsInDomain = 1;
210 48896 : for (unsigned int d = 0; d < dim; d++)
211 36672 : linearNumberCellsInDomain *= numberCellsInDomain[d];
212 :
213 12224 : tarch::la::Vector<dim, unsigned int> divFactor{1};
214 36672 : for (unsigned int d = 1; d < dim; d++)
215 24448 : divFactor[d] = divFactor[d - 1] * (numberCellsInDomain[d - 1]);
216 12224 : divisionFactor = divFactor;
217 12224 : }
218 :
219 : /**
220 : * Checks if a given index is within the domain between this the lower and
221 : * upper bound of this kind of indexing.
222 : *
223 : * @param CellIndex index index to be checked
224 : * @return true iff i is within the domain of this CellIndex template
225 : * specialisation.
226 : */
227 3910352 : static bool contains(const coupling::indexing::BaseIndex<dim>& index) {
228 10562976 : for (unsigned int d = 0; d < dim; d++) {
229 17323968 : if (index.get()[d] < CellIndex::lowerBoundary.get()[d])
230 : return false;
231 15314864 : if (index.get()[d] > CellIndex::upperBoundary.get()[d])
232 : return false;
233 : }
234 :
235 : return true;
236 : }
237 :
238 0 : tarch::la::Vector<dim, double> getCellMidPoint() const {
239 0 : BaseIndex<dim> globalIndex{*this};
240 0 : tarch::la::Vector<dim, double> cellMidPoint(IndexingService<dim>::getInstance().getGlobalMDDomainOffset() -
241 0 : 0.5 * IndexingService<dim>::getInstance().getCouplingCellSize());
242 0 : for (unsigned int d = 0; d < dim; d++)
243 0 : cellMidPoint[d] += ((double)(globalIndex.get()[d])) * IndexingService<dim>::getInstance().getCouplingCellSize()[d];
244 0 : return cellMidPoint;
245 : }
246 :
247 : /**
248 : * Defines where this type of cell index starts counting.
249 : * Read inclusively, e.g.: lowerBoundary = {1,2,3} means {1,2,3} is the first
250 : * index contained in this type of cell index' domain.
251 : */
252 : static BaseIndex<dim> lowerBoundary;
253 : /**
254 : * Defines where this type of cell index stops counting.
255 : * Read inclusively, e.g.: upperBoundary = {4,5,6} means {4,5,6} is the last
256 : * index contained in this type of cell index' domain.
257 : */
258 : static BaseIndex<dim> upperBoundary;
259 :
260 : /**
261 : * Number of cells in this indexing's domain. The above declared boundaries
262 : * are inclusive. Initialised in setDomainParameters().
263 : */
264 : static tarch::la::Vector<dim, unsigned int> numberCellsInDomain;
265 :
266 : /**
267 : * Number of cells in this indexing's domain. Same as numberCellsInDomain, but
268 : * expressed as a scalar.
269 : */
270 : static unsigned int linearNumberCellsInDomain;
271 :
272 : /**
273 : * Used in scalar -> vector indexing conversion functions
274 : * Initialised in setDomainParameters().
275 : */
276 : static tarch::la::Vector<dim, unsigned int> divisionFactor;
277 :
278 : class IndexIterator {
279 : public:
280 696 : IndexIterator(CellIndex x) : _idx(x) {}
281 65092 : IndexIterator(const IndexIterator& a) : _idx(a._idx) {}
282 :
283 262456 : const CellIndex& operator*() const { return _idx; }
284 : const CellIndex* operator->() const { return &_idx; }
285 :
286 : // Prefix increment
287 812552 : IndexIterator& operator++() {
288 464484 : ++_idx;
289 766140 : return *this;
290 : }
291 :
292 : // Postfix increment
293 : IndexIterator operator++(int) {
294 : IndexIterator tmp = *this;
295 : ++(*this);
296 : return tmp;
297 : }
298 :
299 : friend bool operator==(const IndexIterator& a, const IndexIterator& b) { return a._idx == b._idx; };
300 767112 : friend bool operator!=(const IndexIterator& a, const IndexIterator& b) { return a._idx != b._idx; };
301 :
302 : private:
303 : CellIndex _idx;
304 : };
305 1036 : static IndexIterator begin() {
306 : // if constexpr is neccessary here (otherwise benchmark performance breaks, due to runtime if)
307 : if constexpr (std::is_same_v<unsigned int, value_T>) {
308 : // If this CellIndex is scalar and if the domain is empty, then lowerBoundary can not be converted to CellIndex
309 : // Thus we have to return something else as an iterator for the empty domain
310 468 : if (linearNumberCellsInDomain == 0)
311 16 : return IndexIterator(CellIndex{});
312 : else
313 452 : return IndexIterator(lowerBoundary);
314 : } else
315 1096 : return IndexIterator(lowerBoundary);
316 : }
317 8988 : static IndexIterator end() {
318 : // Todo: This if-else branching is quite slow? (~factor 2)
319 8988 : if (linearNumberCellsInDomain == 0)
320 16 : return begin();
321 : else
322 26012 : return ++IndexIterator(upperBoundary);
323 : }
324 :
325 : static const char TNAME[];
326 :
327 : private:
328 : value_T _index;
329 :
330 : template <coupling::indexing::IndexTrait T1, coupling::indexing::IndexTrait... other_traits> CellIndex<dim, other_traits...> getScalarCellIndex(int value) {
331 : static_assert(T1 == coupling::indexing::IndexTrait::vector);
332 : return CellIndex<dim, other_traits...>(value);
333 : }
334 : };
335 :
336 : template <unsigned int dim, coupling::indexing::IndexTrait... lhs_traits, coupling::indexing::IndexTrait... rhs_traits>
337 1440 : bool operator==(const coupling::indexing::CellIndex<dim, lhs_traits...>& lhs, const coupling::indexing::CellIndex<dim, rhs_traits...>& rhs) {
338 : using Base = coupling::indexing::BaseIndex<dim>;
339 :
340 8640 : return Base(lhs).get() == Base(rhs).get();
341 : }
342 :
343 : template <unsigned int dim, coupling::indexing::IndexTrait... lhs_traits, coupling::indexing::IndexTrait... rhs_traits>
344 : bool operator!=(const coupling::indexing::CellIndex<dim, lhs_traits...>& lhs, const coupling::indexing::CellIndex<dim, rhs_traits...>& rhs) {
345 : return not(lhs == rhs);
346 : }
347 :
348 : // Include implementation
349 : #include "CellIndex.cpph"
|
