DUNE-FEM (unstable)

yaspgridentity.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_GRID_YASPGRIDENTITY_HH
6#define DUNE_GRID_YASPGRIDENTITY_HH
7
8#include <dune/common/math.hh>
9#include <dune/geometry/type.hh>
10
24//========================================================================
25
26
27
28
29namespace Dune {
30
31 namespace Yasp {
32
33#ifndef DOXYGEN
34
35 // table for quick evaluation of binomial coefficients
36 template<int n>
37 struct BinomialTable
38 {
39 // evaluation - note that in general d!=n, n is only the
40 // maximum value of d (in our case dimworld)
41 static constexpr int evaluate(int d, int c)
42 {
43 return _values[_offsets[d] + c];
44 }
45
46 private:
47 // prevent construction
48 BinomialTable() = delete;
49
50 // compute binomial(r, c) and advance row `r` and column `c`
51 static constexpr int nextValue(int& r, int& c)
52 {
53 const auto result = Dune::binomial(r, c);
54
55 c += 1;
56 if (c > r) {
57 r += 1;
58 c = 0;
59 }
60
61 return result;
62 }
63
64 template<std::size_t... I>
65 static constexpr std::array<int, sizeof...(I)> computeValues(std::index_sequence<I...>)
66 {
67 int r = 0, c = 0;
68 return {{ ((void)I, nextValue(r, c))... }};
69 }
70
71 template<std::size_t... I>
72 static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
73 { return {{ (I*(I+1)/2)... }}; }
74
75 static constexpr std::array<int,(n+1)*(n+2)/2> _values = computeValues(std::make_index_sequence<(n+1)*(n+2)/2>{});
76 static constexpr std::array<int,n+1> _offsets = computeOffsets(std::make_index_sequence<n+1>{});
77 };
78
85 template<int dimworld>
86 constexpr int subEnt(int d, int c)
87 {
88 return (d < c ? 0 : BinomialTable<dimworld>::evaluate(d,c) << c);
89 }
90
91 // Make a table mapping all subentities of a codim 0 entity to a value.
92 // F is the functor to be evaluated.
93 template<typename F, int dim>
94 struct EntityShiftTable
95 {
96 typedef std::bitset<dim> value_type;
97
98 static value_type evaluate(int i, int codim)
99 {
100 return {_values[_offsets[codim] + i]};
101 }
102
103 private:
104
105 // prevent construction
106 EntityShiftTable() = delete;
107
108 // compute offset of codimension `codim` entities and advance `offset`
109 static constexpr int nextOffset(int& offset, int codim)
110 {
111 if (codim == 0) {
112 offset = 0;
113 return 0;
114 }
115
116 offset += subEnt<dim>(dim, codim-1);
117 return offset;
118 }
119
120 template<std::size_t... I>
121 static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
122 {
123 int offset = 0;
124 return {{ (nextOffset(offset, I))... }};
125 }
126
127 // compute shift table entry for (`codim`, `i`) and advance `codim`, `i`
128 static constexpr unsigned char nextValue(int& codim, int& i)
129 {
130 const auto result = F::evaluate(i, codim);
131
132 i += 1;
133 if (i >= subEnt<dim>(dim, codim)) {
134 codim += 1;
135 i = 0;
136 }
137
138 return result;
139 }
140
141 template<std::size_t... I>
142 static constexpr std::array<unsigned char, sizeof...(I)> computeValues(std::index_sequence<I...>)
143 {
144 int codim = 0, i = 0;
145 return {{ ((void)I, nextValue(codim, i))... }};
146 }
147
148 static constexpr std::array<int,dim+1> _offsets = computeOffsets(std::make_index_sequence<dim+1>{});
149 static constexpr std::array<unsigned char,Dune::power(3,dim)> _values = computeValues(std::make_index_sequence<Dune::power(3,dim)>{});
150
151 };
152
153 // functor for doing the actual entity shift calculation
154 template<int dim>
155 struct calculate_entity_shift
156 {
157 static constexpr unsigned long long evaluate(int index, int cc)
158 {
159 auto result = 0ull;
160 for (int d = dim; d>0; d--)
161 {
162 if (cc == d)
163 return result;
164 if (index < subEnt<dim>(d-1,cc))
165 result |= 1ull << (d-1);
166 else
167 {
168 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
169 cc--;
170 }
171 }
172 return result;
173 }
174 };
175
184 template<int dim>
185 std::bitset<dim> entityShift(int index, int cc)
186 {
187 return EntityShiftTable<calculate_entity_shift<dim>,dim>::evaluate(index,cc);
188 }
189
190 // functor for doing the actual entity move calculation
191 template<int dim>
192 struct calculate_entity_move
193 {
194 static constexpr unsigned long long evaluate(int index, int cc)
195 {
196 auto result = 0ull;
197 for (int d = dim; d>0; d--)
198 {
199 if (d == cc)
200 {
201 // result[d-1] = index & (1<<(d-1));
202 result &= ~(1ull << (d-1));
203 result |= index & (1ull << (d-1));
204
205 index &= ~(1<<(d-1));
206 }
207 if (index >= subEnt<dim>(d-1,cc))
208 {
209 if ((index - subEnt<dim>(d-1,cc)) / subEnt<dim>(d-1,cc-1) == 1)
210 {
211 result |= 1ull << (d-1);
212 }
213 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
214 cc--;
215 }
216 }
217 return result;
218 }
219
220 };
221
229 template<int dim>
230 std::bitset<dim> entityMove(int index, int cc)
231 {
232 return EntityShiftTable<calculate_entity_move<dim>,dim>::evaluate(index,cc);
233 }
234
235#endif //DOXYGEN
236
237 } // namespace Yasp.
238
239 template<int codim, int dim, class GridImp>
240 class YaspEntity
241 : public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
242 {
243
244 template<int, PartitionIteratorType, typename>
245 friend class YaspLevelIterator;
246
247 public:
248 typedef typename GridImp::ctype ctype;
249
250 typedef typename GridImp::template Codim<codim>::Geometry Geometry;
251 typedef typename GridImp::Traits::template Codim<codim>::GeometryImpl GeometryImpl;
252
253 typedef typename GridImp::template Codim<codim>::EntitySeed EntitySeed;
254
256 int level () const
257 {
258 return _g->level();
259 }
260
264 EntitySeed seed() const
265 {
266 return EntitySeed(YaspEntitySeed<codim,GridImp>(_g->level(), _it.coord(), _it.which()));
267 }
268
270 Geometry geometry () const
271 {
272 GeometryImpl _geometry(_it.lowerleft(),_it.upperright(),_it.shift());
273 return Geometry(_geometry);
274 }
275
279 constexpr GeometryType type () const
280 {
281 return GeometryTypes::cube(Geometry::mydimension);
282 }
283
288 unsigned int subEntities (unsigned int cc) const
289 {
290 return Dune::Yasp::subEnt<dim>(dim-codim,cc-codim);
291 }
292
294 PartitionType partitionType () const
295 {
296 if (_g->interior[codim].inside(_it.coord(),_it.shift()))
297 return InteriorEntity;
298 if (_g->interiorborder[codim].inside(_it.coord(),_it.shift()))
299 return BorderEntity;
300 if (_g->overlap[codim].inside(_it.coord(),_it.shift()))
301 return OverlapEntity;
302 if (_g->overlapfront[codim].inside(_it.coord(),_it.shift()))
303 return FrontEntity;
304 return GhostEntity;
305 }
306
307 typedef typename GridImp::YGridLevelIterator YGLI;
308 typedef typename GridImp::YGrid::Iterator I;
309 YaspEntity ()
310 {}
311
312 YaspEntity (const YGLI& g, const I& it)
313 : _it(it), _g(g)
314 {}
315
316 YaspEntity (YGLI&& g, const I&& it)
317 : _it(std::move(it)), _g(std::move(g))
318 {}
319
321 bool equals (const YaspEntity& e) const
322 {
323 return _it == e._it && _g == e._g;
324 }
325
326 // IndexSets needs access to the private index methods
327 friend class Dune::YaspIndexSet<GridImp,true>;
328 friend class Dune::YaspIndexSet<GridImp,false>;
329 friend class Dune::YaspGlobalIdSet<GridImp>;
330 typedef typename GridImp::PersistentIndexType PersistentIndexType;
331
333 PersistentIndexType persistentIndex () const
334 {
335 // get size of global grid (in elements)
336 std::array<int,dim> size;
337
338 for (int i=0; i<dim; i++)
339 {
340 // correct size according to shift
341 size[i] = _g->mg->levelSize(_g->level(), i);
342 if (!_it.shift(i))
343 size[i]++;
344 }
345
346 // encode codim
347 PersistentIndexType id(_it.shift().to_ulong());
348
349 // encode level
350 id = id << yaspgrid_level_bits;
351 id = id+PersistentIndexType(_g->level());
352
353 // encode coordinates
354 for (int i=dim-1; i>=0; i--)
355 {
356 id = id << yaspgrid_dim_bits;
357 id = id+PersistentIndexType(_it.coord(i));
358 }
359
360 return id;
361 }
362
364 int compressedIndex () const
365 {
366 return _it.superindex();
367 }
368
370 int subCompressedIndex (int i, unsigned int cc) const
371 {
372 // get the shift of the entity and the subentity
373 // the subentity shift is only available in the space spanned by the entity
374 std::bitset<dim-codim> subent_shift = Dune::Yasp::entityShift<dim-codim>(i,cc-codim);
375 std::bitset<dim-codim> subent_move = Dune::Yasp::entityMove<dim-codim>(i,cc-codim);
376
377 std::bitset<dim> shift = _it.shift();
378 std::array<int, dim> coord = _it.coord();
379 for (int j=0, k=0; j<dim; j++)
380 {
381 if (!shift[j])
382 continue;
383
384 coord[j] += subent_move[k];
385 shift[j] = subent_shift[k];
386 k++;
387 }
388
389 int which = _g->overlapfront[cc].shiftmapping(shift);
390 return _g->overlapfront[cc].superindex(coord,which);
391 }
392 public:
393 const I& transformingsubiterator() const { return _it; }
394 const YGLI& gridlevel() const { return _g; }
395 I& transformingsubiterator() { return _it; }
396 YGLI& gridlevel() { return _g; }
397 const GridImp * yaspgrid() const { return _g->mg; }
398 protected:
399 I _it = {}; // position in the grid level
400 YGLI _g = {}; // access to grid level
401 };
402
403
404 // specialization for codim=0
405 template<int dim, class GridImp>
406 class YaspEntity<0,dim,GridImp>
407 : public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
408 {
409 constexpr static int dimworld = GridImp::dimensionworld;
410
411 typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
412
413 template<int, PartitionIteratorType, typename>
414 friend class YaspLevelIterator;
415
416 template<typename>
417 friend class YaspHierarchicIterator;
418
419 public:
420 typedef typename GridImp::ctype ctype;
421
422 typedef typename GridImp::YGridLevelIterator YGLI;
423 typedef typename GridImp::YGrid::Iterator I;
424
425 typedef typename GridImp::template Codim< 0 >::Geometry Geometry;
426 typedef typename GridImp::template Codim< 0 >::LocalGeometry LocalGeometry;
427
428 template <int cd>
429 struct Codim
430 {
431 typedef typename GridImp::template Codim<cd>::Entity Entity;
432 };
433
434 typedef typename GridImp::template Codim<0>::Entity Entity;
435 typedef typename GridImp::template Codim<0>::EntitySeed EntitySeed;
436 typedef typename GridImp::LevelIntersectionIterator IntersectionIterator;
437 typedef typename GridImp::LevelIntersectionIterator LevelIntersectionIterator;
438 typedef typename GridImp::LeafIntersectionIterator LeafIntersectionIterator;
439 typedef typename GridImp::HierarchicIterator HierarchicIterator;
440
442 typedef typename GridImp::PersistentIndexType PersistentIndexType;
443
445 typedef typename GridImp::YGrid::iTupel iTupel;
446
447 // constructor
448 YaspEntity ()
449 {}
450
451 YaspEntity (const YGLI& g, const I& it)
452 : _it(it), _g(g)
453 {}
454
455 YaspEntity (const YGLI& g, I&& it)
456 : _it(std::move(it)), _g(g)
457 {}
458
459 YaspEntity (YGLI&& g, I&& it)
460 : _it(std::move(it)), _g(std::move(g))
461 {}
462
464 bool equals (const YaspEntity& e) const
465 {
466 return _it == e._it && _g == e._g;
467 }
468
470 int level () const { return _g->level(); }
471
475 EntitySeed seed () const {
476 return EntitySeed(YaspEntitySeed<0,GridImp>(_g->level(), _it.coord()));
477 }
478
480 PartitionType partitionType () const
481 {
482 if (_g->interior[0].inside(_it.coord(),_it.shift()))
483 return InteriorEntity;
484 if (_g->overlap[0].inside(_it.coord(),_it.shift()))
485 return OverlapEntity;
486 DUNE_THROW(GridError, "Impossible GhostEntity");
487 return GhostEntity;
488 }
489
491 Geometry geometry () const {
492 // the element geometry
493 auto ll = _it.lowerleft();
494 auto ur = _it.upperright();
495
496 // If on periodic overlap, transform coordinates by domain size
497 for (int i=0; i<dimworld; i++) {
498 if (gridlevel()->mg->isPeriodic(i)) {
499 int coord = transformingsubiterator().coord(i);
500 if (coord < 0) {
501 auto size = _g->mg->domainSize()[i];
502 ll[i] += size;
503 ur[i] += size;
504 } else if (coord + 1 > gridlevel()->mg->levelSize(gridlevel()->level(),i)) {
505 auto size = _g->mg->domainSize()[i];
506 ll[i] -= size;
507 ur[i] -= size;
508 }
509 }
510 }
511
512 GeometryImpl _geometry(ll,ur);
513 return Geometry( _geometry );
514 }
515
519 constexpr GeometryType type () const
520 {
521 return GeometryTypes::cube(Geometry::mydimension);
522 }
523
528 template<int cc> int count () const
529 {
530 return Dune::Yasp::subEnt<dim>(dim,cc);
531 }
532
537 unsigned int subEntities (unsigned int codim) const
538 {
539 return Dune::Yasp::subEnt<dim>(dim,codim);
540 }
541
544 template<int cc>
545 typename Codim<cc>::Entity subEntity (int i) const
546 {
547 // calculate move bitset
548 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
549
550 // get the coordinate and modify it
551 iTupel coord = _it.coord();
552 for (int j=0; j<dim; j++)
553 if (move[j])
554 coord[j]++;
555
556 int which = _g->overlapfront[cc].shiftmapping(Dune::Yasp::entityShift<dim>(i,cc));
557 return typename Codim<cc>::Entity(YaspEntity<cc,GridImp::dimension,GridImp>(_g,_g->overlapfront[cc].begin(coord, which)));
558 }
559
561 Entity father () const
562 {
563 // check if coarse level exists
564 if (_g->level()<=0)
565 DUNE_THROW(GridError, "tried to call father on level 0");
566
567 // yes, get iterator to it
568 YGLI cg(_g);
569 --cg;
570
571 // coordinates of the cell
572 iTupel coord = _it.coord();
573
574 // get coordinates on next coarser level
575 for (int k=0; k<dim; k++) coord[k] = coord[k]/2;
576
577 return Entity(YaspEntity<0,GridImp::dimension,GridImp>(cg,cg->overlap[0].begin(coord)));
578 }
579
581 bool hasFather () const
582 {
583 return (_g->level()>0);
584 }
585
588 LocalGeometry geometryInFather () const
589 {
590 // configure one of the 2^dim transformations
591 FieldVector<ctype,dim> ll(0.0),ur(0.5);
592
593 for (int k=0; k<dim; k++)
594 {
595 if (_it.coord(k)%2)
596 {
597 ll[k] = 0.5;
598 ur[k] = 1.0;
599 }
600 }
601
602 return LocalGeometry( YaspGeometry<dim,dim,GridImp>(ll,ur) );
603 }
604
605 const I& transformingsubiterator () const { return _it; }
606 const YGLI& gridlevel () const { return _g; }
607 I& transformingsubiterator() { return _it; }
608 YGLI& gridlevel() { return _g; }
609 const GridImp* yaspgrid () const { return _g->mg; }
610
611 bool isLeaf() const
612 {
613 return (_g->level() == yaspgrid()->maxLevel());
614 }
615
618 bool isNew () const { return yaspgrid()->adaptRefCount > 0 && yaspgrid()->maxLevel() < _g->level() + yaspgrid()->adaptRefCount; }
619
622 bool mightVanish () const { return false; }
623
625 IntersectionIterator ibegin () const
626 {
627 return YaspIntersectionIterator<GridImp>(*this,false);
628 }
629
631 LeafIntersectionIterator ileafbegin () const
632 {
633 // only if entity is leaf this iterator delivers intersections
634 return YaspIntersectionIterator<GridImp>(*this, ! isLeaf() );
635 }
636
638 LevelIntersectionIterator ilevelbegin () const
639 {
640 return ibegin();
641 }
642
644 IntersectionIterator iend () const
645 {
646 return YaspIntersectionIterator<GridImp>(*this,true);
647 }
648
650 LeafIntersectionIterator ileafend () const
651 {
652 return iend();
653 }
654
656 LevelIntersectionIterator ilevelend () const
657 {
658 return iend();
659 }
660
665 HierarchicIterator hbegin (int maxlevel) const
666 {
667 return YaspHierarchicIterator<GridImp>(_g,_it,maxlevel);
668 }
669
671 HierarchicIterator hend (int /* maxlevel */) const
672 {
673 return YaspHierarchicIterator<GridImp>(_g,_it,_g->level());
674 }
675
676 private:
677 // IndexSets needs access to the private index methods
678 friend class Dune::YaspIndexSet<GridImp,true>;
679 friend class Dune::YaspIndexSet<GridImp,false>;
680 friend class Dune::YaspGlobalIdSet<GridImp>;
681
683 PersistentIndexType persistentIndex () const
684 {
685 // encode codim
686 PersistentIndexType id(_it.shift().to_ulong());
687
688 // encode level
689 id = id << yaspgrid_level_bits;
690 id = id+PersistentIndexType(_g->level());
691
692
693 // encode coordinates
694 for (int i=dim-1; i>=0; i--)
695 {
696 id = id << yaspgrid_dim_bits;
697 id = id+PersistentIndexType(_it.coord(i));
698 }
699
700 return id;
701 }
702
704 int compressedIndex () const
705 {
706 return _it.superindex();
707 }
708
710 PersistentIndexType subPersistentIndex (int i, int cc) const
711 {
712 // calculate shift and move bitsets
713 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
714 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
715
716 int trailing = (cc == dim) ? 1000 : 0;
717
718 std::array<int,dim> size = _g->mg->levelSize(_g->level());
719 std::array<int, dim> coord = _it.coord();
720 for (int j=0; j<dim; j++)
721 {
722 // correct size according to shift
723 if (!shift[j])
724 size[j]++;
725
726 // move the coordinates to the cell on which the entity lives
727 if (move[j])
728 coord[j]++;
729 }
730
731 for (int j=0; j<dim; j++)
732 {
733 // in the codim==dim case, count trailing zeroes.
734 if (cc == dim)
735 {
736 int zeroes = 0;
737 for (int k=0; k<_g->level(); k++)
738 if (coord[j] & (1<<k))
739 break;
740 else
741 zeroes++;
742 trailing = std::min(trailing,zeroes);
743 }
744 }
745
746 // encode codim
747 PersistentIndexType id(shift.to_ulong());
748
749 // encode level
750 id = id << yaspgrid_level_bits;
751 id = id+PersistentIndexType(_g->level()-trailing);
752
753 // encode coordinates
754 for (int j=dim-1; j>=0; j--)
755 {
756 id = id << yaspgrid_dim_bits;
757 id = id+PersistentIndexType(coord[j]>>trailing);
758 }
759
760 return id;
761 }
762
764 int subCompressedIndex (int i, int cc) const
765 {
766 // get shift and move of the subentity in question
767 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
768 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
769
770 std::array<int, dim> coord = _it.coord();
771 for (int j=0; j<dim; j++)
772 coord[j] += move[j];
773
774 int which = _g->overlapfront[cc].shiftmapping(shift);
775 return _g->overlapfront[cc].superindex(coord,which);
776 }
777
778 I _it = {}; // position in the grid level
779 YGLI _g = {}; // access to grid level
780 };
781
782
783 // specialization for codim=dim (vertex)
784 template<int dim, class GridImp>
785 class YaspEntity<dim,dim,GridImp>
786 : public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
787 {
788 constexpr static int dimworld = GridImp::dimensionworld;
789
790 template<int, PartitionIteratorType, typename>
791 friend class YaspLevelIterator;
792
793 typedef typename GridImp::Traits::template Codim<dim>::GeometryImpl GeometryImpl;
794
795 public:
796 typedef typename GridImp::ctype ctype;
797
798 typedef typename GridImp::YGridLevelIterator YGLI;
799 typedef typename GridImp::YGrid::Iterator I;
800
801 typedef typename GridImp::template Codim<dim>::Geometry Geometry;
802
803 typedef typename GridImp::template Codim<dim>::EntitySeed EntitySeed;
804
806 typedef typename GridImp::PersistentIndexType PersistentIndexType;
807
809 typedef typename GridImp::YGrid::iTupel iTupel;
810
811 // constructor
812 YaspEntity ()
813 {}
814
815 YaspEntity (const YGLI& g, const I& it)
816 : _it(it), _g(g)
817 {}
818
819 YaspEntity (YGLI&& g, I&& it)
820 : _it(std::move(it)), _g(std::move(g))
821 {}
822
824 bool equals (const YaspEntity& e) const
825 {
826 return _it == e._it && _g == e._g;
827 }
828
830 int level () const {return _g->level();}
831
835 EntitySeed seed () const {
836 return EntitySeed(YaspEntitySeed<dim,GridImp>(_g->level(), _it.coord(), _it.which()));
837 }
838
843 unsigned int subEntities (unsigned int cc) const
844 {
845 return Dune::Yasp::subEnt<dim>(dim-dim,cc-dim);
846 }
847
849 Geometry geometry () const {
850 GeometryImpl _geometry((_it).lowerleft());
851 return Geometry( _geometry );
852 }
853
857 constexpr GeometryType type () const
858 {
859 return GeometryTypes::cube(Geometry::mydimension);
860 }
861
863 PartitionType partitionType () const
864 {
865 if (_g->interior[dim].inside(_it.coord(),_it.shift()))
866 return InteriorEntity;
867 if (_g->interiorborder[dim].inside(_it.coord(),_it.shift()))
868 return BorderEntity;
869 if (_g->overlap[dim].inside(_it.coord(),_it.shift()))
870 return OverlapEntity;
871 if (_g->overlapfront[dim].inside(_it.coord(),_it.shift()))
872 return FrontEntity;
873 return GhostEntity;
874 }
875
877 int subCompressedIndex (int, unsigned int ) const
878 {
879 return compressedIndex();
880 }
881
882 private:
883 // IndexSets needs access to the private index methods
884 friend class Dune::YaspIndexSet<GridImp,true>;
885 friend class Dune::YaspIndexSet<GridImp,false>;
886 friend class Dune::YaspGlobalIdSet<GridImp>;
887
889 PersistentIndexType persistentIndex () const
890 {
891 // get coordinate and size of global grid
892 iTupel size = _g->mg->levelSize(_g->level());
893
894 for (int i=0; i<dim; i++)
895 {
896 // we have vertices, add 1 size to all directions
897 size[i]++;
898 }
899
900 // determine min number of trailing zeroes
901 int trailing = 1000;
902 for (int i=0; i<dim; i++)
903 {
904 // count trailing zeros
905 int zeros = 0;
906 for (int j=0; j<_g->level(); j++)
907 if (_it.coord(i)&(1<<j))
908 break;
909 else
910 zeros++;
911 trailing = std::min(trailing,zeros);
912 }
913
914 // determine the level of this vertex
915 int level = _g->level()-trailing;
916
917 // encode codim: shift vector of vertices is 0.
918 PersistentIndexType id(0);
919
920 // encode level
921 id = id << yaspgrid_level_bits;
922 id = id+PersistentIndexType(level);
923
924 // encode coordinates
925 for (int i=dim-1; i>=0; i--)
926 {
927 id = id << yaspgrid_dim_bits;
928 id = id+PersistentIndexType(_it.coord(i)>>trailing);
929 }
930
931 return id;
932 }
933
935 int compressedIndex () const { return _it.superindex();}
936
937 public:
938 const I& transformingsubiterator() const { return _it; }
939 const YGLI& gridlevel() const { return _g; }
940 I& transformingsubiterator() { return _it; }
941 YGLI& gridlevel() { return _g; }
942
943 const GridImp * yaspgrid() const { return _g->mg; }
944 protected:
945 I _it = {}; // position in the grid level
946 YGLI _g = {}; // access to grid level
947 };
948
949} // namespace Dune
950
951#endif // DUNE_GRID_YASPGRIDENTITY_HH
Dune::Geometry::mydimension
static constexpr int mydimension
geometry dimension
Definition: geometry.hh:94
Dune::YaspGlobalIdSet
persistent, globally unique Ids
Definition: yaspgrididset.hh:25
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:44
Dune::YaspIndexSet
Implementation of Level- and LeafIndexSets for YaspGrid.
Definition: yaspgridindexsets.hh:25
Dune::VTK::GeometryType
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
DUNE_THROW
#define DUNE_THROW(E,...)
Definition: exceptions.hh:314
Dune::PartitionType
PartitionType
Attributes used in the generic overlap model.
Definition: gridenums.hh:30
Dune::FrontEntity
@ FrontEntity
on boundary between overlap and ghost
Definition: gridenums.hh:34
Dune::InteriorEntity
@ InteriorEntity
all interior entities
Definition: gridenums.hh:31
Dune::GhostEntity
@ GhostEntity
ghost entities
Definition: gridenums.hh:35
Dune::BorderEntity
@ BorderEntity
on boundary between interior and overlap
Definition: gridenums.hh:32
Dune::OverlapEntity
@ OverlapEntity
all entities lying in the overlap zone
Definition: gridenums.hh:33
Dune::GeometryTypes::cube
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
math.hh
Some useful basic math stuff.
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
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:113
Dune::power
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
std
STL namespace.
type.hh
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Nov 3, 23:36, 2025)