Dune Core Modules (2.6.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 
24 namespace 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
persistent, globally unique Ids
Definition: yaspgrididset.hh:23
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
Implementation of Level- and LeafIndexSets for YaspGrid.
Definition: yaspgridindexsets.hh:23
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
PartitionType
Attributes used in the generic overlap model.
Definition: gridenums.hh:28
@ FrontEntity
on boundary between overlap and ghost
Definition: gridenums.hh:32
@ InteriorEntity
all interior entities
Definition: gridenums.hh:29
@ GhostEntity
ghost entities
Definition: gridenums.hh:33
@ BorderEntity
on boundary between interior and overlap
Definition: gridenums.hh:30
@ OverlapEntity
all entities lying in the overlap zone
Definition: gridenums.hh:31
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:435
int binomial(int upper, int lower)
calculate
Definition: simplex.cc:294
Dune namespace.
Definition: alignedallocator.hh:10
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 6, 22:30, 2024)