4 #ifndef DUNE_ALBERTA_TREEITERATOR_HH
5 #define DUNE_ALBERTA_TREEITERATOR_HH
7 #include <dune/common/hybridutilities.hh>
8 #include <dune/common/std/utility.hh>
9 #include <dune/common/typetraits.hh>
30 template<
int dim,
int dimworld >
31 class AlbertaMarkerVector
33 typedef AlbertaMarkerVector< dim, dimworld > This;
35 typedef AlbertaGrid< dim, dimworld > Grid;
41 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
42 typedef Alberta::ElementInfo< dimension > ElementInfo;
45 struct NoMarkSubEntities;
47 struct MarkSubEntities;
51 explicit AlbertaMarkerVector (
const DofNumbering &dofNumbering )
52 : dofNumbering_( dofNumbering )
54 for(
int codim = 0; codim <= dimension; ++codim )
58 AlbertaMarkerVector (
const This &other )
59 : dofNumbering_( other.dofNumbering_ )
61 for(
int codim = 0; codim <= dimension; ++codim )
65 ~AlbertaMarkerVector ()
71 This &operator= (
const This & );
76 bool subEntityOnElement (
const ElementInfo &elementInfo,
int subEntity )
const;
78 template<
int firstCodim,
class Iterator >
79 void markSubEntities (
const Iterator &begin,
const Iterator &end );
83 for(
int codim = 0; codim <= dimension; ++codim )
85 if( marker_[ codim ] != 0 )
86 delete[] marker_[ codim ];
94 return (marker_[ dimension ] != 0);
98 void print ( std::ostream &out = std::cout )
const;
101 const DofNumbering &dofNumbering_;
102 int *marker_[ dimension+1 ];
110 template<
int dim,
int dimworld >
112 struct AlbertaMarkerVector< dim, dimworld >::NoMarkSubEntities
114 template<
int firstCodim,
class Iterator >
115 static void mark (
const DofNumbering &dofNumbering,
int *(&marker)[ dimension + 1 ],
116 const Iterator &begin,
const Iterator &end )
125 template<
int dim,
int dimworld >
127 struct AlbertaMarkerVector< dim, dimworld >::MarkSubEntities
129 template<
int codim >
132 static const int numSubEntities = Alberta::NumSubEntities< dimension, codim >::value;
134 typedef Alberta::ElementInfo< dimension > ElementInfo;
136 static void apply (
const DofNumbering &dofNumbering,
137 int *(&marker)[ dimension + 1 ],
138 const ElementInfo &elementInfo )
140 int *array = marker[ codim ];
142 const int index = dofNumbering( elementInfo, 0, 0 );
143 for(
int i = 0; i < numSubEntities; ++i )
145 int &mark = array[ dofNumbering( elementInfo, codim, i ) ];
146 mark = std::max( index, mark );
151 template<
int firstCodim,
class Iterator >
152 static void mark (
const DofNumbering &dofNumbering,
int *(&marker)[ dimension + 1 ],
153 const Iterator &begin,
const Iterator &end )
155 for(
int codim = firstCodim; codim <= dimension; ++codim )
157 const int size = dofNumbering.size( codim );
158 marker[ codim ] =
new int[ size ];
160 int *array = marker[ codim ];
161 for(
int i = 0; i < size; ++i )
165 for( Iterator it = begin; it != end; ++it )
167 const ElementInfo &elementInfo = Grid::getRealImplementation( *it ).elementInfo();
168 Hybrid::forEach( Std::make_index_sequence< dimension+1-firstCodim >{},
169 [ & ](
auto i ){ Codim< i+firstCodim >::apply( dofNumbering, marker, elementInfo ); } );
182 template<
int codim,
class Gr
idImp,
bool leafIterator >
183 class AlbertaGridTreeIterator
185 typedef AlbertaGridTreeIterator< codim, GridImp, leafIterator > This;
188 static const int dimension = GridImp::dimension;
189 static const int codimension = codim;
190 static const int dimensionworld = GridImp::dimensionworld;
193 friend class AlbertaGrid< dimension, dimensionworld >;
195 static const int numSubEntities
196 = Alberta::NumSubEntities< dimension, codimension >::value;
199 typedef Alberta::MeshPointer< dimension > MeshPointer;
200 typedef typename MeshPointer::MacroIterator MacroIterator;
202 typedef typename GridImp::template Codim< codim >::Entity Entity;
203 typedef MakeableInterfaceObject< Entity > EntityObject;
204 typedef typename EntityObject::ImplementationType EntityImp;
205 typedef typename EntityImp::ElementInfo ElementInfo;
207 typedef AlbertaMarkerVector< dimension, dimensionworld > MarkerVector;
209 AlbertaGridTreeIterator ();
212 AlbertaGridTreeIterator (
const This &other );
215 This &operator= (
const This &other );
218 AlbertaGridTreeIterator (
const GridImp &grid,
int travLevel );
221 AlbertaGridTreeIterator (
const GridImp &grid,
222 const MarkerVector *marker,
226 bool equals (
const This &other )
const
228 return entityImp().equals( other.entityImp() );
232 Entity &dereference ()
const
240 return entityImp().level();
248 EntityImp &entityImp ()
250 return GridImp::getRealImplementation( entity_ );
254 const EntityImp &entityImp ()
const
256 return GridImp::getRealImplementation( entity_ );
260 const GridImp &grid ()
const
262 return entityImp().grid();
266 void nextElement ( ElementInfo &elementInfo );
267 void nextElementStop (ElementInfo &elementInfo );
268 bool stopAtElement (
const ElementInfo &elementInfo )
const;
270 void goNext ( ElementInfo &elementInfo );
271 void goNext (
const std::integral_constant< int, 0 > cdVariable,
272 ElementInfo &elementInfo );
273 void goNext (
const std::integral_constant< int, 1 > cdVariable,
274 ElementInfo &elementInfo );
276 void goNext (
const std::integral_constant< int, cd > cdVariable,
277 ElementInfo &elementInfo );
279 mutable Entity entity_;
287 MacroIterator macroIterator_;
290 const MarkerVector *marker_;
298 template<
int dim,
int dimworld >
299 template<
int codim >
300 inline bool AlbertaMarkerVector< dim, dimworld >
301 ::subEntityOnElement (
const ElementInfo &elementInfo,
int subEntity )
const
303 assert( marker_[ codim ] != 0 );
305 const int subIndex = dofNumbering_( elementInfo, codim, subEntity );
306 const int markIndex = marker_[ codim ][ subIndex ];
307 assert( (markIndex >= 0) );
309 const int index = dofNumbering_( elementInfo, 0, 0 );
310 return (markIndex == index);
314 template<
int dim,
int dimworld >
315 template<
int firstCodim,
class Iterator >
316 inline void AlbertaMarkerVector< dim, dimworld >
317 ::markSubEntities (
const Iterator &begin,
const Iterator &end )
320 std::conditional< (firstCodim <= dimension), MarkSubEntities<true>, NoMarkSubEntities<false> >::type
321 ::template mark< firstCodim, Iterator >( dofNumbering_, marker_, begin, end );
325 template<
int dim,
int dimworld >
326 inline void AlbertaMarkerVector< dim, dimworld >::print ( std::ostream &out )
const
328 for(
int codim = 1; codim <= dimension; ++codim )
330 int *marker = marker_[ codim ];
333 const int size = dofNumbering_.size( codim );
335 out <<
"Codimension " << codim <<
" (" << size <<
" entries)" << std::endl;
336 for(
int i = 0; i < size; ++i )
337 out <<
"subentity " << i <<
" visited on Element " << marker[ i ] << std::endl;
347 template<
int codim,
class Gr
idImp,
bool leafIterator >
348 inline AlbertaGridTreeIterator< codim, GridImp, leafIterator >
349 ::AlbertaGridTreeIterator ()
357 template<
int codim,
class Gr
idImp,
bool leafIterator >
358 inline AlbertaGridTreeIterator< codim, GridImp, leafIterator >
359 ::AlbertaGridTreeIterator (
const GridImp &grid,
360 const MarkerVector *marker,
362 : entity_( EntityImp( grid ) ),
364 subEntity_( (codim == 0 ? 0 : -1) ),
365 macroIterator_( grid.meshPointer().begin() ),
368 ElementInfo elementInfo = *macroIterator_;
369 nextElementStop( elementInfo );
371 goNext( elementInfo );
373 entityImp().setElement( elementInfo, subEntity_ );
378 template<
int codim,
class Gr
idImp,
bool leafIterator >
379 inline AlbertaGridTreeIterator< codim, GridImp, leafIterator >
380 ::AlbertaGridTreeIterator (
const GridImp &grid,
382 : entity_( EntityImp( grid ) ),
385 macroIterator_( grid.meshPointer().end() ),
391 template<
int codim,
class Gr
idImp,
bool leafIterator >
392 inline AlbertaGridTreeIterator< codim, GridImp, leafIterator >
393 ::AlbertaGridTreeIterator(
const This &other )
394 : entity_( other.entity_ ),
395 level_( other.level_ ),
396 subEntity_( other.subEntity_ ),
397 macroIterator_( other.macroIterator_ ),
398 marker_( other.marker_ )
403 template<
int codim,
class Gr
idImp,
bool leafIterator >
404 inline typename AlbertaGridTreeIterator< codim, GridImp, leafIterator >::This &
405 AlbertaGridTreeIterator< codim, GridImp, leafIterator >::operator= (
const This &other )
407 entity_ = other.entity_;
408 level_ = other.level_;
409 subEntity_ = other.subEntity_;
410 macroIterator_ = other.macroIterator_;
411 marker_ = other.marker_;
417 template<
int codim,
class Gr
idImp,
bool leafIterator >
418 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >::increment ()
420 ElementInfo elementInfo = entityImp().elementInfo_;
421 goNext ( elementInfo );
423 entityImp().setElement( elementInfo, subEntity_ );
427 template<
int codim,
class Gr
idImp,
bool leafIterator >
428 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >
429 ::nextElement ( ElementInfo &elementInfo )
431 if( elementInfo.isLeaf() || (elementInfo.level() >= level_) )
433 while( (elementInfo.level() > 0) && (elementInfo.indexInFather() == 1) )
434 elementInfo = elementInfo.father();
435 if( elementInfo.level() == 0 )
438 elementInfo = *macroIterator_;
441 elementInfo = elementInfo.father().child( 1 );
444 elementInfo = elementInfo.child( 0 );
448 template<
int codim,
class Gr
idImp,
bool leafIterator >
449 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >
450 ::nextElementStop ( ElementInfo &elementInfo )
452 while( !(!elementInfo || stopAtElement( elementInfo )) )
453 nextElement( elementInfo );
457 template<
int codim,
class Gr
idImp,
bool leafIterator >
458 inline bool AlbertaGridTreeIterator< codim, GridImp, leafIterator >
459 ::stopAtElement (
const ElementInfo &elementInfo )
const
463 return (leafIterator ? elementInfo.isLeaf() : (level_ == elementInfo.level()));
467 template<
int codim,
class Gr
idImp,
bool leafIterator >
468 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >
469 ::goNext ( ElementInfo &elementInfo )
471 std::integral_constant< int, codim > codimVariable;
472 goNext( codimVariable, elementInfo );
475 template<
int codim,
class Gr
idImp,
bool leafIterator >
476 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >
477 ::goNext (
const std::integral_constant< int, 0 > cdVariable,
478 ElementInfo &elementInfo )
480 assert( stopAtElement( elementInfo ) );
482 nextElement( elementInfo );
483 nextElementStop( elementInfo );
486 template<
int codim,
class Gr
idImp,
bool leafIterator >
487 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >
488 ::goNext (
const std::integral_constant< int, 1 > cdVariable,
489 ElementInfo &elementInfo )
491 assert( stopAtElement( elementInfo ) );
494 if( subEntity_ >= numSubEntities )
497 nextElement( elementInfo );
498 nextElementStop( elementInfo );
505 const int face = (dimension == 1 ? (numSubEntities-1)-subEntity_ : subEntity_);
507 const ALBERTA EL *neighbor = elementInfo.elInfo().neigh[ face ];
508 if( (neighbor != NULL) && !elementInfo.isBoundary( face ) )
511 const int elIndex = grid().dofNumbering() ( elementInfo, 0, 0 );
512 const int nbIndex = grid().dofNumbering() ( neighbor, 0, 0 );
513 if( elIndex < nbIndex )
514 goNext( cdVariable, elementInfo );
521 assert( marker_ != 0 );
522 if( !marker_->template subEntityOnElement< 1 >( elementInfo, subEntity_ ) )
523 goNext( cdVariable, elementInfo );
527 template<
int codim,
class Gr
idImp,
bool leafIterator >
529 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >
530 ::goNext (
const std::integral_constant< int, cd > cdVariable,
531 ElementInfo &elementInfo )
533 assert( stopAtElement( elementInfo ) );
536 if( subEntity_ >= numSubEntities )
539 nextElement( elementInfo );
540 nextElementStop( elementInfo );
545 assert( marker_ != 0 );
546 if( !marker_->template subEntityOnElement< cd >( elementInfo, subEntity_ ) )
547 goNext( cdVariable, elementInfo );
552 #endif // #if HAVE_ALBERTA
554 #endif // #ifndef DUNE_ALBERTA_TREEITERATOR_HH