5#ifndef DUNE_GRID_YASPGRIDENTITY_HH
6#define DUNE_GRID_YASPGRIDENTITY_HH
41 static constexpr int evaluate(
int d,
int c)
43 return _values[_offsets[d] + c];
48 BinomialTable() =
delete;
51 static constexpr int nextValue(
int& r,
int& c)
64 template<std::size_t... I>
65 static constexpr std::array<int,
sizeof...(I)> computeValues(std::index_sequence<I...>)
68 return {{ ((void)I, nextValue(r, c))... }};
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)... }}; }
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>{});
85 template<
int dimworld>
86 constexpr int subEnt(
int d,
int c)
88 return (d < c ? 0 : BinomialTable<dimworld>::evaluate(d,c) << c);
93 template<
typename F,
int dim>
94 struct EntityShiftTable
96 typedef std::bitset<dim> value_type;
98 static value_type evaluate(
int i,
int codim)
100 return {_values[_offsets[codim] + i]};
106 EntityShiftTable() =
delete;
109 static constexpr int nextOffset(
int& offset,
int codim)
116 offset += subEnt<dim>(dim, codim-1);
120 template<std::size_t... I>
121 static constexpr std::array<int,
sizeof...(I)> computeOffsets(std::index_sequence<I...>)
124 return {{ (nextOffset(offset, I))... }};
128 static constexpr unsigned char nextValue(
int& codim,
int& i)
130 const auto result = F::evaluate(i, codim);
133 if (i >= subEnt<dim>(dim, codim)) {
141 template<std::size_t... I>
142 static constexpr std::array<
unsigned char,
sizeof...(I)> computeValues(std::index_sequence<I...>)
144 int codim = 0, i = 0;
145 return {{ ((void)I, nextValue(codim, i))... }};
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)>{});
155 struct calculate_entity_shift
157 static constexpr unsigned long long evaluate(
int index,
int cc)
160 for (
int d = dim; d>0; d--)
164 if (index < subEnt<dim>(d-1,cc))
165 result |= 1ull << (d-1);
168 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
185 std::bitset<dim> entityShift(
int index,
int cc)
187 return EntityShiftTable<calculate_entity_shift<dim>,dim>::evaluate(index,cc);
192 struct calculate_entity_move
194 static constexpr unsigned long long evaluate(
int index,
int cc)
197 for (
int d = dim; d>0; d--)
202 result &= ~(1ull << (d-1));
203 result |= index & (1ull << (d-1));
205 index &= ~(1<<(d-1));
207 if (index >= subEnt<dim>(d-1,cc))
209 if ((index - subEnt<dim>(d-1,cc)) / subEnt<dim>(d-1,cc-1) == 1)
211 result |= 1ull << (d-1);
213 index = (index - subEnt<dim>(d-1, cc)) % subEnt<dim>(d-1,cc-1);
230 std::bitset<dim> entityMove(
int index,
int cc)
232 return EntityShiftTable<calculate_entity_move<dim>,dim>::evaluate(index,cc);
239 template<
int codim,
int dim,
class Gr
idImp>
241 :
public EntityDefaultImplementation <codim,dim,GridImp,YaspEntity>
244 template<
int, PartitionIteratorType,
typename>
245 friend class YaspLevelIterator;
248 typedef typename GridImp::ctype ctype;
250 typedef typename GridImp::template Codim<codim>::Geometry Geometry;
251 typedef typename GridImp::Traits::template Codim<codim>::GeometryImpl GeometryImpl;
253 typedef typename GridImp::template Codim<codim>::EntitySeed EntitySeed;
264 EntitySeed seed()
const
266 return EntitySeed(YaspEntitySeed<codim,GridImp>(_g->level(), _it.coord(), _it.which()));
270 Geometry geometry ()
const
272 GeometryImpl _geometry(_it.lowerleft(),_it.upperright(),_it.shift());
273 return Geometry(_geometry);
288 unsigned int subEntities (
unsigned int cc)
const
290 return Dune::Yasp::subEnt<dim>(dim-codim,cc-codim);
296 if (_g->interior[codim].inside(_it.coord(),_it.shift()))
298 if (_g->interiorborder[codim].inside(_it.coord(),_it.shift()))
300 if (_g->overlap[codim].inside(_it.coord(),_it.shift()))
302 if (_g->overlapfront[codim].inside(_it.coord(),_it.shift()))
307 typedef typename GridImp::YGridLevelIterator YGLI;
308 typedef typename GridImp::YGrid::Iterator I;
312 YaspEntity (
const YGLI& g,
const I& it)
316 YaspEntity (YGLI&& g,
const I&& it)
317 : _it(
std::move(it)), _g(
std::move(g))
321 bool equals (
const YaspEntity& e)
const
323 return _it == e._it && _g == e._g;
330 typedef typename GridImp::PersistentIndexType PersistentIndexType;
333 PersistentIndexType persistentIndex ()
const
336 std::array<int,dim>
size;
338 for (
int i=0; i<dim; i++)
341 size[i] = _g->mg->levelSize(_g->level(), i);
347 PersistentIndexType
id(_it.shift().to_ulong());
350 id =
id << yaspgrid_level_bits;
351 id =
id+PersistentIndexType(_g->level());
354 for (
int i=dim-1; i>=0; i--)
356 id =
id << yaspgrid_dim_bits;
357 id =
id+PersistentIndexType(_it.coord(i));
364 int compressedIndex ()
const
366 return _it.superindex();
370 int subCompressedIndex (
int i,
unsigned int cc)
const
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);
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++)
384 coord[j] += subent_move[k];
385 shift[j] = subent_shift[k];
389 int which = _g->overlapfront[cc].shiftmapping(shift);
390 return _g->overlapfront[cc].superindex(coord,which);
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; }
405 template<
int dim,
class Gr
idImp>
406 class YaspEntity<0,dim,GridImp>
407 :
public EntityDefaultImplementation <0,dim,GridImp,YaspEntity>
409 constexpr static int dimworld = GridImp::dimensionworld;
411 typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
413 template<
int, PartitionIteratorType,
typename>
414 friend class YaspLevelIterator;
417 friend class YaspHierarchicIterator;
420 typedef typename GridImp::ctype ctype;
422 typedef typename GridImp::YGridLevelIterator YGLI;
423 typedef typename GridImp::YGrid::Iterator I;
425 typedef typename GridImp::template Codim< 0 >::Geometry Geometry;
426 typedef typename GridImp::template Codim< 0 >::LocalGeometry LocalGeometry;
431 typedef typename GridImp::template Codim<cd>::Entity Entity;
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;
442 typedef typename GridImp::PersistentIndexType PersistentIndexType;
445 typedef typename GridImp::YGrid::iTupel iTupel;
451 YaspEntity (
const YGLI& g,
const I& it)
455 YaspEntity (
const YGLI& g, I&& it)
456 : _it(
std::move(it)), _g(g)
459 YaspEntity (YGLI&& g, I&& it)
460 : _it(
std::move(it)), _g(
std::move(g))
464 bool equals (
const YaspEntity& e)
const
466 return _it == e._it && _g == e._g;
470 int level ()
const {
return _g->level(); }
475 EntitySeed seed ()
const {
476 return EntitySeed(YaspEntitySeed<0,GridImp>(_g->level(), _it.coord()));
482 if (_g->interior[0].inside(_it.coord(),_it.shift()))
484 if (_g->overlap[0].inside(_it.coord(),_it.shift()))
486 DUNE_THROW(GridError,
"Impossible GhostEntity");
491 Geometry geometry ()
const {
493 auto ll = _it.lowerleft();
494 auto ur = _it.upperright();
497 for (
int i=0; i<dimworld; i++) {
498 if (gridlevel()->mg->isPeriodic(i)) {
499 int coord = transformingsubiterator().coord(i);
501 auto size = _g->mg->domainSize()[i];
504 }
else if (coord + 1 > gridlevel()->mg->levelSize(gridlevel()->level(),i)) {
505 auto size = _g->mg->domainSize()[i];
512 GeometryImpl _geometry(ll,ur);
513 return Geometry( _geometry );
528 template<
int cc>
int count ()
const
530 return Dune::Yasp::subEnt<dim>(dim,cc);
537 unsigned int subEntities (
unsigned int codim)
const
539 return Dune::Yasp::subEnt<dim>(dim,codim);
545 typename Codim<cc>::Entity subEntity (
int i)
const
548 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
551 iTupel coord = _it.coord();
552 for (
int j=0; j<dim; j++)
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)));
561 Entity father ()
const
565 DUNE_THROW(GridError,
"tried to call father on level 0");
572 iTupel coord = _it.coord();
575 for (
int k=0; k<dim; k++) coord[k] = coord[k]/2;
577 return Entity(YaspEntity<0,GridImp::dimension,GridImp>(cg,cg->overlap[0].begin(coord)));
581 bool hasFather ()
const
583 return (_g->level()>0);
588 LocalGeometry geometryInFather ()
const
591 FieldVector<ctype,dim> ll(0.0),ur(0.5);
593 for (
int k=0; k<dim; k++)
602 return LocalGeometry( YaspGeometry<dim,dim,GridImp>(ll,ur) );
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; }
613 return (_g->level() == yaspgrid()->maxLevel());
618 bool isNew ()
const {
return yaspgrid()->adaptRefCount > 0 && yaspgrid()->maxLevel() < _g->level() + yaspgrid()->adaptRefCount; }
622 bool mightVanish ()
const {
return false; }
625 IntersectionIterator ibegin ()
const
627 return YaspIntersectionIterator<GridImp>(*
this,
false);
631 LeafIntersectionIterator ileafbegin ()
const
634 return YaspIntersectionIterator<GridImp>(*
this, ! isLeaf() );
638 LevelIntersectionIterator ilevelbegin ()
const
644 IntersectionIterator iend ()
const
646 return YaspIntersectionIterator<GridImp>(*
this,
true);
650 LeafIntersectionIterator ileafend ()
const
656 LevelIntersectionIterator ilevelend ()
const
665 HierarchicIterator hbegin (
int maxlevel)
const
667 return YaspHierarchicIterator<GridImp>(_g,_it,maxlevel);
671 HierarchicIterator hend (
int )
const
673 return YaspHierarchicIterator<GridImp>(_g,_it,_g->level());
683 PersistentIndexType persistentIndex ()
const
686 PersistentIndexType
id(_it.shift().to_ulong());
689 id =
id << yaspgrid_level_bits;
690 id =
id+PersistentIndexType(_g->level());
694 for (
int i=dim-1; i>=0; i--)
696 id =
id << yaspgrid_dim_bits;
697 id =
id+PersistentIndexType(_it.coord(i));
704 int compressedIndex ()
const
706 return _it.superindex();
710 PersistentIndexType subPersistentIndex (
int i,
int cc)
const
713 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
714 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
716 int trailing = (cc == dim) ? 1000 : 0;
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++)
731 for (
int j=0; j<dim; j++)
737 for (
int k=0; k<_g->level(); k++)
738 if (coord[j] & (1<<k))
742 trailing =
std::min(trailing,zeroes);
747 PersistentIndexType
id(shift.to_ulong());
750 id =
id << yaspgrid_level_bits;
751 id =
id+PersistentIndexType(_g->level()-trailing);
754 for (
int j=dim-1; j>=0; j--)
756 id =
id << yaspgrid_dim_bits;
757 id =
id+PersistentIndexType(coord[j]>>trailing);
764 int subCompressedIndex (
int i,
int cc)
const
767 std::bitset<dim> shift = Dune::Yasp::entityShift<dim>(i,cc);
768 std::bitset<dim> move = Dune::Yasp::entityMove<dim>(i,cc);
770 std::array<int, dim> coord = _it.coord();
771 for (
int j=0; j<dim; j++)
774 int which = _g->overlapfront[cc].shiftmapping(shift);
775 return _g->overlapfront[cc].superindex(coord,which);
784 template<
int dim,
class Gr
idImp>
785 class YaspEntity<dim,dim,GridImp>
786 :
public EntityDefaultImplementation <dim,dim,GridImp,YaspEntity>
788 constexpr static int dimworld = GridImp::dimensionworld;
790 template<
int, PartitionIteratorType,
typename>
791 friend class YaspLevelIterator;
793 typedef typename GridImp::Traits::template Codim<dim>::GeometryImpl GeometryImpl;
796 typedef typename GridImp::ctype ctype;
798 typedef typename GridImp::YGridLevelIterator YGLI;
799 typedef typename GridImp::YGrid::Iterator I;
801 typedef typename GridImp::template Codim<dim>::Geometry Geometry;
803 typedef typename GridImp::template Codim<dim>::EntitySeed EntitySeed;
806 typedef typename GridImp::PersistentIndexType PersistentIndexType;
809 typedef typename GridImp::YGrid::iTupel iTupel;
815 YaspEntity (
const YGLI& g,
const I& it)
819 YaspEntity (YGLI&& g, I&& it)
820 : _it(
std::move(it)), _g(
std::move(g))
824 bool equals (
const YaspEntity& e)
const
826 return _it == e._it && _g == e._g;
830 int level ()
const {
return _g->level();}
835 EntitySeed seed ()
const {
836 return EntitySeed(YaspEntitySeed<dim,GridImp>(_g->level(), _it.coord(), _it.which()));
843 unsigned int subEntities (
unsigned int cc)
const
845 return Dune::Yasp::subEnt<dim>(dim-dim,cc-dim);
849 Geometry geometry ()
const {
850 GeometryImpl _geometry((_it).lowerleft());
851 return Geometry( _geometry );
865 if (_g->interior[dim].inside(_it.coord(),_it.shift()))
867 if (_g->interiorborder[dim].inside(_it.coord(),_it.shift()))
869 if (_g->overlap[dim].inside(_it.coord(),_it.shift()))
871 if (_g->overlapfront[dim].inside(_it.coord(),_it.shift()))
877 int subCompressedIndex (
int,
unsigned int )
const
879 return compressedIndex();
889 PersistentIndexType persistentIndex ()
const
892 iTupel
size = _g->mg->levelSize(_g->level());
894 for (
int i=0; i<dim; i++)
902 for (
int i=0; i<dim; i++)
906 for (
int j=0; j<_g->level(); j++)
907 if (_it.coord(i)&(1<<j))
911 trailing =
std::min(trailing,zeros);
915 int level = _g->level()-trailing;
918 PersistentIndexType
id(0);
921 id =
id << yaspgrid_level_bits;
922 id =
id+PersistentIndexType(level);
925 for (
int i=dim-1; i>=0; i--)
927 id =
id << yaspgrid_dim_bits;
928 id =
id+PersistentIndexType(_it.coord(i)>>trailing);
935 int compressedIndex ()
const {
return _it.superindex();}
938 const I& transformingsubiterator()
const {
return _it; }
939 const YGLI& gridlevel()
const {
return _g; }
940 I& transformingsubiterator() {
return _it; }
941 YGLI& gridlevel() {
return _g; }
943 const GridImp * yaspgrid()
const {
return _g->mg; }
static constexpr int mydimension
geometry dimension
Definition: geometry.hh:94
persistent, globally unique Ids
Definition: yaspgrididset.hh:25
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
Implementation of Level- and LeafIndexSets for YaspGrid.
Definition: yaspgridindexsets.hh:25
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:587
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
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
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
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
A unique label for each type of element that can occur in a grid.