1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=4 sw=2 sts=2: 3 #ifndef DUNE_GRID_YASPGRID_YGRID_HH 4 #define DUNE_GRID_YASPGRID_YGRID_HH 5 6 #include <array> 7 #include <vector> 8 #include <bitset> 9 #include <deque> 10 11 #include <dune/common/fvector.hh> 12 #include <dune/common/power.hh> 13 #include <dune/common/streamoperators.hh> 14 15 /** \file 16 \brief This provides a YGrid, the elemental component of the yaspgrid implementation 17 */ 18 19 namespace Dune { 20 21 namespace Yasp { 22 /** @returns an array containing the sizes of the grids associated with vectors in given array. 23 * Needed in this form due to the need of such functionality in class initializer lists. 24 * @param v the array of vectors to examine 25 */ 26 template<int d, typename ct> sizeArray(const std::array<std::vector<ct>,d> & v)27 std::array<int,d> sizeArray(const std::array<std::vector<ct>,d>& v) 28 { 29 std::array<int,d> tmp; 30 for (int i=0; i<d; ++i) 31 tmp[i] = v[i].size() - 1; 32 return tmp; 33 } 34 } //namespace Yasp 35 36 /** 37 The YGrid considered here describes a finite set \f$d\f$-tupels of the form 38 \f[ G = \{ (k_0,\ldots,k_{d-1}) | o_i \leq k_i < o_i+s_i \} \f] 39 40 together with an affine mapping. 41 42 A YGrid is characterized by the following quantities: 43 44 - The origin \f$ o=(o_0,\ldots,o_{d-1}) \in Z^d\f$, 45 - the size \f$ s=(s_0,\ldots,s_{d-1}) \in Z^d\f$, 46 - The shift \f$ r=(r_0,\ldots,r_{d-1}) \in R^d\f$. 47 - a coordinate container, that gives the mapping of the index to global coordinates (see coordinates.hh) 48 49 The shift can be used to interpret the points of a grid as midpoints of cells, faces, edges, etc. 50 51 Here is a graphical illustration of a grid: 52 53 \image html grid.png "A YGrid." 54 \image latex grid.eps "A YGrid." width=\textwidth 55 56 A YGrid allows to iterate over all its cells with an Iterator class. 57 58 A YGrid is always considered as being embedded in a larger grid. 59 This embedding is characterized by an offset and an enclosing grid as 60 shown in the following picture: 61 62 \image html subgrid.png "The SubYGrid is shown in red, blue is the enclosing grid." 63 \image latex subgrid.eps "The SubYGrid is shown in red, blue is the enclosing grid." width=\textwidth 64 65 The iterator provides also a mapping to the consecutive index in the enclosing grid. 66 67 Note: as of november 2013 there are only YGrid and YGrid::Iterator. These represent 68 the functionality of former SubYGrid and SubYGrid::TransformingSubIterator. All other 69 classes in the hierarchy have not been used. 70 */ 71 template<class Coordinates> 72 class YGridComponent 73 { 74 public: 75 //extract coordinate type and dimension from the coordinate container 76 typedef typename Coordinates::ctype ct; 77 static const int d = Coordinates::dimension; 78 79 typedef std::array<int, d> iTupel; 80 typedef FieldVector<ct,d> fTupel; 81 82 //! make uninitialized ygrid YGridComponent()83 YGridComponent () : _shift(0ULL) 84 { 85 std::fill(_origin.begin(), _origin.end(), 0); 86 std::fill(_offset.begin(), _offset.end(), 0); 87 std::fill(_size.begin(), _size.end(), 0); 88 } 89 90 /** @brief make ygrid without coordinate information 91 * @param origin origin of the grid in global coordinates 92 * @param size size of the grid 93 * Such grid has no coordinate information stored but can be 94 * used to determine an intersection with a grid with coordinate 95 * information. This avoids sending coordinates in the parallel case. 96 */ YGridComponent(iTupel origin,iTupel size)97 YGridComponent(iTupel origin, iTupel size) 98 : _origin(origin), _size(size) 99 {} 100 101 /** @brief make a subgrid by taking coordinates from a larger grid 102 * @param origin origin of the grid to be constructed 103 * @param size size of the grid to be constructed 104 * @param enclosing the grid to take coordinates and shift vector from 105 */ YGridComponent(iTupel origin,iTupel size,const YGridComponent<Coordinates> & enclosing)106 YGridComponent (iTupel origin, iTupel size, const YGridComponent<Coordinates>& enclosing) 107 : _origin(origin), _shift(enclosing.shift()), _coords(enclosing.getCoords()), _size(size), _supersize(enclosing.supersize()) 108 { 109 for (int i=0; i<d; i++) 110 _offset[i] = origin[i] - enclosing.origin(i) + enclosing.offset(i); 111 112 // compute superincrements 113 int inc = 1; 114 for (int i=0; i<d; ++i) 115 { 116 _superincrement[i] = inc; 117 inc *= _supersize[i]; 118 } 119 } 120 121 /** @brief Make YGridComponent by giving all parameters 122 * @param origin the origin of the grid in global coordinates 123 * @param shift the shift vector 124 * @param coords the coordinate vectors to be used 125 * @param size the size vector 126 * @param offset the offset in the enclosing grid 127 * @param supersize size of the enclosing grid 128 */ YGridComponent(iTupel origin,std::bitset<d> shift,Coordinates * coords,iTupel size,iTupel offset,iTupel supersize)129 YGridComponent (iTupel origin, std::bitset<d> shift, Coordinates* coords, iTupel size, iTupel offset, iTupel supersize) 130 : _origin(origin), _shift(shift), _coords(coords), _size(size), _offset(offset), _supersize(supersize) 131 { 132 // compute superincrements 133 int inc = 1; 134 for (int i=0; i<d; ++i) 135 { 136 _superincrement[i] = inc; 137 inc *= _supersize[i]; 138 } 139 } 140 141 //! Return origin in direction i origin(int i) const142 int origin (int i) const 143 { 144 return _origin[i]; 145 } 146 147 //! return reference to origin origin() const148 const iTupel& origin () const 149 { 150 return _origin; 151 } 152 153 //! Return shift in direction i shift(int i) const154 bool shift (int i) const 155 { 156 return _shift[i]; 157 } 158 159 //! Return shift tupel shift() const160 const std::bitset<d>& shift () const 161 { 162 return _shift; 163 } 164 getCoords() const165 Coordinates* getCoords() const 166 { 167 return _coords; 168 } 169 170 //! Return offset to origin of enclosing grid offset(int i) const171 int offset (int i) const 172 { 173 return _offset[i]; 174 } 175 176 //! Return offset to origin of enclosing grid offset() const177 const iTupel & offset () const 178 { 179 return _offset; 180 } 181 182 //! return size of enclosing grid supersize(int i) const183 int supersize (int i) const 184 { 185 return _supersize[i]; 186 } 187 188 //! return size of enclosing grid supersize() const189 const iTupel & supersize () const 190 { 191 return _supersize; 192 } 193 194 //! return size in direction i size(int i) const195 int size (int i) const 196 { 197 return _size[i]; 198 } 199 200 //! retrun size size() const201 iTupel size () const 202 { 203 return _size; 204 } 205 206 //! Return total size of index set which is the product of all size per direction. totalsize() const207 int totalsize () const 208 { 209 int s=1; 210 for (int i=0; i<d; ++i) 211 s *= size(i); 212 return s; 213 } 214 215 //! Return minimum index in direction i min(int i) const216 int min (int i) const 217 { 218 return _origin[i]; 219 } 220 221 //! Return maximum index in direction i max(int i) const222 int max (int i) const 223 { 224 return _origin[i] + size(i) - 1; 225 } 226 227 //! Return true if YGrid is empty, i.e. has size 0 in all directions. empty() const228 bool empty () const 229 { 230 for (int i=0; i<d; ++i) 231 { 232 if (size(i) == 0) 233 return true; 234 } 235 return false; 236 } 237 238 //! given a coordinate, return true if it is in the grid inside(const iTupel & coord) const239 bool inside (const iTupel& coord) const 240 { 241 for (int i=0; i<d; i++) 242 { 243 if ((coord[i]<_origin[i]) || (coord[i]>=_origin[i]+_size[i])) 244 return false; 245 } 246 return true; 247 } 248 249 //! given a tupel compute its index in the lexicographic numbering index(const iTupel & coord) const250 int index (const iTupel& coord) const 251 { 252 int index = (coord[d-1]-_origin[d-1]); 253 254 for (int i=d-2; i>=0; i--) 255 index = index*_size[i] + (coord[i]-_origin[i]); 256 257 return index; 258 } 259 260 //! return grid moved by the vector v move(iTupel v) const261 YGridComponent<Coordinates> move (iTupel v) const 262 { 263 for (int i=0; i<d; i++) 264 v[i] += _origin[i]; 265 return YGridComponent<Coordinates>(v,_size,*this); 266 } 267 268 //! Return YGridComponent of supergrid of self which is the intersection of self and another YGridComponent intersection(const YGridComponent<Coordinates> & r) const269 YGridComponent<Coordinates> intersection (const YGridComponent<Coordinates>& r) const 270 { 271 for (int i=0; i<d; i++) 272 { 273 //empty coordinate vectors result in empty intersections 274 if (empty() || r.empty()) 275 return YGridComponent<Coordinates>(); 276 } 277 278 iTupel neworigin; 279 iTupel newsize; 280 for (int i=0; i<d; ++i) 281 { 282 neworigin[i] = std::max(origin(i),r.origin(i)); 283 newsize[i] = std::min(max(i),r.max(i)) - neworigin[i] + 1; 284 } 285 286 return YGridComponent<Coordinates>(neworigin,newsize,*this); 287 } 288 289 290 /** Iterator class allows one to run over all cells of a grid. 291 * The cells of the grid to iterate over are numbered consecutively starting 292 * with zero. Via the index() method the iterator provides a mapping of the 293 * cells of the grid to a one-dimensional array. The number of entries 294 * in this array must be the size of the grid. 295 */ 296 class Iterator { 297 public: 298 // default constructor 299 Iterator () = default; 300 301 //! Make iterator pointing to first cell in a grid. Iterator(const YGridComponent<Coordinates> & r)302 Iterator (const YGridComponent<Coordinates>& r) : _grid(&r) 303 { 304 iTupel coord(r.origin()); 305 reinit(r,coord); 306 } 307 308 //! Make iterator pointing to given cell in a grid. Iterator(const YGridComponent<Coordinates> & r,const iTupel & coord)309 Iterator (const YGridComponent<Coordinates>& r, const iTupel& coord) 310 { 311 reinit(r,coord); 312 } 313 314 //! reinitialize iterator to given position reinit(const YGridComponent<Coordinates> & r,const iTupel & coord)315 void reinit (const YGridComponent<Coordinates>& r, const iTupel& coord) 316 { 317 // initialize to given position in index set 318 for (int i=0; i<d; ++i) 319 _coord[i] = coord[i]; 320 321 // move superindex to first cell in subgrid 322 _superindex = 0; 323 for (int i=0; i<d; ++i) 324 _superindex += (r.offset(i)+coord[i]-r.origin(i))*r.superincrement(i); 325 326 _grid = &r; 327 } 328 329 //! Return true when two iterators over the same grid are equal (!). operator ==(const Iterator & i) const330 bool operator== (const Iterator& i) const 331 { 332 return _superindex == i._superindex; 333 } 334 335 //! Return true when two iterators over the same grid are not equal (!). operator !=(const Iterator & i) const336 bool operator!= (const Iterator& i) const 337 { 338 return _superindex != i._superindex; 339 } 340 341 //! Return consecutive index in enclosing grid superindex() const342 int superindex () const 343 { 344 return _superindex; 345 } 346 347 //! Return coordinate of the cell in direction i. coord(int i) const348 int coord (int i) const 349 { 350 return _coord[i]; 351 } 352 353 //! Return coordinate of the cell as reference (do not modify). coord() const354 const iTupel& coord () const 355 { 356 return _coord; 357 } 358 359 //! move this iterator dist cells in direction i move(int i,int dist)360 void move (int i, int dist) 361 { 362 _coord[i] += dist; 363 _superindex += dist*_grid->superincrement(i); 364 } 365 366 //! move this iterator dist cells in direction i move(const iTupel & dist)367 void move (const iTupel& dist) 368 { 369 for (int i = 0; i < d; ++i) 370 { 371 _coord[i] += dist[i]; 372 _superindex += dist[i]*_grid->superincrement(i); 373 } 374 } 375 376 //! Increment iterator to next cell with position. operator ++()377 Iterator& operator++ () 378 { 379 for (int i=0; i<d; i++) // check for wrap around 380 { 381 _superindex += _grid->superincrement(i); // move on cell in direction i 382 if (++_coord[i] <= _grid->max(i)) 383 return *this; 384 else 385 { 386 _coord[i] = _grid->origin(i); // move back to origin in direction i 387 _superindex -= _grid->size(i) * _grid->superincrement(i); 388 } 389 } 390 // if we wrapped around, back to to begin(), we must put the iterator to end() 391 if (_coord == _grid->origin()) 392 { 393 for (int i=0; i<d; i++) 394 _superindex += (_grid->size(i)-1) * _grid->superincrement(i); 395 _superindex += _grid->superincrement(0); 396 } 397 return *this; 398 } 399 400 //! Return ith component of lower left corner of the entity associated with the current coordinates and shift. lowerleft(int i) const401 ct lowerleft(int i) const 402 { 403 return _grid->getCoords()->coordinate(i,_coord[i]); 404 } 405 406 //! Return lower left corner of the entity associated with the current coordinates and shift. lowerleft() const407 fTupel lowerleft() const 408 { 409 fTupel ll; 410 for (int i=0; i<d; i++) 411 ll[i] = lowerleft(i); 412 return ll; 413 } 414 415 //! Return ith component of upper right corder of the entity associated with the current coordinates and shift. upperright(int i) const416 ct upperright(int i) const 417 { 418 int coord = _coord[i]; 419 if (shift(i)) 420 coord++; 421 return _grid->getCoords()->coordinate(i,coord); 422 } 423 424 //! Return upper right corder of the entity associated with the current coordinates and shift. upperright() const425 fTupel upperright() const 426 { 427 fTupel ur; 428 for (int i=0; i<d; i++) 429 ur[i] = upperright(i); 430 return ur; 431 } 432 433 //! Return meshsize in direction i meshsize(int i) const434 ct meshsize (int i) const 435 { 436 return _grid->getCoords()->meshsize(i,_coord[i]); 437 } 438 439 //! Return meshsize of current cell as reference. meshsize() const440 fTupel meshsize () const 441 { 442 fTupel h; 443 for (int i=0; i<d; i++) 444 h[i] = meshsize(i); 445 return h; 446 } 447 shift(int i) const448 bool shift (int i) const 449 { 450 return _grid->shift(i); 451 } 452 shift() const453 std::bitset<d> shift() const 454 { 455 return _grid->shift(); 456 } 457 coordCont() const458 Coordinates* coordCont() const 459 { 460 return _grid->getCoords(); 461 } 462 463 protected: 464 iTupel _coord; //!< current position in index set 465 int _superindex = 0; //!< consecutive index in enclosing grid 466 const YGridComponent<Coordinates>* _grid = nullptr; 467 }; 468 469 superindex(iTupel coord) const470 int superindex(iTupel coord) const 471 { 472 // move superindex to first cell in subgrid 473 int si = 0; 474 for (int i=0; i<d; ++i) 475 si += (offset(i)+coord[i]-origin(i))*_superincrement[i]; 476 return si; 477 } 478 superincrement(int i) const479 int superincrement(int i) const 480 { 481 return _superincrement[i]; 482 } 483 484 //! return iterator to first element of index set begin() const485 Iterator begin () const 486 { 487 return Iterator(*this); 488 } 489 490 //! return iterator to given element of index set begin(const iTupel & co) const491 Iterator begin (const iTupel& co) const 492 { 493 return Iterator(*this,co); 494 } 495 496 //! return subiterator to last element of index set end() const497 Iterator end () const 498 { 499 iTupel last; 500 for (int i=0; i<d; i++) 501 last[i] = max(i); 502 last[0] += 1; 503 return Iterator(*this,last); 504 } 505 506 private: 507 iTupel _origin; 508 std::bitset<d> _shift; 509 Coordinates* _coords; 510 iTupel _size; 511 iTupel _offset; //!< offset to origin of the enclosing grid 512 iTupel _supersize; //!< size of the enclosing grid 513 iTupel _superincrement; //!< moves consecutive index by one in this direction in supergrid 514 515 }; 516 517 518 //! Output operator for ygrids 519 template <class Coordinates> operator <<(std::ostream & s,YGridComponent<Coordinates> e)520 inline std::ostream& operator<< (std::ostream& s, YGridComponent<Coordinates> e) 521 { 522 s << "Printing YGridComponent structure:" << std::endl; 523 s << "Origin: " << e.origin() << std::endl; 524 s << "Shift: " << e.shift() << std::endl; 525 s << "Size: " << e.size() << std::endl; 526 s << "Offset: " << e.offset() << std::endl; 527 s << "Supersize: " << e.supersize() << std::endl; 528 return s; 529 } 530 531 //! Output operator for ygrids 532 template <class Coordinates> operator <<(std::ostream & s,typename YGridComponent<Coordinates>::Iterator & e)533 inline std::ostream& operator<< (std::ostream& s, typename YGridComponent<Coordinates>::Iterator& e) 534 { 535 s << "Printing YGridComponent Iterator:" << std::endl << "Iterator at " << e.coord() << " (index "; 536 s << e.index() << "), position " << e.position(); 537 return s; 538 } 539 540 /** \brief implements a collection of YGridComponents which form a codimension 541 * Entities of given codimension c need to be represented by d choose c YgridComponents. 542 * All entities in one such component share the same set of spanning unit vectors. 543 * A YGrid is used to iterate over the entire set of components the codimension 544 * consists of. It doesn't hold any data, but instead holds an iterator range into 545 * an array of components (which is owned by YGridLevel). 546 */ 547 template<class Coordinates> 548 class YGrid 549 { 550 public: 551 static const int dim = Coordinates::dimension; 552 553 // define data array iterator 554 typedef YGridComponent<Coordinates>* DAI; 555 556 typedef typename std::array<int, dim> iTupel; 557 558 //! set start iterator in the data array setBegin(DAI begin)559 void setBegin(DAI begin) 560 { 561 _begin = begin; 562 } 563 564 //! get which component belongs to a given shift vector shiftmapping(const std::bitset<dim> & shift) const565 int shiftmapping(const std::bitset<dim>& shift) const 566 { 567 return _shiftmapping[shift.to_ulong()]; 568 } 569 570 //! get start iterator in the data array dataBegin() const571 DAI dataBegin() const 572 { 573 return _begin; 574 } 575 576 //! get end iterator in the data array dataEnd() const577 DAI dataEnd() const 578 { 579 return _end; 580 } 581 582 //! decide whether a coordinate is in the grid (depending on the component) inside(const iTupel & coord,const std::bitset<dim> & shift=std::bitset<dim> ()) const583 bool inside(const iTupel& coord, const std::bitset<dim>& shift = std::bitset<dim>()) const 584 { 585 return (_begin+_shiftmapping[shift.to_ulong()])->inside(coord); 586 } 587 588 /** \brief Iterator over a collection o YGrids 589 * A YGrid::Iterator is the heart of an entity in YaspGrid. 590 */ 591 class Iterator 592 { 593 public: 594 595 //! default constructor 596 Iterator () = default; 597 598 //! construct an iterator from coordinates and component Iterator(const YGrid<Coordinates> & yg,const std::array<int,dim> & coords,int which=0)599 Iterator (const YGrid<Coordinates>& yg, const std::array<int,dim>& coords, int which = 0) 600 : _which(which), _yg(&yg) 601 { 602 _it = typename YGridComponent<Coordinates>::Iterator(*(_yg->dataBegin()+which),coords); 603 } 604 605 //! create an iterator to start or end of the codimension Iterator(const YGrid<Coordinates> & yg,bool end=false)606 Iterator (const YGrid<Coordinates>& yg, bool end=false) : _yg(&yg) 607 { 608 if (end) 609 { 610 _it = _yg->_itends.back(); 611 _which = _yg->_itends.size() - 1; 612 } 613 else 614 { 615 _it = _yg->_itbegins[0]; 616 _which = 0; 617 } 618 } 619 620 //! reinitializes an iterator, as if it was just constructed. reinit(const YGrid<Coordinates> & yg,const std::array<int,dim> & coords,int which=0)621 void reinit(const YGrid<Coordinates>& yg, const std::array<int,dim>& coords, int which = 0) 622 { 623 _yg = &yg; 624 _which = which; 625 _it = typename YGridComponent<Coordinates>::Iterator(*(_yg->dataBegin()+which),coords); 626 } 627 628 //! return coordinate at the current position (direction i) coord(int i) const629 int coord (int i) const 630 { 631 return _it.coord(i); 632 } 633 634 //! return coordinate array at the current position coord() const635 const std::array<int, dim>& coord () const 636 { 637 return _it.coord(); 638 } 639 lowerleft(int i) const640 typename Coordinates::ctype lowerleft(int i) const 641 { 642 return _it.lowerleft(i); 643 } 644 lowerleft() const645 Dune::FieldVector<typename Coordinates::ctype,dim> lowerleft() const 646 { 647 return _it.lowerleft(); 648 } 649 upperright(int i) const650 typename Coordinates::ctype upperright(int i) const 651 { 652 return _it.upperright(i); 653 } 654 upperright() const655 Dune::FieldVector<typename Coordinates::ctype,dim> upperright() const 656 { 657 return _it.upperright(); 658 } 659 660 //! return the current meshsize in direction i meshsize(int i) const661 typename Coordinates::ctype meshsize (int i) const 662 { 663 return _it.meshsize(i); 664 } 665 666 //! return the current meshsize vector meshsize() const667 Dune::FieldVector<typename Coordinates::ctype,dim> meshsize() const 668 { 669 return _it.meshsize(); 670 } 671 672 //! return the shift in direction i shift(int i) const673 bool shift (int i) const 674 { 675 return _it.shift(i); 676 } 677 678 //! return the shift vector shift() const679 std::bitset<dim> shift () const 680 { 681 return _it.shift(); 682 } 683 684 //! return the superindex superindex() const685 int superindex() const 686 { 687 // the offset of the current component has to be taken into account 688 return _yg->_indexOffset[_which] + _it.superindex(); 689 } 690 691 //! increment to the next entity jumping to next component if necessary operator ++()692 Iterator& operator++ () 693 { 694 if ((++_it == _yg->_itends[_which]) && (_which < _yg->_itends.size()-1)) 695 _it = _yg->_itbegins[++_which]; 696 return *this; 697 } 698 699 //! compare two iterators: component has to match operator ==(const Iterator & i) const700 bool operator==(const Iterator& i) const 701 { 702 if (_which != i._which) 703 return false; 704 return _it == i._it; 705 } 706 707 //! compare two iterators: component has to match operator !=(const Iterator & i) const708 bool operator!=(const Iterator& i) const 709 { 710 if (_it != i._it) 711 return true; 712 return _which != i._which; 713 } 714 715 //! return the current component number which() const716 int which() const 717 { 718 return _which; 719 } 720 721 //! move the grid, this is only done and needed for codim 0 move(int i,int dist)722 void move(int i, int dist) 723 { 724 _it.move(i,dist); 725 } 726 move(const iTupel & dist)727 void move(const iTupel& dist) 728 { 729 _it.move(dist); 730 } 731 coordCont() const732 Coordinates* coordCont() const 733 { 734 return _it.coordCont(); 735 } 736 737 738 private: 739 unsigned int _which = 0; 740 const YGrid<Coordinates>* _yg = nullptr; 741 typename YGridComponent<Coordinates>::Iterator _it; 742 }; 743 744 //! return begin iterator for the codimension and partition the ygrid represents begin() const745 Iterator begin() const 746 { 747 return Iterator(*this); 748 } 749 750 //! return iterator pointint to a specified position begin(const std::array<int,dim> & coord,int which=0) const751 Iterator begin(const std::array<int, dim>& coord, int which = 0) const 752 { 753 return Iterator(*this, coord, which); 754 } 755 756 //! return end iterator for the codimension and partition the ygrid represents end() const757 Iterator end() const 758 { 759 return Iterator(*this,true); 760 } 761 superindex(const iTupel & coord,int which) const762 int superindex(const iTupel& coord, int which) const 763 { 764 return _indexOffset[which] + (dataBegin()+which)->superindex(coord); 765 } 766 767 768 // finalize the ygrid construction by storing component iterators finalize(const DAI & end,int artificialOffset=0)769 void finalize(const DAI& end, int artificialOffset = 0) 770 { 771 // set the end iterator in the ygrid component array 772 _end = end; 773 774 _indexOffset.push_back(artificialOffset); 775 int k = 0; 776 for (DAI i=_begin; i != _end; ++i, ++k) 777 { 778 //store begin and end iterators 779 _itbegins.push_back(i->begin()); 780 _itends.push_back(i->end()); 781 782 // store index offset 783 _indexOffset.push_back(_indexOffset.back() + i->totalsize()); 784 785 // store shift to component mapping 786 _shiftmapping[i->shift().to_ulong()] = k; 787 } 788 _indexOffset.resize(_itends.size()); 789 } 790 791 private: 792 793 friend class YGrid<Coordinates>::Iterator; 794 DAI _begin; 795 DAI _end; 796 std::array<int,StaticPower<2,dim>::power> _shiftmapping; 797 std::vector<typename YGridComponent<Coordinates>::Iterator> _itbegins; 798 std::vector<typename YGridComponent<Coordinates>::Iterator> _itends; 799 std::vector<int> _indexOffset; 800 }; 801 802 //! Output operator for ygrids 803 template <class Coordinates> operator <<(std::ostream & s,const YGrid<Coordinates> & e)804 inline std::ostream& operator<< (std::ostream& s, const YGrid<Coordinates>& e) 805 { 806 s << "Printing YGrid structure:" << std::endl; 807 for (auto it = e.dataBegin(); it != e.dataEnd(); ++it) 808 s << *it << std::endl; 809 return s; 810 } 811 812 /** \brief implements a collection of multiple std::deque<Intersection> 813 * Intersections with neighboring processors are stored as std::deque<Intersection>. 814 * Eachsuch intersection only holds one YGridComponent. To do all communication 815 * associated with one codimension, multiple such deques have to be concatenated. 816 * YGridList manges this concatenation. As for YGrids, YGridList doesn't hold any 817 * data, but an iterator range into a data array owned by YGridLevel. 818 */ 819 template<class Coordinates> 820 class YGridList 821 { 822 public: 823 static const int dim = Coordinates::dimension; 824 825 /** \brief type describing an intersection with a neighboring processor */ 826 struct Intersection 827 { 828 /** \brief The intersection as a subgrid of the local grid */ 829 YGridComponent<Coordinates> grid; 830 /** \brief Rank of the process where the other grid is stored */ 831 int rank; 832 /** \brief Manhattan distance to the other grid */ 833 int distance; 834 /** \brief a YGrid stub, that acts wraps above YGrid Component and handels the index offset */ 835 YGrid<Coordinates> yg; 836 }; 837 838 // define data array iterator type 839 typedef typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator DAI; 840 841 // iterator that allows to iterate over a concatenation of deques. namely those 842 // that belong to the same codimension. 843 class Iterator 844 { 845 public: 846 847 //! return iterator to begin and end of the container Iterator(const YGridList<Coordinates> & ygl,bool end=false)848 Iterator(const YGridList<Coordinates>& ygl, bool end=false) : _end(ygl.dataEnd()), _which(ygl.dataBegin()) 849 { 850 _it = _which->begin(); 851 852 // advance the iterator to the first element that exists. 853 // some deques might be empty and should be skipped 854 while ((_which != _end) && (_it == _which->end())) 855 { 856 ++_which; 857 if (_which != _end) 858 _it = _which->begin(); 859 } 860 // the iterator is at the end if and only if _which==_end 861 if (end) 862 { 863 _which = _end; 864 } 865 } 866 867 //! increment iterator operator ++()868 Iterator& operator++ () 869 { 870 ++_it; 871 // advance the iterator to the next element that exists. 872 // some deques might be empty and should be skipped 873 while ((_which != _end) && (_it == _which->end())) 874 { 875 ++_which; 876 if (_which != _end) 877 _it = _which->begin(); 878 } 879 return *this; 880 } 881 882 //! dereference iterator operator ->() const883 typename std::deque<Intersection>::iterator operator->() const 884 { 885 return _it; 886 } 887 888 //! dereference iterator operator *() const889 typename std::deque<Intersection>::iterator operator*() const 890 { 891 return _it; 892 } 893 894 //! compare two iterators operator ==(const Iterator & i) const895 bool operator== (const Iterator& i) const 896 { 897 if (_which != i._which) 898 return false; 899 if (_which == _end) 900 return true; 901 return _it == i._it; 902 } 903 904 //! compare two iterators operator !=(const Iterator & i) const905 bool operator!= (const Iterator& i) const 906 { 907 if (_which != i._which) 908 return true; 909 if (_which == _end) 910 return false; 911 return _it != i._it; 912 } 913 914 private: 915 typename std::deque<Intersection>::iterator _it; 916 DAI _end; 917 DAI _which; 918 }; 919 920 //! return iterator pointing to the begin of the container begin() const921 Iterator begin() const 922 { 923 return Iterator(*this); 924 } 925 926 //! return iterator pointing to the end of the container end() const927 Iterator end() const 928 { 929 return Iterator(*this,true); 930 } 931 932 //! set start iterator in the data array setBegin(typename std::array<std::deque<Intersection>,StaticPower<2,dim>::power>::iterator begin)933 void setBegin(typename std::array<std::deque<Intersection>, StaticPower<2,dim>::power>::iterator begin) 934 { 935 _begin = begin; 936 } 937 938 //! get start iterator in the data array dataBegin() const939 DAI dataBegin() const 940 { 941 return _begin; 942 } 943 944 //! get end iterator in the data array dataEnd() const945 DAI dataEnd() const 946 { 947 return _end; 948 } 949 950 //! return the size of the container, this is the sum of the sizes of all deques size() const951 int size() const 952 { 953 int count = 0; 954 for (DAI it = _begin; it != _end; ++it) 955 count += it->size(); 956 return count; 957 } 958 959 //! finalize the YGridLIst finalize(DAI end,const YGrid<Coordinates> & ygrid)960 void finalize(DAI end, const YGrid<Coordinates>& ygrid) 961 { 962 // Instead of directly iterating over the intersection deques, this code 963 // iterates over the components of an associated ygrid and works its way 964 // through the list of intersection deques in parallel. 965 // The reason for this convoluted iteration technique is that there are not 966 // necessarily intersections for all possible shifts, but we have to make 967 // sure that we stop at each shift to update the per-component index shift. 968 // This is ensured by iterating over the ygrid, which is guaranteed to have 969 // a component for each shift vector. 970 971 // set end iterator in the data array 972 _end = end; 973 974 //! set offsets allow the YGridComponents in the Intersctions to behave as YGrids 975 int offset = 0; 976 977 DAI i = _begin; 978 979 // make sure that we have a valid deque (i.e. a non-empty one) 980 while (i != _end && i->begin() == i->end()) 981 ++i; 982 983 for (auto yit = ygrid.dataBegin(); yit != ygrid.dataEnd(); ++yit) 984 { 985 if (i == _end) 986 break; 987 auto it = i->begin(); 988 if (it->grid.shift() == yit->shift()) 989 { 990 // iterate over the intersections in the deque and set the offset 991 for (; it != i->end(); ++it) 992 { 993 it->yg.setBegin(&(it->grid)); 994 it->yg.finalize(&(it->grid)+1, offset); 995 } 996 997 // advance to next non-empty deque 998 ++i; 999 while (i != _end && i->begin() == i->end()) 1000 ++i; 1001 } 1002 1003 // update the offset from the ygrid component 1004 int add = 1; 1005 for (int j=0; j<dim; j++) 1006 add *= yit->supersize(j); 1007 offset += add; 1008 } 1009 assert (i == end); 1010 } 1011 1012 private: 1013 DAI _begin; 1014 DAI _end; 1015 }; 1016 1017 } // namespace Dune 1018 1019 #endif 1020