1 #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH
2 #define DUNE_ALBERTAGRIDINDEXSETS_HH
4 #include <dune/common/stdstreams.hh>
33 template<
int dim,
int dimworld >
35 :
public IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int >
61 struct CreateEntityNumbers;
64 struct RefineNumbering;
67 struct CoarsenNumbering;
75 template<
class Entity >
86 IndexType index (
const typename Traits::template Codim< cc >::Entity &entity )
const
90 return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
95 IndexType subIndex (
const typename Traits::template Codim< cc >::Entity &entity,
int i,
unsigned int codim )
const
103 const GenericReferenceElement< Alberta::Real, dimension > &refElement
105 k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
108 const int j = entityImp.grid().generic2alberta( codim, k );
109 return subIndex( entityImp.elementInfo(), j, codim );
115 return (type.isSimplex() ?
size(
dimension - type.dim() ) : 0);
121 assert( (codim >= 0) && (codim <=
dimension) );
122 return indexStack_[ codim ].
size();
126 const std::vector< GeometryType > &
geomTypes(
int codim )
const
128 assert( (codim >= 0) && (codim <=
dimension) );
129 return geomTypes_[ codim ];
134 assert( !elementInfo ==
false );
148 assert( (subIndex >= 0) && (subIndex <
size( codim )) );
170 void read (
const std::string &filename );
171 bool write (
const std::string &filename )
const;
180 template<
int codim >
181 static IndexStack &getIndexStack (
const IndexVectorPointer &dofVector )
185 indexStack = dofVector.template getAdaptationData< IndexStack >();
188 assert( indexStack != 0 );
199 IndexVectorPointer entityNumbers_[
dimension+1 ];
202 std::vector< GeometryType > geomTypes_[
dimension+1 ];
210 template<
int dim,
int dimworld >
217 : indexStack_( indexStack )
220 void operator() (
int &dof )
231 template<
int dim,
int dimworld >
232 template<
int codim >
240 static void apply (
const std::string &filename,
250 template<
int dim,
int dimworld >
251 template<
int codim >
252 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::RefineNumbering
255 static const int codimension = codim;
258 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
260 explicit RefineNumbering (
const IndexVectorPointer &dofVector )
261 : indexStack_( getIndexStack< codimension >( dofVector ) ),
262 dofVector_( dofVector ),
263 dofAccess_( dofVector.dofSpace() )
269 typedef Alberta::Patch< dimension > Patch;
270 static void interpolateVector (
const IndexVectorPointer &dofVector,
271 const Patch &patch );
275 IndexVectorPointer dofVector_;
276 DofAccess dofAccess_;
284 template<
int dim,
int dimworld >
285 template<
int codim >
286 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CoarsenNumbering
289 static const int codimension = codim;
292 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
294 explicit CoarsenNumbering (
const IndexVectorPointer &dofVector )
295 : indexStack_( getIndexStack< codimension >( dofVector ) ),
296 dofVector_( dofVector ),
297 dofAccess_( dofVector.dofSpace() )
303 typedef Alberta::Patch< dimension > Patch;
304 static void restrictVector (
const IndexVectorPointer &dofVector,
305 const Patch &patch );
308 IndexVectorPointer dofVector_;
309 DofAccess dofAccess_;
317 template<
int dim,
int dimworld >
318 class AlbertaGridIndexSet
319 :
public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int >
321 typedef AlbertaGridIndexSet< dim, dimworld > This;
322 typedef IndexSet< AlbertaGrid< dim, dimworld >, This,
int > Base;
337 template<
int codim >
342 : dofNumbering_( dofNumbering )
344 for(
int codim = 0; codim <=
dimension; ++codim )
346 indices_[ codim ] = 0;
349 geomTypes_[ codim ].push_back( type );
355 for(
int codim = 0; codim <=
dimension; ++codim )
356 delete[] indices_[ codim ];
359 template<
class Entity >
368 const IndexType *
const array = indices_[ codim ];
371 return (subIndex >= 0);
379 IndexType index (
const typename Traits::template Codim< cc >::Entity &entity )
const
383 return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
388 IndexType subIndex (
const typename Traits::template Codim< cc >::Entity &entity,
int i,
unsigned int codim )
const
396 const GenericReferenceElement< Alberta::Real, dimension > &refElement
398 k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
401 const int j = entityImp.grid().generic2alberta( codim, k );
402 return subIndex( entityImp.elementInfo(), j, codim );
407 return (type.isSimplex() ?
size(
dimension - type.dim() ) : 0);
412 assert( (codim >= 0) && (codim <=
dimension) );
413 return size_[ codim ];
416 const std::vector< GeometryType > &
geomTypes(
int codim )
const
418 assert( (codim >= 0) && (codim <=
dimension) );
419 return geomTypes_[ codim ];
422 template<
class Iterator >
423 void update (
const Iterator &begin,
const Iterator &end )
425 for(
int codim = 0; codim <=
dimension; ++codim )
427 delete[] indices_[ codim ];
429 const unsigned int dofSize = dofNumbering_.
size( codim );
430 indices_[ codim ] =
new IndexType[ dofSize ];
431 for(
unsigned int i = 0; i < dofSize; ++i )
432 indices_[ codim ][ i ] = -1;
437 for( Iterator it = begin; it != end; ++it )
442 ForLoop< Insert, 0, dimension >::apply( element, *
this );
449 assert( !elementInfo ==
false );
450 return subIndex( elementInfo.element(), i, codim );
461 const IndexType *
const array = indices_[ codim ];
463 assert( (subIndex >= 0) && (subIndex <
size( codim )) );
477 std::vector< GeometryType > geomTypes_[
dimension+1 ];
485 template<
int dim,
int dimworld >
486 template<
int codim >
487 struct AlbertaGridIndexSet< dim, dimworld >::Insert
490 AlbertaGridIndexSet< dim, dimworld > &indexSet )
492 int *
const array = indexSet.indices_[ codim ];
495 for(
int i = 0; i < Alberta::NumSubEntities< dim, codim >::value; ++i )
497 int &
index = array[ indexSet.dofNumbering_( element, codim, i ) ];
510 template<
int dim,
int dimworld >
511 class AlbertaGridIdSet
512 :
public IdSet< AlbertaGrid< dim, dimworld >, AlbertaGridIdSet< dim, dimworld >, unsigned int >
514 typedef AlbertaGridIdSet< dim, dimworld > This;
515 typedef IdSet< AlbertaGrid< dim, dimworld >, This,
unsigned int > Base;
532 : hIndexSet_( hIndexSet )
537 template<
class Entity >
541 return id< codim >( e );
545 template<
int codim >
546 IdType id (
const typename Grid::template Codim< codim >::Entity &e )
const
548 assert( (codim >= 0) && (codim <= dimension) );
550 return (index << 2) |
IdType( codim );
554 IdType subId (
const typename Grid::template Codim< 0 >::Entity &e,
int i,
unsigned int subcodim )
const
556 assert(
int( subcodim ) <= dimension );
558 return (index << 2) |
IdType( subcodim );
561 template<
int codim >
562 IdType subId (
const typename Grid::template Codim< codim >::Entity &e,
int i,
unsigned int subcodim )
const
564 assert( (codim >= 0) && (codim <= dimension) && (
int( codim + subcodim ) <= dimension) );
566 return (index << 2) |
IdType( codim + subcodim );
569 template<
class Entity >
572 return subId< Entity::codimension >( e, i, subcodim );
579 const HierarchicIndexSet &hIndexSet_;
584 #endif // #if HAVE_ALBERTA
586 #endif // #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH