3 #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH
4 #define DUNE_ALBERTAGRIDINDEXSETS_HH
8 #include <dune/common/hybridutilities.hh>
9 #include <dune/common/stdstreams.hh>
10 #include <dune/common/std/utility.hh>
31 extern IndexStack *currentIndexStack;
39 template<
int dim,
int dimworld >
40 class AlbertaGridHierarchicIndexSet
41 :
public IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int, std::array< GeometryType, 1 > >
43 typedef AlbertaGridHierarchicIndexSet< dim, dimworld > This;
44 typedef IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int, std::array< GeometryType, 1 > > Base;
46 friend class AlbertaGrid< dim, dimworld >;
49 typedef AlbertaGrid< dim, dimworld > Grid;
50 typedef AlbertaGridFamily< dim, dimworld > GridFamily;
52 typedef typename Base::IndexType IndexType;
54 typedef typename Base::Types Types;
56 static const int dimension = GridFamily::dimension;
58 typedef Alberta::ElementInfo< dimension > ElementInfo;
59 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
62 typedef typename GridFamily::Traits Traits;
64 typedef Alberta::DofVectorPointer< IndexType > IndexVectorPointer;
66 class InitEntityNumber;
69 struct CreateEntityNumbers;
72 struct RefineNumbering;
75 struct CoarsenNumbering;
77 explicit AlbertaGridHierarchicIndexSet (
const DofNumbering &dofNumbering );
80 typedef Alberta::IndexStack IndexStack;
83 template<
class Entity >
84 bool contains (
const Entity & )
const
94 IndexType index (
const typename Traits::template Codim< cc >::Entity &entity )
const
96 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
97 const EntityImp &entityImp = entity.impl();
98 return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
103 IndexType subIndex (
const typename Traits::template Codim< cc >::Entity &entity,
int i,
unsigned int codim )
const
105 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
106 const EntityImp &entityImp = entity.impl();
111 auto refElement = ReferenceElements< Alberta::Real, dimension >::simplex();
112 k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
115 const int j = entityImp.grid().generic2alberta( codim, k );
116 return subIndex( entityImp.elementInfo(), j, codim );
122 return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
126 IndexType size (
int codim )
const
128 assert( (codim >= 0) && (codim <= dimension) );
129 return indexStack_[ codim ].size();
132 Types types (
int codim )
const
134 assert( (codim >= 0) && (codim <= dimension) );
135 return {{ GeometryTypes::simplex( dimension - codim ) }};
139 const std::vector< GeometryType > &geomTypes(
int codim )
const
141 assert( (codim >= 0) && (codim <= dimension) );
142 return geomTypes_[ codim ];
145 IndexType subIndex (
const ElementInfo &elementInfo,
int i,
unsigned int codim )
const
147 assert( !elementInfo ==
false );
148 return subIndex( elementInfo.element(), i, codim );
157 IndexType subIndex (
const Alberta::Element *element,
int i,
unsigned int codim )
const
159 IndexType *array = (IndexType *)entityNumbers_[ codim ];
160 const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
161 assert( (subIndex >= 0) && (subIndex < size( codim )) );
168 if( !IndexVectorPointer::supportsAdaptationData )
170 assert( Alberta::currentIndexStack == 0 );
171 Alberta::currentIndexStack = indexStack_;
178 if( !IndexVectorPointer::supportsAdaptationData )
179 Alberta::currentIndexStack = 0;
183 void read (
const std::string &filename );
184 bool write (
const std::string &filename )
const;
188 for(
int i = 0; i <= dimension; ++i )
189 entityNumbers_[ i ].release();
193 template<
int codim >
194 static IndexStack &getIndexStack (
const IndexVectorPointer &dofVector )
196 IndexStack *indexStack;
197 if( IndexVectorPointer::supportsAdaptationData )
198 indexStack = dofVector.template getAdaptationData< IndexStack >();
200 indexStack = &Alberta::currentIndexStack[ codim ];
201 assert( indexStack != 0 );
206 const DofNumbering &dofNumbering_;
209 IndexStack indexStack_[ dimension+1 ];
212 IndexVectorPointer entityNumbers_[ dimension+1 ];
215 std::vector< GeometryType > geomTypes_[ dimension+1 ];
223 template<
int dim,
int dimworld >
224 class AlbertaGridHierarchicIndexSet< dim, dimworld >::InitEntityNumber
226 IndexStack &indexStack_;
229 InitEntityNumber ( IndexStack &indexStack )
230 : indexStack_( indexStack )
233 void operator() (
int &dof )
235 dof = indexStack_.getIndex();
244 template<
int dim,
int dimworld >
245 template<
int codim >
246 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CreateEntityNumbers
248 static void setup ( AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
250 static void apply (
const Alberta::HierarchyDofNumbering< dimension > &dofNumbering,
251 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
253 static void apply (
const std::string &filename,
254 const Alberta::MeshPointer< dimension > &mesh,
255 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
263 template<
int dim,
int dimworld >
264 template<
int codim >
265 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::RefineNumbering
267 static const int dimension = dim;
268 static const int codimension = codim;
271 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
273 explicit RefineNumbering (
const IndexVectorPointer &dofVector )
274 : indexStack_( getIndexStack< codimension >( dofVector ) ),
275 dofVector_( dofVector ),
276 dofAccess_( dofVector.dofSpace() )
280 void operator() (
const Alberta::Element *child,
int subEntity );
282 typedef Alberta::Patch< dimension > Patch;
283 static void interpolateVector (
const IndexVectorPointer &dofVector,
284 const Patch &patch );
287 IndexStack &indexStack_;
288 IndexVectorPointer dofVector_;
289 DofAccess dofAccess_;
297 template<
int dim,
int dimworld >
298 template<
int codim >
299 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CoarsenNumbering
301 static const int dimension = dim;
302 static const int codimension = codim;
305 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
307 explicit CoarsenNumbering (
const IndexVectorPointer &dofVector )
308 : indexStack_( getIndexStack< codimension >( dofVector ) ),
309 dofVector_( dofVector ),
310 dofAccess_( dofVector.dofSpace() )
314 void operator() (
const Alberta::Element *child,
int subEntity );
316 typedef Alberta::Patch< dimension > Patch;
317 static void restrictVector (
const IndexVectorPointer &dofVector,
318 const Patch &patch );
320 IndexStack &indexStack_;
321 IndexVectorPointer dofVector_;
322 DofAccess dofAccess_;
330 template<
int dim,
int dimworld >
331 class AlbertaGridIndexSet
332 :
public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >
334 typedef AlbertaGridIndexSet< dim, dimworld > This;
335 typedef IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > > Base;
338 typedef AlbertaGrid< dim, dimworld > Grid;
340 typedef typename Base::IndexType IndexType;
342 typedef typename Base::Types Types;
346 typedef Alberta::ElementInfo< dimension > ElementInfo;
347 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
350 typedef typename Grid::Traits Traits;
352 template<
int codim >
356 explicit AlbertaGridIndexSet (
const DofNumbering &dofNumbering )
357 : dofNumbering_( dofNumbering )
359 for(
int codim = 0; codim <= dimension; ++codim )
361 indices_[ codim ] = 0;
362 geomTypes_[ codim ].push_back( GeometryTypes::simplex( dimension - codim ) );
366 ~AlbertaGridIndexSet ()
368 for(
int codim = 0; codim <= dimension; ++codim )
369 delete[] indices_[ codim ];
372 template<
class Entity >
373 bool contains (
const Entity &entity )
const
377 const AlbertaGridEntity< codim, dim, const Grid > &entityImp
379 const Alberta::Element *element = entityImp.elementInfo().el();
381 const IndexType *
const array = indices_[ codim ];
382 const IndexType subIndex = array[ dofNumbering_( element, codim, entityImp.subEntity() ) ];
384 return (subIndex >= 0);
388 using Base::subIndex;
392 IndexType index (
const typename Traits::template Codim< cc >::Entity &entity )
const
394 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
395 const EntityImp &entityImp = entity.impl();
396 return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
401 IndexType subIndex (
const typename Traits::template Codim< cc >::Entity &entity,
int i,
unsigned int codim )
const
403 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
404 const EntityImp &entityImp = entity.impl();
409 auto refElement = ReferenceElements< Alberta::Real, dimension >::simplex();
410 k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
413 const int j = entityImp.grid().generic2alberta( codim, k );
414 return subIndex( entityImp.elementInfo(), j, codim );
419 return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
422 IndexType size (
int codim )
const
424 assert( (codim >= 0) && (codim <= dimension) );
425 return size_[ codim ];
428 Types types (
int codim )
const
430 assert( (codim >= 0) && (codim <= dimension) );
431 return {{ GeometryTypes::simplex( dimension - codim ) }};
434 const std::vector< GeometryType > &geomTypes(
int codim )
const
436 assert( (codim >= 0) && (codim <= dimension) );
437 return geomTypes_[ codim ];
440 template<
class Iterator >
441 void update (
const Iterator &begin,
const Iterator &end )
443 for(
int codim = 0; codim <= dimension; ++codim )
445 delete[] indices_[ codim ];
447 const unsigned int dofSize = dofNumbering_.size( codim );
448 indices_[ codim ] =
new IndexType[ dofSize ];
449 for(
unsigned int i = 0; i < dofSize; ++i )
450 indices_[ codim ][ i ] = -1;
455 for( Iterator it = begin; it != end; ++it )
457 const AlbertaGridEntity< 0, dim, const Grid > &entityImp
459 const Alberta::Element *element = entityImp.elementInfo().el();
460 Hybrid::forEach( Std::make_index_sequence< dimension+1 >{},
461 [ & ](
auto i ){ Insert< i >::apply( element, *
this ); } );
466 IndexType subIndex (
const ElementInfo &elementInfo,
int i,
unsigned int codim )
const
468 assert( !elementInfo ==
false );
469 return subIndex( elementInfo.element(), i, codim );
478 IndexType subIndex (
const Alberta::Element *element,
int i,
unsigned int codim )
const
480 const IndexType *
const array = indices_[ codim ];
481 const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
482 assert( (subIndex >= 0) && (subIndex < size( codim )) );
487 const DofNumbering &dofNumbering_;
490 IndexType *indices_[ dimension+1 ];
493 IndexType size_[ dimension+1 ];
496 std::vector< GeometryType > geomTypes_[ dimension+1 ];
504 template<
int dim,
int dimworld >
505 template<
int codim >
506 struct AlbertaGridIndexSet< dim, dimworld >::Insert
508 static void apply (
const Alberta::Element *
const element,
509 AlbertaGridIndexSet< dim, dimworld > &indexSet )
511 int *
const array = indexSet.indices_[ codim ];
512 IndexType &size = indexSet.size_[ codim ];
514 for(
int i = 0; i < Alberta::NumSubEntities< dim, codim >::value; ++i )
516 int &index = array[ indexSet.dofNumbering_( element, codim, i ) ];
529 template<
int dim,
int dimworld >
530 class AlbertaGridIdSet
531 :
public IdSet< AlbertaGrid< dim, dimworld >, AlbertaGridIdSet< dim, dimworld >, unsigned int >
533 typedef AlbertaGridIdSet< dim, dimworld > This;
534 typedef IdSet< AlbertaGrid< dim, dimworld >, This,
unsigned int > Base;
536 friend class AlbertaGrid< dim, dimworld >;
540 typedef typename Base::IdType IdType;
543 typedef AlbertaGrid< dim, dimworld > Grid;
547 typedef typename Grid::HierarchicIndexSet HierarchicIndexSet;
550 AlbertaGridIdSet (
const HierarchicIndexSet &hIndexSet )
551 : hIndexSet_( hIndexSet )
556 template<
class Entity >
557 IdType id (
const Entity &e )
const
560 return id< codim >( e );
564 template<
int codim >
565 IdType id (
const typename Grid::template Codim< codim >::Entity &e )
const
567 assert( (codim >= 0) && (codim <= dimension) );
568 const IdType index = hIndexSet_.index( e );
569 return (index << 2) | IdType( codim );
573 IdType subId (
const typename Grid::template Codim< 0 >::Entity &e,
int i,
unsigned int subcodim )
const
575 assert(
int( subcodim ) <= dimension );
576 const IdType index = hIndexSet_.subIndex( e, i, subcodim );
577 return (index << 2) | IdType( subcodim );
580 template<
int codim >
581 IdType subId (
const typename Grid::template Codim< codim >::Entity &e,
int i,
unsigned int subcodim )
const
583 assert( (codim >= 0) && (codim <= dimension) && (
int( codim + subcodim ) <= dimension) );
584 const IdType index = hIndexSet_.subIndex( e, i, subcodim );
585 return (index << 2) | IdType( codim + subcodim );
588 template<
class Entity >
589 IdType subId (
const Entity &e,
int i,
unsigned int subcodim )
const
591 return subId< Entity::codimension >( e, i, subcodim );
596 AlbertaGridIdSet (
const This & );
598 const HierarchicIndexSet &hIndexSet_;
603 #endif // #if HAVE_ALBERTA
605 #endif // #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH