dune-grid
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=8 sw=2 sts=2: 00003 00004 #ifndef DUNE_ALBERTA_TREEITERATOR_HH 00005 #define DUNE_ALBERTA_TREEITERATOR_HH 00006 00007 #include <dune/common/typetraits.hh> 00008 00009 #include <dune/grid/albertagrid/meshpointer.hh> 00010 #include <dune/grid/albertagrid/entitypointer.hh> 00011 00012 #if HAVE_ALBERTA 00013 00014 namespace Dune 00015 { 00016 00017 // AlbertaMarkerVector 00018 // ------------------- 00019 00028 template< int dim, int dimworld > 00029 class AlbertaMarkerVector 00030 { 00031 typedef AlbertaMarkerVector< dim, dimworld > This; 00032 00033 typedef AlbertaGrid< dim, dimworld > Grid; 00034 00035 //friend class AlbertaGrid< dim, dimworld >; 00036 00037 static const int dimension = Grid::dimension; 00038 00039 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering; 00040 typedef Alberta::ElementInfo< dimension > ElementInfo; 00041 00042 template< bool > 00043 struct NoMarkSubEntities; 00044 template< bool > 00045 struct MarkSubEntities; 00046 00047 public: 00049 explicit AlbertaMarkerVector ( const DofNumbering &dofNumbering ) 00050 : dofNumbering_( dofNumbering ) 00051 { 00052 for( int codim = 0; codim <= dimension; ++codim ) 00053 marker_[ codim ] = 0; 00054 } 00055 00056 AlbertaMarkerVector ( const This &other ) 00057 : dofNumbering_( other.dofNumbering_ ) 00058 { 00059 for( int codim = 0; codim <= dimension; ++codim ) 00060 marker_[ codim ] = 0; 00061 } 00062 00063 ~AlbertaMarkerVector () 00064 { 00065 clear(); 00066 } 00067 00068 private: 00069 This &operator= ( const This & ); 00070 00071 public: 00073 template< int codim > 00074 bool subEntityOnElement ( const ElementInfo &elementInfo, int subEntity ) const; 00075 00076 template< int firstCodim, class Iterator > 00077 void markSubEntities ( const Iterator &begin, const Iterator &end ); 00078 00079 void clear () 00080 { 00081 for( int codim = 0; codim <= dimension; ++codim ) 00082 { 00083 if( marker_[ codim ] != 0 ) 00084 delete[] marker_[ codim ]; 00085 marker_[ codim ] = 0; 00086 } 00087 } 00088 00090 bool up2Date () const 00091 { 00092 return (marker_[ dimension ] != 0); 00093 } 00094 00096 void print ( std::ostream &out = std::cout ) const; 00097 00098 private: 00099 const DofNumbering &dofNumbering_; 00100 int *marker_[ dimension+1 ]; 00101 }; 00102 00103 00104 00105 // AlbertaMarkerVector::NoMarkSubEntities 00106 // -------------------------------------- 00107 00108 template< int dim, int dimworld > 00109 template< bool > 00110 struct AlbertaMarkerVector< dim, dimworld >::NoMarkSubEntities 00111 { 00112 template< int firstCodim, class Iterator > 00113 static void mark ( const DofNumbering &dofNumbering, int *(&marker)[ dimension + 1 ], 00114 const Iterator &begin, const Iterator &end ) 00115 {} 00116 }; 00117 00118 00119 00120 // AlbertaMarkerVector::MarkSubEntities 00121 // ------------------------------------ 00122 00123 template< int dim, int dimworld > 00124 template< bool > 00125 struct AlbertaMarkerVector< dim, dimworld >::MarkSubEntities 00126 { 00127 template< int codim > 00128 struct Codim 00129 { 00130 static const int numSubEntities = Alberta::NumSubEntities< dimension, codim >::value; 00131 00132 typedef Alberta::ElementInfo< dimension > ElementInfo; 00133 00134 static void apply ( const DofNumbering &dofNumbering, 00135 int *(&marker)[ dimension + 1 ], 00136 const ElementInfo &elementInfo ) 00137 { 00138 int *array = marker[ codim ]; 00139 00140 const int index = dofNumbering( elementInfo, 0, 0 ); 00141 for( int i = 0; i < numSubEntities; ++i ) 00142 { 00143 int &mark = array[ dofNumbering( elementInfo, codim, i ) ]; 00144 mark = std::max( index, mark ); 00145 } 00146 } 00147 }; 00148 00149 template< int firstCodim, class Iterator > 00150 static void mark ( const DofNumbering &dofNumbering, int *(&marker)[ dimension + 1 ], 00151 const Iterator &begin, const Iterator &end ) 00152 { 00153 for( int codim = firstCodim; codim <= dimension; ++codim ) 00154 { 00155 const int size = dofNumbering.size( codim ); 00156 marker[ codim ] = new int[ size ]; 00157 00158 int *array = marker[ codim ]; 00159 for( int i = 0; i < size; ++i ) 00160 array[ i ] = -1; 00161 } 00162 00163 for( Iterator it = begin; it != end; ++it ) 00164 { 00165 const ElementInfo &elementInfo = Grid::getRealImplementation( *it ).elementInfo(); 00166 ForLoop< Codim, firstCodim, dimension >::apply( dofNumbering, marker, elementInfo ); 00167 } 00168 } 00169 }; 00170 00171 00172 00173 // AlbertaGridTreeIterator 00174 // ----------------------- 00175 00179 template< int codim, class GridImp, bool leafIterator > 00180 class AlbertaGridTreeIterator 00181 : public AlbertaGridEntityPointer< codim, GridImp > 00182 { 00183 typedef AlbertaGridTreeIterator< codim, GridImp, leafIterator > This; 00184 typedef AlbertaGridEntityPointer< codim, GridImp > Base; 00185 00186 public: 00187 static const int dimension = GridImp::dimension; 00188 static const int codimension = codim; 00189 static const int dimensionworld = GridImp::dimensionworld; 00190 00191 private: 00192 friend class AlbertaGrid< dimension, dimensionworld >; 00193 00194 static const int numSubEntities 00195 = Alberta::NumSubEntities< dimension, codimension >::value; 00196 00197 public: 00198 typedef typename Base::ElementInfo ElementInfo; 00199 typedef Alberta::MeshPointer< dimension > MeshPointer; 00200 typedef typename MeshPointer::MacroIterator MacroIterator; 00201 00202 typedef typename GridImp::template Codim< codim >::Entity Entity; 00203 typedef MakeableInterfaceObject< Entity > EntityObject; 00204 typedef typename EntityObject::ImplementationType EntityImp; 00205 00206 typedef AlbertaMarkerVector< dimension, dimensionworld > MarkerVector; 00207 00209 AlbertaGridTreeIterator ( const This &other ); 00210 00212 This &operator= ( const This &other ); 00213 00215 AlbertaGridTreeIterator ( const GridImp &grid, int travLevel ); 00216 00218 AlbertaGridTreeIterator ( const GridImp &grid, 00219 const MarkerVector *marker, 00220 int travLevel ); 00221 00223 void increment(); 00224 00225 protected: 00226 using Base::entityImp; 00227 using Base::grid; 00228 00229 private: 00230 void nextElement ( ElementInfo &elementInfo ); 00231 void nextElementStop (ElementInfo &elementInfo ); 00232 bool stopAtElement ( const ElementInfo &elementInfo ) const; 00233 00234 void goNext ( ElementInfo &elementInfo ); 00235 void goNext ( const integral_constant< int, 0 > cdVariable, 00236 ElementInfo &elementInfo ); 00237 void goNext ( const integral_constant< int, 1 > cdVariable, 00238 ElementInfo &elementInfo ); 00239 template< int cd > 00240 void goNext ( const integral_constant< int, cd > cdVariable, 00241 ElementInfo &elementInfo ); 00242 00244 int level_; 00245 00247 int subEntity_; 00248 00249 MacroIterator macroIterator_; 00250 00251 // knows on which element a point,edge,face is viewed 00252 const MarkerVector *marker_; 00253 }; 00254 00255 00256 00257 // Implementation of AlbertaMarkerVector 00258 // ------------------------------------- 00259 00260 template< int dim, int dimworld > 00261 template< int codim > 00262 inline bool AlbertaMarkerVector< dim, dimworld > 00263 ::subEntityOnElement ( const ElementInfo &elementInfo, int subEntity ) const 00264 { 00265 assert( marker_[ codim ] != 0 ); 00266 00267 const int subIndex = dofNumbering_( elementInfo, codim, subEntity ); 00268 const int markIndex = marker_[ codim ][ subIndex ]; 00269 assert( (markIndex >= 0) ); 00270 00271 const int index = dofNumbering_( elementInfo, 0, 0 ); 00272 return (markIndex == index); 00273 } 00274 00275 00276 template< int dim, int dimworld > 00277 template< int firstCodim, class Iterator > 00278 inline void AlbertaMarkerVector< dim, dimworld > 00279 ::markSubEntities ( const Iterator &begin, const Iterator &end ) 00280 { 00281 clear(); 00282 SelectType< (firstCodim <= dimension), MarkSubEntities<true>, NoMarkSubEntities<false> >::Type 00283 ::template mark< firstCodim, Iterator >( dofNumbering_, marker_, begin, end ); 00284 } 00285 00286 00287 template< int dim, int dimworld > 00288 inline void AlbertaMarkerVector< dim, dimworld >::print ( std::ostream &out ) const 00289 { 00290 for( int codim = 1; codim <= dimension; ++codim ) 00291 { 00292 int *marker = marker_[ codim ]; 00293 if( marker != 0 ) 00294 { 00295 const int size = dofNumbering_.size( codim ); 00296 out << std::endl; 00297 out << "Codimension " << codim << " (" << size << " entries)" << std::endl; 00298 for( int i = 0; i < size; ++i ) 00299 out << "subentity " << i << " visited on Element " << marker[ i ] << std::endl; 00300 } 00301 } 00302 } 00303 00304 00305 00306 // Implementation of AlbertaGridTreeIterator 00307 // ----------------------------------------- 00308 00309 template< int codim, class GridImp, bool leafIterator > 00310 inline AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00311 ::AlbertaGridTreeIterator ( const GridImp &grid, 00312 const MarkerVector *marker, 00313 int travLevel ) 00314 : Base( grid ), 00315 level_( travLevel ), 00316 subEntity_( (codim == 0 ? 0 : -1) ), 00317 macroIterator_( grid.meshPointer().begin() ), 00318 marker_( marker ) 00319 { 00320 ElementInfo elementInfo = *macroIterator_; 00321 nextElementStop( elementInfo ); 00322 if( codim > 0 ) 00323 goNext( elementInfo ); 00324 // it is ok to set the invalid ElementInfo 00325 entityImp().setElement( elementInfo, subEntity_ ); 00326 } 00327 00328 00329 // Make LevelIterator with point to element from previous iterations 00330 template< int codim, class GridImp, bool leafIterator > 00331 inline AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00332 ::AlbertaGridTreeIterator ( const GridImp &grid, 00333 int travLevel ) 00334 : Base( grid ), 00335 level_( travLevel ), 00336 subEntity_( -1 ), 00337 macroIterator_( grid.meshPointer().end() ), 00338 marker_( 0 ) 00339 {} 00340 00341 00342 // Make LevelIterator with point to element from previous iterations 00343 template< int codim, class GridImp, bool leafIterator > 00344 inline AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00345 ::AlbertaGridTreeIterator( const This &other ) 00346 : Base( other ), 00347 level_( other.level_ ), 00348 subEntity_( other.subEntity_ ), 00349 macroIterator_( other.macroIterator_ ), 00350 marker_( other.marker_ ) 00351 {} 00352 00353 00354 // Make LevelIterator with point to element from previous iterations 00355 template< int codim, class GridImp, bool leafIterator > 00356 inline typename AlbertaGridTreeIterator< codim, GridImp, leafIterator >::This & 00357 AlbertaGridTreeIterator< codim, GridImp, leafIterator >::operator= ( const This &other ) 00358 { 00359 Base::operator=( other ); 00360 00361 level_ = other.level_; 00362 subEntity_ = other.subEntity_; 00363 macroIterator_ = other.macroIterator_; 00364 marker_ = other.marker_; 00365 00366 return *this; 00367 } 00368 00369 00370 template< int codim, class GridImp, bool leafIterator > 00371 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator >::increment () 00372 { 00373 ElementInfo elementInfo = entityImp().elementInfo_; 00374 goNext ( elementInfo ); 00375 // it is ok to set the invalid ElementInfo 00376 entityImp().setElement( elementInfo, subEntity_ ); 00377 } 00378 00379 00380 template< int codim, class GridImp, bool leafIterator > 00381 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00382 ::nextElement ( ElementInfo &elementInfo ) 00383 { 00384 if( elementInfo.isLeaf() || (elementInfo.level() >= level_) ) 00385 { 00386 while( (elementInfo.level() > 0) && (elementInfo.indexInFather() == 1) ) 00387 elementInfo = elementInfo.father(); 00388 if( elementInfo.level() == 0 ) 00389 { 00390 ++macroIterator_; 00391 elementInfo = *macroIterator_; 00392 } 00393 else 00394 elementInfo = elementInfo.father().child( 1 ); 00395 } 00396 else 00397 elementInfo = elementInfo.child( 0 ); 00398 } 00399 00400 00401 template< int codim, class GridImp, bool leafIterator > 00402 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00403 ::nextElementStop ( ElementInfo &elementInfo ) 00404 { 00405 while( !(!elementInfo || stopAtElement( elementInfo )) ) 00406 nextElement( elementInfo ); 00407 } 00408 00409 00410 template< int codim, class GridImp, bool leafIterator > 00411 inline bool AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00412 ::stopAtElement ( const ElementInfo &elementInfo ) const 00413 { 00414 if( !elementInfo ) 00415 return true; 00416 return (leafIterator ? elementInfo.isLeaf() : (level_ == elementInfo.level())); 00417 } 00418 00419 00420 template< int codim, class GridImp, bool leafIterator > 00421 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00422 ::goNext ( ElementInfo &elementInfo ) 00423 { 00424 integral_constant< int, codim > codimVariable; 00425 goNext( codimVariable, elementInfo ); 00426 } 00427 00428 template< int codim, class GridImp, bool leafIterator > 00429 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00430 ::goNext ( const integral_constant< int, 0 > cdVariable, 00431 ElementInfo &elementInfo ) 00432 { 00433 assert( stopAtElement( elementInfo ) ); 00434 00435 nextElement( elementInfo ); 00436 nextElementStop( elementInfo ); 00437 } 00438 00439 template< int codim, class GridImp, bool leafIterator > 00440 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00441 ::goNext ( const integral_constant< int, 1 > cdVariable, 00442 ElementInfo &elementInfo ) 00443 { 00444 assert( stopAtElement( elementInfo ) ); 00445 00446 ++subEntity_; 00447 if( subEntity_ >= numSubEntities ) 00448 { 00449 subEntity_ = 0; 00450 nextElement( elementInfo ); 00451 nextElementStop( elementInfo ); 00452 if( !elementInfo ) 00453 return; 00454 } 00455 00456 if( leafIterator ) 00457 { 00458 const int face = (dimension == 1 ? (numSubEntities-1)-subEntity_ : subEntity_); 00459 00460 const ALBERTA EL *neighbor = elementInfo.elInfo().neigh[ face ]; 00461 if( (neighbor != NULL) && !elementInfo.isBoundary( face ) ) 00462 { 00463 // face is reached from element with largest number 00464 const int elIndex = grid().dofNumbering()( elementInfo, 0, 0 ); 00465 const int nbIndex = grid().dofNumbering()( neighbor, 0, 0 ); 00466 if( elIndex < nbIndex ) 00467 goNext( cdVariable, elementInfo ); 00468 } 00469 // uncomment this assertion only if codimension 1 entities are marked 00470 // assert( marker_->template subEntityOnElement< 1 >( elementInfo, subEntity_ ) ); 00471 } 00472 else 00473 { 00474 assert( marker_ != 0 ); 00475 if( !marker_->template subEntityOnElement< 1 >( elementInfo, subEntity_ ) ) 00476 goNext( cdVariable, elementInfo ); 00477 } 00478 } 00479 00480 template< int codim, class GridImp, bool leafIterator > 00481 template< int cd > 00482 inline void AlbertaGridTreeIterator< codim, GridImp, leafIterator > 00483 ::goNext ( const integral_constant< int, cd > cdVariable, 00484 ElementInfo &elementInfo ) 00485 { 00486 assert( stopAtElement( elementInfo ) ); 00487 00488 ++subEntity_; 00489 if( subEntity_ >= numSubEntities ) 00490 { 00491 subEntity_ = 0; 00492 nextElement( elementInfo ); 00493 nextElementStop( elementInfo ); 00494 if( !elementInfo ) 00495 return; 00496 } 00497 00498 assert( marker_ != 0 ); 00499 if( !marker_->template subEntityOnElement< cd >( elementInfo, subEntity_ ) ) 00500 goNext( cdVariable, elementInfo ); 00501 } 00502 00503 } 00504 00505 #endif // #if HAVE_ALBERTA 00506 00507 #endif // #ifndef DUNE_ALBERTA_TREEITERATOR_HH