Dune Core Modules (2.9.0)

yaspgridentity.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (C) 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 // the actual implementation
47 [[deprecated("Use binomial from dune-common's math.hh")]]
48 static constexpr int binomial(int d, int c)
49 {
50 long binomial=1;
51 for (int i=d-c+1; i<=d; i++)
52 binomial *= i;
53 for (long i=2; i<=c; i++)
54 binomial /= i;
55 return binomial;
56 }
57
58 private:
59 // prevent construction
60 BinomialTable() = delete;
61
62 // compute binomial(r, c) and advance row `r` and column `c`
63 static constexpr int nextValue(int& r, int& c)
64 {
65 const auto result = Dune::binomial(r, c);
66
67 c += 1;
68 if (c > r) {
69 r += 1;
70 c = 0;
71 }
72
73 return result;
74 }
75
76 template<std::size_t... I>
77 static constexpr std::array<int, sizeof...(I)> computeValues(std::index_sequence<I...>)
78 {
79 int r = 0, c = 0;
80 return {{ ((void)I, nextValue(r, c))... }};
81 }
82
83 template<std::size_t... I>
84 static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
85 { return {{ (I*(I+1)/2)... }}; }
86
87 static constexpr std::array<int,(n+1)*(n+2)/2> _values = computeValues(std::make_index_sequence<(n+1)*(n+2)/2>{});
88 static constexpr std::array<int,n+1> _offsets = computeOffsets(std::make_index_sequence<n+1>{});
89 };
90
91#if __cplusplus < 201703L
92 template<int n>
93 constexpr std::array<int,(n+1)*(n+2)/2> BinomialTable<n>::_values;
94 template<int n>
95 constexpr std::array<int,n+1> BinomialTable<n>::_offsets;
96#endif
97
104 template<int dimworld>
105 constexpr int subEnt(int d, int c)
106 {
107 return (d < c ? 0 : BinomialTable<dimworld>::evaluate(d,c) << c);
108 }
109
110 // Make a table mapping all subentities of a codim 0 entity to a value.
111 // F is the functor to be evaluated.
112 template<typename F, int dim>
113 struct EntityShiftTable
114 {
115 typedef std::bitset<dim> value_type;
116
117 static value_type evaluate(int i, int codim)
118 {
119 return {_values[_offsets[codim] + i]};
120 }
121
122 private:
123
124 // prevent construction
125 EntityShiftTable() = delete;
126
127 // compute offset of codimension `codim` entities and advance `offset`
128 static constexpr int nextOffset(int& offset, int codim)
129 {
130 if (codim == 0) {
131 offset = 0;
132 return 0;
133 }
134
135 offset += subEnt<dim>(dim, codim-1);
136 return offset;
137 }
138
139 template<std::size_t... I>
140 static constexpr std::array<int, sizeof...(I)> computeOffsets(std::index_sequence<I...>)
141 {
142 int offset = 0;
143 return {{ (nextOffset(offset, I))... }};
144 }
145
146 // compute shift table entry for (`codim`, `i`) and advance `codim`, `i`
147 static constexpr unsigned char nextValue(int& codim, int& i)
148 {
149 const auto result = F::evaluate(i, codim);
150
151 i += 1;
152 if (i >= subEnt<dim>(dim, codim)) {
153 codim += 1;
154 i = 0;
155 }
156
157 return result;
158 }
159
160 template<std::size_t... I>
161 static constexpr std::array<unsigned char, sizeof...(I)> computeValues(std::index_sequence<I...>)
162 {
163 int codim = 0, i = 0;
164 return {{ ((void)I, nextValue(codim, i))... }};
165 }
166
167 static constexpr std::array<int,dim+1> _offsets = computeOffsets(std::make_index_sequence<dim+1>{});
168 static constexpr std::array<unsigned char,Dune::power(3,dim)> _values = computeValues(std::make_index_sequence<Dune::power(3,dim)>{});
169
170 };
171
172#if __cplusplus < 201703L
173 template<typename F, int dim>
174 constexpr std::array<int,dim+1> EntityShiftTable<F, dim>::_offsets;
175 template<typename F, int dim>
176 constexpr std::array<unsigned char,Dune::power(3,dim)> EntityShiftTable<F, dim>::_values;
177#endif
178
179 // functor for doing the actual entity shift calculation
180 template<int dim>
181 struct calculate_entity_shift
182 {
183 static constexpr unsigned long long evaluate(int index, int cc)
184 {
185 auto result = 0ull;
186 for (int d = dim; d>0; d--)
187 {
188 if (cc == d)
189 return result;
190 if (index < subEnt<dim>(d-1,cc))
191 result |= 1ull << (d-1);
192 else
193 {
194 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
195 cc--;
196 }
197 }
198 return result;
199 }
200 };
201
210 template<int dim>
211 std::bitset<dim> entityShift(int index, int cc)
212 {
213 return EntityShiftTable<calculate_entity_shift<dim>,dim>::evaluate(index,cc);
214 }
215
216 // functor for doing the actual entity move calculation
217 template<int dim>
218 struct calculate_entity_move
219 {
220 static constexpr unsigned long long evaluate(int index, int cc)
221 {
222 auto result = 0ull;
223 for (int d = dim; d>0; d--)
224 {
225 if (d == cc)
226 {
227 // result[d-1] = index & (1<<(d-1));
228 result &= ~(1ull << (d-1));
229 result |= index & (1ull << (d-1));
230
231 index &= ~(1<<(d-1));
232 }
233 if (index >= subEnt<dim>(d-1,cc))
234 {
235 if ((index - subEnt<dim>(d-1,cc)) / subEnt<dim>(d-1,cc-1) == 1)
236 {
237 result |= 1ull << (d-1);
238 }
239 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
240 cc--;
241 }
242 }
243 return result;
244 }
245
246 };
247
255 template<int dim>
256 std::bitset<dim> entityMove(int index, int cc)
257 {
258 return EntityShiftTable<calculate_entity_move<dim>,dim>::evaluate(index,cc);
259 }
260
261#endif //DOXYGEN
262
263 } // namespace Yasp.
264
265 template<int codim, int dim, class GridImp>
266 class YaspEntity
267 : public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
268 {
269
270 template<int, PartitionIteratorType, typename>
271 friend class YaspLevelIterator;
272
273 public:
274 typedef typename GridImp::ctype ctype;
275
276 typedef typename GridImp::template Codim<codim>::Geometry Geometry;
277 typedef typename GridImp::Traits::template Codim<codim>::GeometryImpl GeometryImpl;
278
279 typedef typename GridImp::template Codim<codim>::EntitySeed EntitySeed;
280
282 int level () const
283 {
284 return _g->level();
285 }
286
290 EntitySeed seed() const
291 {
292 return EntitySeed(YaspEntitySeed<codim,GridImp>(_g->level(), _it.coord(), _it.which()));
293 }
294
296 Geometry geometry () const
297 {
298 GeometryImpl _geometry(_it.lowerleft(),_it.upperright(),_it.shift());
299 return Geometry(_geometry);
300 }
301
305 constexpr GeometryType type () const
306 {
307 return GeometryTypes::cube(Geometry::mydimension);
308 }
309
314 unsigned int subEntities (unsigned int cc) const
315 {
316 return Dune::Yasp::subEnt<dim>(dim-codim,cc-codim);
317 }
318
320 PartitionType partitionType () const
321 {
322 if (_g->interior[codim].inside(_it.coord(),_it.shift()))
323 return InteriorEntity;
324 if (_g->interiorborder[codim].inside(_it.coord(),_it.shift()))
325 return BorderEntity;
326 if (_g->overlap[codim].inside(_it.coord(),_it.shift()))
327 return OverlapEntity;
328 if (_g->overlapfront[codim].inside(_it.coord(),_it.shift()))
329 return FrontEntity;
330 return GhostEntity;
331 }
332
333 typedef typename GridImp::YGridLevelIterator YGLI;
334 typedef typename GridImp::YGrid::Iterator I;
335 YaspEntity ()
336 {}
337
338 YaspEntity (const YGLI& g, const I& it)
339 : _it(it), _g(g)
340 {}
341
342 YaspEntity (YGLI&& g, const I&& it)
343 : _it(std::move(it)), _g(std::move(g))
344 {}
345
347 bool equals (const YaspEntity& e) const
348 {
349 return _it == e._it && _g == e._g;
350 }
351
352 // IndexSets needs access to the private index methods
353 friend class Dune::YaspIndexSet<GridImp,true>;
354 friend class Dune::YaspIndexSet<GridImp,false>;
355 friend class Dune::YaspGlobalIdSet<GridImp>;
356 typedef typename GridImp::PersistentIndexType PersistentIndexType;
357
359 PersistentIndexType persistentIndex () const
360 {
361 // get size of global grid (in elements)
362 std::array<int,dim> size;
363
364 for (int i=0; i<dim; i++)
365 {
366 // correct size according to shift
367 size[i] = _g->mg->levelSize(_g->level(), i);
368 if (!_it.shift(i))
369 size[i]++;
370 }
371
372 // encode codim
373 PersistentIndexType id(_it.shift().to_ulong());
374
375 // encode level
376 id = id << yaspgrid_level_bits;
377 id = id+PersistentIndexType(_g->level());
378
379 // encode coordinates
380 for (int i=dim-1; i>=0; i--)
381 {
382 id = id << yaspgrid_dim_bits;
383 id = id+PersistentIndexType(_it.coord(i));
384 }
385
386 return id;
387 }
388
390 int compressedIndex () const
391 {
392 return _it.superindex();
393 }
394
396 int subCompressedIndex (int i, unsigned int cc) const
397 {
398 // get the shift of the entity and the subentity
399 // the subentity shift is only available in the space spanned by the entity
400 std::bitset<dim-codim> subent_shift = Dune::Yasp::entityShift<dim-codim>(i,cc-codim);
401 std::bitset<dim-codim> subent_move = Dune::Yasp::entityMove<dim-codim>(i,cc-codim);
402
403 std::bitset<dim> shift = _it.shift();
404 std::array<int, dim> coord = _it.coord();
405 for (int j=0, k=0; j<dim; j++)
406 {
407 if (!shift[j])
408 continue;
409
410 coord[j] += subent_move[k];
411 shift[j] = subent_shift[k];
412 k++;
413 }
414
415 int which = _g->overlapfront[cc].shiftmapping(shift);
416 return _g->overlapfront[cc].superindex(coord,which);
417 }
418 public:
419 const I& transformingsubiterator() const { return _it; }
420 const YGLI& gridlevel() const { return _g; }
421 I& transformingsubiterator() { return _it; }
422 YGLI& gridlevel() { return _g; }
423 const GridImp * yaspgrid() const { return _g->mg; }
424 protected:
425 I _it = {}; // position in the grid level
426 YGLI _g = {}; // access to grid level
427 };
428
429
430 // specialization for codim=0
431 template<int dim, class GridImp>
432 class YaspEntity<0,dim,GridImp>
433 : public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
434 {
435 constexpr static int dimworld = GridImp::dimensionworld;
436
437 typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
438
439 template<int, PartitionIteratorType, typename>
440 friend class YaspLevelIterator;
441
442 template<typename>
443 friend class YaspHierarchicIterator;
444
445 public:
446 typedef typename GridImp::ctype ctype;
447
448 typedef typename GridImp::YGridLevelIterator YGLI;
449 typedef typename GridImp::YGrid::Iterator I;
450
451 typedef typename GridImp::template Codim< 0 >::Geometry Geometry;
452 typedef typename GridImp::template Codim< 0 >::LocalGeometry LocalGeometry;
453
454 template <int cd>
455 struct Codim
456 {
457 typedef typename GridImp::template Codim<cd>::Entity Entity;
458 };
459
460 typedef typename GridImp::template Codim<0>::Entity Entity;
461 typedef typename GridImp::template Codim<0>::EntitySeed EntitySeed;
462 typedef typename GridImp::LevelIntersectionIterator IntersectionIterator;
463 typedef typename GridImp::LevelIntersectionIterator LevelIntersectionIterator;
464 typedef typename GridImp::LeafIntersectionIterator LeafIntersectionIterator;
465 typedef typename GridImp::HierarchicIterator HierarchicIterator;
466
468 typedef typename GridImp::PersistentIndexType PersistentIndexType;
469
471 typedef typename GridImp::YGrid::iTupel iTupel;
472
473 // constructor
474 YaspEntity ()
475 {}
476
477 YaspEntity (const YGLI& g, const I& it)
478 : _it(it), _g(g)
479 {}
480
481 YaspEntity (const YGLI& g, I&& it)
482 : _it(std::move(it)), _g(g)
483 {}
484
485 YaspEntity (YGLI&& g, I&& it)
486 : _it(std::move(it)), _g(std::move(g))
487 {}
488
490 bool equals (const YaspEntity& e) const
491 {
492 return _it == e._it && _g == e._g;
493 }
494
496 int level () const { return _g->level(); }
497
501 EntitySeed seed () const {
502 return EntitySeed(YaspEntitySeed<0,GridImp>(_g->level(), _it.coord()));
503 }
504
506 PartitionType partitionType () const
507 {
508 if (_g->interior[0].inside(_it.coord(),_it.shift()))
509 return InteriorEntity;
510 if (_g->overlap[0].inside(_it.coord(),_it.shift()))
511 return OverlapEntity;
512 DUNE_THROW(GridError, "Impossible GhostEntity");
513 return GhostEntity;
514 }
515
517 Geometry geometry () const {
518 // the element geometry
519 auto ll = _it.lowerleft();
520 auto ur = _it.upperright();
521
522 // If on periodic overlap, transform coordinates by domain size
523 for (int i=0; i<dimworld; i++) {
524 if (gridlevel()->mg->isPeriodic(i)) {
525 int coord = transformingsubiterator().coord(i);
526 if (coord < 0) {
527 auto size = _g->mg->domainSize()[i];
528 ll[i] += size;
529 ur[i] += size;
530 } else if (coord + 1 > gridlevel()->mg->levelSize(gridlevel()->level(),i)) {
531 auto size = _g->mg->domainSize()[i];
532 ll[i] -= size;
533 ur[i] -= size;
534 }
535 }
536 }
537
538 GeometryImpl _geometry(ll,ur);
539 return Geometry( _geometry );
540 }
541
545 constexpr GeometryType type () const
546 {
547 return GeometryTypes::cube(Geometry::mydimension);
548 }
549
554 template<int cc> int count () const
555 {
556 return Dune::Yasp::subEnt<dim>(dim,cc);
557 }
558
563 unsigned int subEntities (unsigned int codim) const
564 {
565 return Dune::Yasp::subEnt<dim>(dim,codim);
566 }
567
570 template<int cc>
571 typename Codim<cc>::Entity subEntity (int i) const
572 {
573 // calculate move bitset
574 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
575
576 // get the coordinate and modify it
577 iTupel coord = _it.coord();
578 for (int j=0; j<dim; j++)
579 if (move[j])
580 coord[j]++;
581
582 int which = _g->overlapfront[cc].shiftmapping(Dune::Yasp::entityShift<dim>(i,cc));
583 return typename Codim<cc>::Entity(YaspEntity<cc,GridImp::dimension,GridImp>(_g,_g->overlapfront[cc].begin(coord, which)));
584 }
585
587 Entity father () const
588 {
589 // check if coarse level exists
590 if (_g->level()<=0)
591 DUNE_THROW(GridError, "tried to call father on level 0");
592
593 // yes, get iterator to it
594 YGLI cg(_g);
595 --cg;
596
597 // coordinates of the cell
598 iTupel coord = _it.coord();
599
600 // get coordinates on next coarser level
601 for (int k=0; k<dim; k++) coord[k] = coord[k]/2;
602
603 return Entity(YaspEntity<0,GridImp::dimension,GridImp>(cg,cg->overlap[0].begin(coord)));
604 }
605
607 bool hasFather () const
608 {
609 return (_g->level()>0);
610 }
611
614 LocalGeometry geometryInFather () const
615 {
616 // configure one of the 2^dim transformations
617 FieldVector<ctype,dim> ll(0.0),ur(0.5);
618
619 for (int k=0; k<dim; k++)
620 {
621 if (_it.coord(k)%2)
622 {
623 ll[k] = 0.5;
624 ur[k] = 1.0;
625 }
626 }
627
628 return LocalGeometry( YaspGeometry<dim,dim,GridImp>(ll,ur) );
629 }
630
631 const I& transformingsubiterator () const { return _it; }
632 const YGLI& gridlevel () const { return _g; }
633 I& transformingsubiterator() { return _it; }
634 YGLI& gridlevel() { return _g; }
635 const GridImp* yaspgrid () const { return _g->mg; }
636
637 bool isLeaf() const
638 {
639 return (_g->level() == yaspgrid()->maxLevel());
640 }
641
644 bool isNew () const { return yaspgrid()->adaptRefCount > 0 && yaspgrid()->maxLevel() < _g->level() + yaspgrid()->adaptRefCount; }
645
648 bool mightVanish () const { return false; }
649
651 IntersectionIterator ibegin () const
652 {
653 return YaspIntersectionIterator<GridImp>(*this,false);
654 }
655
657 LeafIntersectionIterator ileafbegin () const
658 {
659 // only if entity is leaf this iterator delivers intersections
660 return YaspIntersectionIterator<GridImp>(*this, ! isLeaf() );
661 }
662
664 LevelIntersectionIterator ilevelbegin () const
665 {
666 return ibegin();
667 }
668
670 IntersectionIterator iend () const
671 {
672 return YaspIntersectionIterator<GridImp>(*this,true);
673 }
674
676 LeafIntersectionIterator ileafend () const
677 {
678 return iend();
679 }
680
682 LevelIntersectionIterator ilevelend () const
683 {
684 return iend();
685 }
686
691 HierarchicIterator hbegin (int maxlevel) const
692 {
693 return YaspHierarchicIterator<GridImp>(_g,_it,maxlevel);
694 }
695
697 HierarchicIterator hend (int /* maxlevel */) const
698 {
699 return YaspHierarchicIterator<GridImp>(_g,_it,_g->level());
700 }
701
702 private:
703 // IndexSets needs access to the private index methods
704 friend class Dune::YaspIndexSet<GridImp,true>;
705 friend class Dune::YaspIndexSet<GridImp,false>;
706 friend class Dune::YaspGlobalIdSet<GridImp>;
707
709 PersistentIndexType persistentIndex () const
710 {
711 // encode codim
712 PersistentIndexType id(_it.shift().to_ulong());
713
714 // encode level
715 id = id << yaspgrid_level_bits;
716 id = id+PersistentIndexType(_g->level());
717
718
719 // encode coordinates
720 for (int i=dim-1; i>=0; i--)
721 {
722 id = id << yaspgrid_dim_bits;
723 id = id+PersistentIndexType(_it.coord(i));
724 }
725
726 return id;
727 }
728
730 int compressedIndex () const
731 {
732 return _it.superindex();
733 }
734
736 PersistentIndexType subPersistentIndex (int i, int cc) const
737 {
738 // calculate shift and move bitsets
739 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
740 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
741
742 int trailing = (cc == dim) ? 1000 : 0;
743
744 std::array<int,dim> size = _g->mg->levelSize(_g->level());
745 std::array<int, dim> coord = _it.coord();
746 for (int j=0; j<dim; j++)
747 {
748 // correct size according to shift
749 if (!shift[j])
750 size[j]++;
751
752 // move the coordinates to the cell on which the entity lives
753 if (move[j])
754 coord[j]++;
755 }
756
757 for (int j=0; j<dim; j++)
758 {
759 // in the codim==dim case, count trailing zeroes.
760 if (cc == dim)
761 {
762 int zeroes = 0;
763 for (int k=0; k<_g->level(); k++)
764 if (coord[j] & (1<<k))
765 break;
766 else
767 zeroes++;
768 trailing = std::min(trailing,zeroes);
769 }
770 }
771
772 // encode codim
773 PersistentIndexType id(shift.to_ulong());
774
775 // encode level
776 id = id << yaspgrid_level_bits;
777 id = id+PersistentIndexType(_g->level()-trailing);
778
779 // encode coordinates
780 for (int j=dim-1; j>=0; j--)
781 {
782 id = id << yaspgrid_dim_bits;
783 id = id+PersistentIndexType(coord[j]>>trailing);
784 }
785
786 return id;
787 }
788
790 int subCompressedIndex (int i, int cc) const
791 {
792 // get shift and move of the subentity in question
793 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
794 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
795
796 std::array<int, dim> coord = _it.coord();
797 for (int j=0; j<dim; j++)
798 coord[j] += move[j];
799
800 int which = _g->overlapfront[cc].shiftmapping(shift);
801 return _g->overlapfront[cc].superindex(coord,which);
802 }
803
804 I _it = {}; // position in the grid level
805 YGLI _g = {}; // access to grid level
806 };
807
808
809 // specialization for codim=dim (vertex)
810 template<int dim, class GridImp>
811 class YaspEntity<dim,dim,GridImp>
812 : public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
813 {
814 constexpr static int dimworld = GridImp::dimensionworld;
815
816 template<int, PartitionIteratorType, typename>
817 friend class YaspLevelIterator;
818
819 typedef typename GridImp::Traits::template Codim<dim>::GeometryImpl GeometryImpl;
820
821 public:
822 typedef typename GridImp::ctype ctype;
823
824 typedef typename GridImp::YGridLevelIterator YGLI;
825 typedef typename GridImp::YGrid::Iterator I;
826
827 typedef typename GridImp::template Codim<dim>::Geometry Geometry;
828
829 typedef typename GridImp::template Codim<dim>::EntitySeed EntitySeed;
830
832 typedef typename GridImp::PersistentIndexType PersistentIndexType;
833
835 typedef typename GridImp::YGrid::iTupel iTupel;
836
837 // constructor
838 YaspEntity ()
839 {}
840
841 YaspEntity (const YGLI& g, const I& it)
842 : _it(it), _g(g)
843 {}
844
845 YaspEntity (YGLI&& g, I&& it)
846 : _it(std::move(it)), _g(std::move(g))
847 {}
848
850 bool equals (const YaspEntity& e) const
851 {
852 return _it == e._it && _g == e._g;
853 }
854
856 int level () const {return _g->level();}
857
861 EntitySeed seed () const {
862 return EntitySeed(YaspEntitySeed<dim,GridImp>(_g->level(), _it.coord(), _it.which()));
863 }
864
869 unsigned int subEntities (unsigned int cc) const
870 {
871 return Dune::Yasp::subEnt<dim>(dim-dim,cc-dim);
872 }
873
875 Geometry geometry () const {
876 GeometryImpl _geometry((_it).lowerleft());
877 return Geometry( _geometry );
878 }
879
883 constexpr GeometryType type () const
884 {
885 return GeometryTypes::cube(Geometry::mydimension);
886 }
887
889 PartitionType partitionType () const
890 {
891 if (_g->interior[dim].inside(_it.coord(),_it.shift()))
892 return InteriorEntity;
893 if (_g->interiorborder[dim].inside(_it.coord(),_it.shift()))
894 return BorderEntity;
895 if (_g->overlap[dim].inside(_it.coord(),_it.shift()))
896 return OverlapEntity;
897 if (_g->overlapfront[dim].inside(_it.coord(),_it.shift()))
898 return FrontEntity;
899 return GhostEntity;
900 }
901
903 int subCompressedIndex (int, unsigned int ) const
904 {
905 return compressedIndex();
906 }
907
908 private:
909 // IndexSets needs access to the private index methods
910 friend class Dune::YaspIndexSet<GridImp,true>;
911 friend class Dune::YaspIndexSet<GridImp,false>;
912 friend class Dune::YaspGlobalIdSet<GridImp>;
913
915 PersistentIndexType persistentIndex () const
916 {
917 // get coordinate and size of global grid
918 iTupel size = _g->mg->levelSize(_g->level());
919
920 for (int i=0; i<dim; i++)
921 {
922 // we have vertices, add 1 size to all directions
923 size[i]++;
924 }
925
926 // determine min number of trailing zeroes
927 int trailing = 1000;
928 for (int i=0; i<dim; i++)
929 {
930 // count trailing zeros
931 int zeros = 0;
932 for (int j=0; j<_g->level(); j++)
933 if (_it.coord(i)&(1<<j))
934 break;
935 else
936 zeros++;
937 trailing = std::min(trailing,zeros);
938 }
939
940 // determine the level of this vertex
941 int level = _g->level()-trailing;
942
943 // encode codim: shift vector of vertices is 0.
944 PersistentIndexType id(0);
945
946 // encode level
947 id = id << yaspgrid_level_bits;
948 id = id+PersistentIndexType(level);
949
950 // encode coordinates
951 for (int i=dim-1; i>=0; i--)
952 {
953 id = id << yaspgrid_dim_bits;
954 id = id+PersistentIndexType(_it.coord(i)>>trailing);
955 }
956
957 return id;
958 }
959
961 int compressedIndex () const { return _it.superindex();}
962
963 public:
964 const I& transformingsubiterator() const { return _it; }
965 const YGLI& gridlevel() const { return _g; }
966 I& transformingsubiterator() { return _it; }
967 YGLI& gridlevel() { return _g; }
968
969 const GridImp * yaspgrid() const { return _g->mg; }
970 protected:
971 I _it = {}; // position in the grid level
972 YGLI _g = {}; // access to grid level
973 };
974
975} // namespace Dune
976
977#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, m)
Definition: exceptions.hh:218
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:472
Dune::Hybrid::equals
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:402
math.hh
Some useful basic math stuff.
Dune
Dune namespace.
Definition: alignedallocator.hh:13
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:131
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  |  generated with Hugo v0.111.3 (Dec 21, 23:30, 2024)