Dune Core Modules (2.8.0)

yaspgridentity.hh
Go to the documentation of this file.
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_YASPGRIDENTITY_HH
4#define DUNE_GRID_YASPGRIDENTITY_HH
5
19//========================================================================
20
21
22
23
24namespace Dune {
25
26 namespace Yasp {
27
28#ifndef DOXYGEN
29
30 // table for quick evaluation of binomial coefficients
31 template<int n>
32 struct BinomialTable
33 {
34 // evaluation - note that in general d!=n, n is only the
35 // maximum value of d (in our case dimworld)
36 static constexpr int evaluate(int d, int c)
37 {
38 return _values[_offsets[d] + c];
39 }
40
41 // the actual implementation
42 static constexpr int binomial(int d, int c)
43 {
44 long binomial=1;
45 for (int i=d-c+1; i<=d; i++)
46 binomial *= i;
47 for (long i=2; i<=c; i++)
48 binomial /= i;
49 return binomial;
50 }
51
52 private:
53 // prevent construction
54 BinomialTable() = delete;
55
56 // compute binomial(r, c) and advance row `r` and column `c`
57 static constexpr int nextValue(int& r, int& c)
58 {
59 const auto result = binomial(r, c);
60
61 c += 1;
62 if (c > r) {
63 r += 1;
64 c = 0;
65 }
66
67 return result;
68 }
69
70 template<std::size_t... I>
71 static constexpr std::array<int, sizeof...(I)> computeValues(std::index_sequence<I...>)
72 {
73 int r = 0, c = 0;
74 return {{ ((void)I, nextValue(r, c))... }};
75 }
76
77 template<std::size_t... I>
78 static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
79 { return {{ (I*(I+1)/2)... }}; }
80
81 static constexpr std::array<int,(n+1)*(n+2)/2> _values = computeValues(std::make_index_sequence<(n+1)*(n+2)/2>{});
82 static constexpr std::array<int,n+1> _offsets = computeOffsets(std::make_index_sequence<n+1>{});
83 };
84
85#if __cplusplus < 201703L
86 template<int n>
87 constexpr std::array<int,(n+1)*(n+2)/2> BinomialTable<n>::_values;
88 template<int n>
89 constexpr std::array<int,n+1> BinomialTable<n>::_offsets;
90#endif
91
98 template<int dimworld>
99 constexpr int subEnt(int d, int c)
100 {
101 return (d < c ? 0 : BinomialTable<dimworld>::evaluate(d,c) << c);
102 }
103
104 // Make a table mapping all subentities of a codim 0 entity to a value.
105 // F is the functor to be evaluated.
106 template<typename F, int dim>
107 struct EntityShiftTable
108 {
109 typedef std::bitset<dim> value_type;
110
111 static value_type evaluate(int i, int codim)
112 {
113 return {_values[_offsets[codim] + i]};
114 }
115
116 private:
117
118 // prevent construction
119 EntityShiftTable() = delete;
120
121 // compute offset of codimension `codim` entities and advance `offset`
122 static constexpr int nextOffset(int& offset, int codim)
123 {
124 if (codim == 0) {
125 offset = 0;
126 return 0;
127 }
128
129 offset += subEnt<dim>(dim, codim-1);
130 return offset;
131 }
132
133 template<std::size_t... I>
134 static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
135 {
136 int offset = 0;
137 return {{ (nextOffset(offset, I))... }};
138 }
139
140 // compute shift table entry for (`codim`, `i`) and advance `codim`, `i`
141 static constexpr unsigned char nextValue(int& codim, int& i)
142 {
143 const auto result = F::evaluate(i, codim);
144
145 i += 1;
146 if (i >= subEnt<dim>(dim, codim)) {
147 codim += 1;
148 i = 0;
149 }
150
151 return result;
152 }
153
154 template<std::size_t... I>
155 static constexpr std::array<unsigned char, sizeof...(I)> computeValues(std::index_sequence<I...>)
156 {
157 int codim = 0, i = 0;
158 return {{ ((void)I, nextValue(codim, i))... }};
159 }
160
161 static constexpr std::array<int,dim+1> _offsets = computeOffsets(std::make_index_sequence<dim+1>{});
162 static constexpr std::array<unsigned char,StaticPower<3,dim>::power> _values = computeValues(std::make_index_sequence<StaticPower<3,dim>::power>{});
163
164 };
165
166#if __cplusplus < 201703L
167 template<typename F, int dim>
168 constexpr std::array<int,dim+1> EntityShiftTable<F, dim>::_offsets;
169 template<typename F, int dim>
170 constexpr std::array<unsigned char,StaticPower<3,dim>::power> EntityShiftTable<F, dim>::_values;
171#endif
172
173 // functor for doing the actual entity shift calculation
174 template<int dim>
175 struct calculate_entity_shift
176 {
177 static constexpr unsigned long long evaluate(int index, int cc)
178 {
179 auto result = 0ull;
180 for (int d = dim; d>0; d--)
181 {
182 if (cc == d)
183 return result;
184 if (index < subEnt<dim>(d-1,cc))
185 result |= 1ull << (d-1);
186 else
187 {
188 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
189 cc--;
190 }
191 }
192 return result;
193 }
194 };
195
204 template<int dim>
205 std::bitset<dim> entityShift(int index, int cc)
206 {
207 return EntityShiftTable<calculate_entity_shift<dim>,dim>::evaluate(index,cc);
208 }
209
210 // functor for doing the actual entity move calculation
211 template<int dim>
212 struct calculate_entity_move
213 {
214 static constexpr unsigned long long evaluate(int index, int cc)
215 {
216 auto result = 0ull;
217 for (int d = dim; d>0; d--)
218 {
219 if (d == cc)
220 {
221 // result[d-1] = index & (1<<(d-1));
222 result &= ~(1ull << (d-1));
223 result |= index & (1ull << (d-1));
224
225 index &= ~(1<<(d-1));
226 }
227 if (index >= subEnt<dim>(d-1,cc))
228 {
229 if ((index - subEnt<dim>(d-1,cc)) / subEnt<dim>(d-1,cc-1) == 1)
230 {
231 result |= 1ull << (d-1);
232 }
233 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
234 cc--;
235 }
236 }
237 return result;
238 }
239
240 };
241
249 template<int dim>
250 std::bitset<dim> entityMove(int index, int cc)
251 {
252 return EntityShiftTable<calculate_entity_move<dim>,dim>::evaluate(index,cc);
253 }
254
255#endif //DOXYGEN
256
257 } // namespace Yasp.
258
259 template<int codim, int dim, class GridImp>
260 class YaspEntity
261 : public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
262 {
263
264 template<int, PartitionIteratorType, typename>
265 friend class YaspLevelIterator;
266
267 public:
268 typedef typename GridImp::ctype ctype;
269
270 typedef typename GridImp::template Codim<codim>::Geometry Geometry;
271 typedef typename GridImp::Traits::template Codim<codim>::GeometryImpl GeometryImpl;
272
273 typedef typename GridImp::template Codim<codim>::EntitySeed EntitySeed;
274
276 int level () const
277 {
278 return _g->level();
279 }
280
284 EntitySeed seed() const
285 {
286 return EntitySeed(YaspEntitySeed<codim,GridImp>(_g->level(), _it.coord(), _it.which()));
287 }
288
290 Geometry geometry () const
291 {
292 GeometryImpl _geometry(_it.lowerleft(),_it.upperright(),_it.shift());
293 return Geometry(_geometry);
294 }
295
300 unsigned int subEntities (unsigned int cc) const
301 {
302 return Dune::Yasp::subEnt<dim>(dim-codim,cc-codim);
303 }
304
306 PartitionType partitionType () const
307 {
308 if (_g->interior[codim].inside(_it.coord(),_it.shift()))
309 return InteriorEntity;
310 if (_g->interiorborder[codim].inside(_it.coord(),_it.shift()))
311 return BorderEntity;
312 if (_g->overlap[codim].inside(_it.coord(),_it.shift()))
313 return OverlapEntity;
314 if (_g->overlapfront[codim].inside(_it.coord(),_it.shift()))
315 return FrontEntity;
316 return GhostEntity;
317 }
318
319 typedef typename GridImp::YGridLevelIterator YGLI;
320 typedef typename GridImp::YGrid::Iterator I;
321 YaspEntity ()
322 {}
323
324 YaspEntity (const YGLI& g, const I& it)
325 : _it(it), _g(g)
326 {}
327
328 YaspEntity (YGLI&& g, const I&& it)
329 : _it(std::move(it)), _g(std::move(g))
330 {}
331
333 bool equals (const YaspEntity& e) const
334 {
335 return _it == e._it && _g == e._g;
336 }
337
338 // IndexSets needs access to the private index methods
339 friend class Dune::YaspIndexSet<GridImp,true>;
340 friend class Dune::YaspIndexSet<GridImp,false>;
341 friend class Dune::YaspGlobalIdSet<GridImp>;
342 typedef typename GridImp::PersistentIndexType PersistentIndexType;
343
345 PersistentIndexType persistentIndex () const
346 {
347 // get size of global grid (in elements)
348 std::array<int,dim> size;
349
350 for (int i=0; i<dim; i++)
351 {
352 // correct size according to shift
353 size[i] = _g->mg->levelSize(_g->level(), i);
354 if (!_it.shift(i))
355 size[i]++;
356 }
357
358 // encode codim
359 PersistentIndexType id(_it.shift().to_ulong());
360
361 // encode level
362 id = id << yaspgrid_level_bits;
363 id = id+PersistentIndexType(_g->level());
364
365 // encode coordinates
366 for (int i=dim-1; i>=0; i--)
367 {
368 id = id << yaspgrid_dim_bits;
369 id = id+PersistentIndexType(_it.coord(i));
370 }
371
372 return id;
373 }
374
376 int compressedIndex () const
377 {
378 return _it.superindex();
379 }
380
382 int subCompressedIndex (int i, unsigned int cc) const
383 {
384 // get the shift of the entity and the subentity
385 // the subentity shift is only available in the space spanned by the entity
386 std::bitset<dim-codim> subent_shift = Dune::Yasp::entityShift<dim-codim>(i,cc-codim);
387 std::bitset<dim-codim> subent_move = Dune::Yasp::entityMove<dim-codim>(i,cc-codim);
388
389 std::bitset<dim> shift = _it.shift();
390 std::array<int, dim> coord = _it.coord();
391 for (int j=0, k=0; j<dim; j++)
392 {
393 if (!shift[j])
394 continue;
395
396 coord[j] += subent_move[k];
397 shift[j] = subent_shift[k];
398 k++;
399 }
400
401 int which = _g->overlapfront[cc].shiftmapping(shift);
402 return _g->overlapfront[cc].superindex(coord,which);
403 }
404 public:
405 const I& transformingsubiterator() const { return _it; }
406 const YGLI& gridlevel() const { return _g; }
407 I& transformingsubiterator() { return _it; }
408 YGLI& gridlevel() { return _g; }
409 const GridImp * yaspgrid() const { return _g->mg; }
410 protected:
411 I _it; // position in the grid level
412 YGLI _g; // access to grid level
413 };
414
415
416 // specialization for codim=0
417 template<int dim, class GridImp>
418 class YaspEntity<0,dim,GridImp>
419 : public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
420 {
421 enum { dimworld = GridImp::dimensionworld };
422
423 typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
424
425 template<int, PartitionIteratorType, typename>
426 friend class YaspLevelIterator;
427
428 template<typename>
429 friend class YaspHierarchicIterator;
430
431 public:
432 typedef typename GridImp::ctype ctype;
433
434 typedef typename GridImp::YGridLevelIterator YGLI;
435 typedef typename GridImp::YGrid::Iterator I;
436
437 typedef typename GridImp::template Codim< 0 >::Geometry Geometry;
438 typedef typename GridImp::template Codim< 0 >::LocalGeometry LocalGeometry;
439
440 template <int cd>
441 struct Codim
442 {
443 typedef typename GridImp::template Codim<cd>::Entity Entity;
444 };
445
446 typedef typename GridImp::template Codim<0>::Entity Entity;
447 typedef typename GridImp::template Codim<0>::EntitySeed EntitySeed;
448 typedef typename GridImp::LevelIntersectionIterator IntersectionIterator;
449 typedef typename GridImp::LevelIntersectionIterator LevelIntersectionIterator;
450 typedef typename GridImp::LeafIntersectionIterator LeafIntersectionIterator;
451 typedef typename GridImp::HierarchicIterator HierarchicIterator;
452
454 typedef typename GridImp::PersistentIndexType PersistentIndexType;
455
457 typedef typename GridImp::YGrid::iTupel iTupel;
458
459 // constructor
460 YaspEntity ()
461 {}
462
463 YaspEntity (const YGLI& g, const I& it)
464 : _it(it), _g(g)
465 {}
466
467 YaspEntity (const YGLI& g, I&& it)
468 : _it(std::move(it)), _g(g)
469 {}
470
471 YaspEntity (YGLI&& g, I&& it)
472 : _it(std::move(it)), _g(std::move(g))
473 {}
474
476 bool equals (const YaspEntity& e) const
477 {
478 return _it == e._it && _g == e._g;
479 }
480
482 int level () const { return _g->level(); }
483
487 EntitySeed seed () const {
488 return EntitySeed(YaspEntitySeed<0,GridImp>(_g->level(), _it.coord()));
489 }
490
492 PartitionType partitionType () const
493 {
494 if (_g->interior[0].inside(_it.coord(),_it.shift()))
495 return InteriorEntity;
496 if (_g->overlap[0].inside(_it.coord(),_it.shift()))
497 return OverlapEntity;
498 DUNE_THROW(GridError, "Impossible GhostEntity");
499 return GhostEntity;
500 }
501
503 Geometry geometry () const {
504 // the element geometry
505 auto ll = _it.lowerleft();
506 auto ur = _it.upperright();
507
508 // If on periodic overlap, transform coordinates by domain size
509 for (int i=0; i<dimworld; i++) {
510 if (gridlevel()->mg->isPeriodic(i)) {
511 int coord = transformingsubiterator().coord(i);
512 if (coord < 0) {
513 auto size = _g->mg->domainSize()[i];
514 ll[i] += size;
515 ur[i] += size;
516 } else if (coord + 1 > gridlevel()->mg->levelSize(gridlevel()->level(),i)) {
517 auto size = _g->mg->domainSize()[i];
518 ll[i] -= size;
519 ur[i] -= size;
520 }
521 }
522 }
523
524 GeometryImpl _geometry(ll,ur);
525 return Geometry( _geometry );
526 }
527
532 template<int cc> int count () const
533 {
534 return Dune::Yasp::subEnt<dim>(dim,cc);
535 }
536
541 unsigned int subEntities (unsigned int codim) const
542 {
543 return Dune::Yasp::subEnt<dim>(dim,codim);
544 }
545
548 template<int cc>
549 typename Codim<cc>::Entity subEntity (int i) const
550 {
551 // calculate move bitset
552 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
553
554 // get the coordinate and modify it
555 iTupel coord = _it.coord();
556 for (int j=0; j<dim; j++)
557 if (move[j])
558 coord[j]++;
559
560 int which = _g->overlapfront[cc].shiftmapping(Dune::Yasp::entityShift<dim>(i,cc));
561 return typename Codim<cc>::Entity(YaspEntity<cc,GridImp::dimension,GridImp>(_g,_g->overlapfront[cc].begin(coord, which)));
562 }
563
565 Entity father () const
566 {
567 // check if coarse level exists
568 if (_g->level()<=0)
569 DUNE_THROW(GridError, "tried to call father on level 0");
570
571 // yes, get iterator to it
572 YGLI cg(_g);
573 --cg;
574
575 // coordinates of the cell
576 iTupel coord = _it.coord();
577
578 // get coordinates on next coarser level
579 for (int k=0; k<dim; k++) coord[k] = coord[k]/2;
580
581 return Entity(YaspEntity<0,GridImp::dimension,GridImp>(cg,cg->overlap[0].begin(coord)));
582 }
583
585 bool hasFather () const
586 {
587 return (_g->level()>0);
588 }
589
592 LocalGeometry geometryInFather () const
593 {
594 // configure one of the 2^dim transformations
595 FieldVector<ctype,dim> ll(0.0),ur(0.5);
596
597 for (int k=0; k<dim; k++)
598 {
599 if (_it.coord(k)%2)
600 {
601 ll[k] = 0.5;
602 ur[k] = 1.0;
603 }
604 }
605
606 return LocalGeometry( YaspGeometry<dim,dim,GridImp>(ll,ur) );
607 }
608
609 const I& transformingsubiterator () const { return _it; }
610 const YGLI& gridlevel () const { return _g; }
611 I& transformingsubiterator() { return _it; }
612 YGLI& gridlevel() { return _g; }
613 const GridImp* yaspgrid () const { return _g->mg; }
614
615 bool isLeaf() const
616 {
617 return (_g->level() == yaspgrid()->maxLevel());
618 }
619
622 bool isNew () const { return yaspgrid()->adaptRefCount > 0 && yaspgrid()->maxLevel() < _g->level() + yaspgrid()->adaptRefCount; }
623
626 bool mightVanish () const { return false; }
627
629 IntersectionIterator ibegin () const
630 {
631 return YaspIntersectionIterator<GridImp>(*this,false);
632 }
633
635 LeafIntersectionIterator ileafbegin () const
636 {
637 // only if entity is leaf this iterator delivers intersections
638 return YaspIntersectionIterator<GridImp>(*this, ! isLeaf() );
639 }
640
642 LevelIntersectionIterator ilevelbegin () const
643 {
644 return ibegin();
645 }
646
648 IntersectionIterator iend () const
649 {
650 return YaspIntersectionIterator<GridImp>(*this,true);
651 }
652
654 LeafIntersectionIterator ileafend () const
655 {
656 return iend();
657 }
658
660 LevelIntersectionIterator ilevelend () const
661 {
662 return iend();
663 }
664
669 HierarchicIterator hbegin (int maxlevel) const
670 {
671 return YaspHierarchicIterator<GridImp>(_g,_it,maxlevel);
672 }
673
675 HierarchicIterator hend (int /* maxlevel */) const
676 {
677 return YaspHierarchicIterator<GridImp>(_g,_it,_g->level());
678 }
679
680 private:
681 // IndexSets needs access to the private index methods
682 friend class Dune::YaspIndexSet<GridImp,true>;
683 friend class Dune::YaspIndexSet<GridImp,false>;
684 friend class Dune::YaspGlobalIdSet<GridImp>;
685
687 PersistentIndexType persistentIndex () const
688 {
689 // encode codim
690 PersistentIndexType id(_it.shift().to_ulong());
691
692 // encode level
693 id = id << yaspgrid_level_bits;
694 id = id+PersistentIndexType(_g->level());
695
696
697 // encode coordinates
698 for (int i=dim-1; i>=0; i--)
699 {
700 id = id << yaspgrid_dim_bits;
701 id = id+PersistentIndexType(_it.coord(i));
702 }
703
704 return id;
705 }
706
708 int compressedIndex () const
709 {
710 return _it.superindex();
711 }
712
714 PersistentIndexType subPersistentIndex (int i, int cc) const
715 {
716 // calculate shift and move bitsets
717 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
718 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
719
720 int trailing = (cc == dim) ? 1000 : 0;
721
722 std::array<int,dim> size = _g->mg->levelSize(_g->level());
723 std::array<int, dim> coord = _it.coord();
724 for (int j=0; j<dim; j++)
725 {
726 // correct size according to shift
727 if (!shift[j])
728 size[j]++;
729
730 // move the coordinates to the cell on which the entity lives
731 if (move[j])
732 coord[j]++;
733 }
734
735 for (int j=0; j<dim; j++)
736 {
737 // in the codim==dim case, count trailing zeroes.
738 if (cc == dim)
739 {
740 int zeroes = 0;
741 for (int k=0; k<_g->level(); k++)
742 if (coord[j] & (1<<k))
743 break;
744 else
745 zeroes++;
746 trailing = std::min(trailing,zeroes);
747 }
748 }
749
750 // encode codim
751 PersistentIndexType id(shift.to_ulong());
752
753 // encode level
754 id = id << yaspgrid_level_bits;
755 id = id+PersistentIndexType(_g->level()-trailing);
756
757 // encode coordinates
758 for (int j=dim-1; j>=0; j--)
759 {
760 id = id << yaspgrid_dim_bits;
761 id = id+PersistentIndexType(coord[j]>>trailing);
762 }
763
764 return id;
765 }
766
768 int subCompressedIndex (int i, int cc) const
769 {
770 // get shift and move of the subentity in question
771 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
772 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
773
774 std::array<int, dim> coord = _it.coord();
775 for (int j=0; j<dim; j++)
776 coord[j] += move[j];
777
778 int which = _g->overlapfront[cc].shiftmapping(shift);
779 return _g->overlapfront[cc].superindex(coord,which);
780 }
781
782 I _it; // position in the grid level
783 YGLI _g; // access to grid level
784 };
785
786
787 // specialization for codim=dim (vertex)
788 template<int dim, class GridImp>
789 class YaspEntity<dim,dim,GridImp>
790 : public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
791 {
792 enum { dimworld = GridImp::dimensionworld };
793
794 template<int, PartitionIteratorType, typename>
795 friend class YaspLevelIterator;
796
797 typedef typename GridImp::Traits::template Codim<dim>::GeometryImpl GeometryImpl;
798
799 public:
800 typedef typename GridImp::ctype ctype;
801
802 typedef typename GridImp::YGridLevelIterator YGLI;
803 typedef typename GridImp::YGrid::Iterator I;
804
805 typedef typename GridImp::template Codim<dim>::Geometry Geometry;
806
807 typedef typename GridImp::template Codim<dim>::EntitySeed EntitySeed;
808
810 typedef typename GridImp::PersistentIndexType PersistentIndexType;
811
813 typedef typename GridImp::YGrid::iTupel iTupel;
814
815 // constructor
816 YaspEntity ()
817 {}
818
819 YaspEntity (const YGLI& g, const I& it)
820 : _it(it), _g(g)
821 {}
822
823 YaspEntity (YGLI&& g, I&& it)
824 : _it(std::move(it)), _g(std::move(g))
825 {}
826
828 bool equals (const YaspEntity& e) const
829 {
830 return _it == e._it && _g == e._g;
831 }
832
834 int level () const {return _g->level();}
835
839 EntitySeed seed () const {
840 return EntitySeed(YaspEntitySeed<dim,GridImp>(_g->level(), _it.coord(), _it.which()));
841 }
842
847 unsigned int subEntities (unsigned int cc) const
848 {
849 return Dune::Yasp::subEnt<dim>(dim-dim,cc-dim);
850 }
851
853 Geometry geometry () const {
854 GeometryImpl _geometry((_it).lowerleft());
855 return Geometry( _geometry );
856 }
857
859 PartitionType partitionType () const
860 {
861 if (_g->interior[dim].inside(_it.coord(),_it.shift()))
862 return InteriorEntity;
863 if (_g->interiorborder[dim].inside(_it.coord(),_it.shift()))
864 return BorderEntity;
865 if (_g->overlap[dim].inside(_it.coord(),_it.shift()))
866 return OverlapEntity;
867 if (_g->overlapfront[dim].inside(_it.coord(),_it.shift()))
868 return FrontEntity;
869 return GhostEntity;
870 }
871
873 int subCompressedIndex (int, unsigned int ) const
874 {
875 return compressedIndex();
876 }
877
878 private:
879 // IndexSets needs access to the private index methods
880 friend class Dune::YaspIndexSet<GridImp,true>;
881 friend class Dune::YaspIndexSet<GridImp,false>;
882 friend class Dune::YaspGlobalIdSet<GridImp>;
883
885 PersistentIndexType persistentIndex () const
886 {
887 // get coordinate and size of global grid
888 iTupel size = _g->mg->levelSize(_g->level());
889
890 for (int i=0; i<dim; i++)
891 {
892 // we have vertices, add 1 size to all directions
893 size[i]++;
894 }
895
896 // determine min number of trailing zeroes
897 int trailing = 1000;
898 for (int i=0; i<dim; i++)
899 {
900 // count trailing zeros
901 int zeros = 0;
902 for (int j=0; j<_g->level(); j++)
903 if (_it.coord(i)&(1<<j))
904 break;
905 else
906 zeros++;
907 trailing = std::min(trailing,zeros);
908 }
909
910 // determine the level of this vertex
911 int level = _g->level()-trailing;
912
913 // encode codim: shift vector of vertices is 0.
914 PersistentIndexType id(0);
915
916 // encode level
917 id = id << yaspgrid_level_bits;
918 id = id+PersistentIndexType(level);
919
920 // encode coordinates
921 for (int i=dim-1; i>=0; i--)
922 {
923 id = id << yaspgrid_dim_bits;
924 id = id+PersistentIndexType(_it.coord(i)>>trailing);
925 }
926
927 return id;
928 }
929
931 int compressedIndex () const { return _it.superindex();}
932
933 public:
934 const I& transformingsubiterator() const { return _it; }
935 const YGLI& gridlevel() const { return _g; }
936 I& transformingsubiterator() { return _it; }
937 YGLI& gridlevel() { return _g; }
938
939 const GridImp * yaspgrid() const { return _g->mg; }
940 protected:
941 I _it; // position in the grid level
942 YGLI _g; // access to grid level
943 };
944
945} // namespace Dune
946
947#endif // DUNE_GRID_YASPGRIDENTITY_HH
Dune::YaspGlobalIdSet
persistent, globally unique Ids
Definition: yaspgrididset.hh:23
Dune::YaspGlobalIdSet::id
IdType id(const typename std::remove_const< GridImp >::type::Traits::template Codim< cd >::Entity &e) const
get id of an entity
Definition: yaspgrididset.hh:42
Dune::YaspIndexSet
Implementation of Level- and LeafIndexSets for YaspGrid.
Definition: yaspgridindexsets.hh:23
Dune::StaticPower::power
static constexpr int power
power stores m^p
Definition: power.hh:26
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::PartitionType
PartitionType
Attributes used in the generic overlap model.
Definition: gridenums.hh:28
Dune::FrontEntity
@ FrontEntity
on boundary between overlap and ghost
Definition: gridenums.hh:32
Dune::InteriorEntity
@ InteriorEntity
all interior entities
Definition: gridenums.hh:29
Dune::GhostEntity
@ GhostEntity
ghost entities
Definition: gridenums.hh:33
Dune::BorderEntity
@ BorderEntity
on boundary between interior and overlap
Definition: gridenums.hh:30
Dune::OverlapEntity
@ OverlapEntity
all entities lying in the overlap zone
Definition: gridenums.hh:31
Dune::Hybrid::equals
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:400
Dune::Simd::Overloads::min
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
Dune
Dune namespace.
Definition: alignedallocator.hh:11
Dune::binomial
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:128
Dune::power
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
std
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)