1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=4 sw=2 sts=2: 3 #ifndef DUNE_ALBERTAGRID_CC 4 #define DUNE_ALBERTAGRID_CC 5 6 //************************************************************************ 7 // 8 // implementation of AlbertaGrid 9 // 10 // namespace Dune 11 // 12 //************************************************************************ 13 #include "geometry.cc" 14 #include "entity.cc" 15 #include "intersection.cc" 16 17 namespace Dune 18 { 19 20 namespace Alberta 21 { 22 static void *adaptationDataHandler_; 23 } 24 25 26 template< int dim, int dimworld > checkAlbertaDimensions()27 static void checkAlbertaDimensions () 28 { 29 // If this check fails, define ALBERTA_DIM accordingly 30 static_assert((dimworld == Alberta::dimWorld), 31 "Template Parameter dimworld does not match " 32 "ALBERTA's DIM_OF_WORLD setting."); 33 } 34 35 36 // AlbertaGrid 37 // ----------- 38 39 template< int dim, int dimworld > AlbertaGrid()40 inline AlbertaGrid < dim, dimworld >::AlbertaGrid () 41 : mesh_(), 42 maxlevel_( 0 ), 43 numBoundarySegments_( 0 ), 44 hIndexSet_( dofNumbering_ ), 45 idSet_( hIndexSet_ ), 46 levelIndexVec_( (size_t)MAXL, 0 ), 47 leafIndexSet_( 0 ), 48 sizeCache_( *this ), 49 leafMarkerVector_( dofNumbering_ ), 50 levelMarkerVector_( (size_t)MAXL, MarkerVector( dofNumbering_ ) ) 51 { 52 checkAlbertaDimensions< dim, dimworld>(); 53 } 54 55 56 template< int dim, int dimworld > 57 template< class Proj, class Impl > 58 inline AlbertaGrid< dim, dimworld > AlbertaGrid(const Alberta::MacroData<dimension> & macroData,const Alberta::ProjectionFactoryInterface<Proj,Impl> & projectionFactory)59 ::AlbertaGrid ( const Alberta::MacroData< dimension> ¯oData, 60 const Alberta::ProjectionFactoryInterface< Proj, Impl > &projectionFactory ) 61 : mesh_(), 62 maxlevel_( 0 ), 63 numBoundarySegments_( 0 ), 64 hIndexSet_( dofNumbering_ ), 65 idSet_( hIndexSet_ ), 66 levelIndexVec_( (size_t)MAXL, 0 ), 67 leafIndexSet_ ( 0 ), 68 sizeCache_( *this ), 69 leafMarkerVector_( dofNumbering_ ), 70 levelMarkerVector_( (size_t)MAXL, MarkerVector( dofNumbering_ ) ) 71 { 72 checkAlbertaDimensions< dim, dimworld >(); 73 74 numBoundarySegments_ = mesh_.create( macroData, projectionFactory ); 75 if( !mesh_ ) 76 DUNE_THROW( AlbertaError, "Invalid macro data structure." ); 77 78 setup(); 79 hIndexSet_.create(); 80 81 calcExtras(); 82 } 83 84 85 template< int dim, int dimworld > 86 inline AlbertaGrid< dim, dimworld > AlbertaGrid(const Alberta::MacroData<dimension> & macroData,const std::shared_ptr<DuneBoundaryProjection<dimensionworld>> & projection)87 ::AlbertaGrid ( const Alberta::MacroData< dimension> ¯oData, 88 const std::shared_ptr< DuneBoundaryProjection< dimensionworld > > &projection ) 89 : mesh_(), 90 maxlevel_( 0 ), 91 numBoundarySegments_( 0 ), 92 hIndexSet_( dofNumbering_ ), 93 idSet_( hIndexSet_ ), 94 levelIndexVec_( (size_t)MAXL, 0 ), 95 leafIndexSet_ ( 0 ), 96 sizeCache_( *this ), 97 leafMarkerVector_( dofNumbering_ ), 98 levelMarkerVector_( (size_t)MAXL, MarkerVector( dofNumbering_ ) ) 99 { 100 checkAlbertaDimensions< dim, dimworld >(); 101 102 if( projection != 0 ) 103 { 104 Alberta::DuneGlobalBoundaryProjectionFactory< dimension > projectionFactory( projection ); 105 numBoundarySegments_ = mesh_.create( macroData, projectionFactory ); 106 } 107 else 108 numBoundarySegments_ = mesh_.create( macroData ); 109 if( !mesh_ ) 110 DUNE_THROW( AlbertaError, "Invalid macro data structure." ); 111 112 setup(); 113 hIndexSet_.create(); 114 115 calcExtras(); 116 } 117 118 119 template < int dim, int dimworld > 120 inline AlbertaGrid< dim, dimworld > AlbertaGrid(const std::string & macroGridFileName)121 ::AlbertaGrid ( const std::string ¯oGridFileName ) 122 : mesh_(), 123 maxlevel_( 0 ), 124 hIndexSet_( dofNumbering_ ), 125 idSet_( hIndexSet_ ), 126 levelIndexVec_( (size_t)MAXL, 0 ), 127 leafIndexSet_ ( 0 ), 128 sizeCache_( *this ), 129 leafMarkerVector_( dofNumbering_ ), 130 levelMarkerVector_( (size_t)MAXL, MarkerVector( dofNumbering_ ) ) 131 { 132 checkAlbertaDimensions< dim, dimworld >(); 133 134 numBoundarySegments_ = mesh_.create( macroGridFileName ); 135 if( !mesh_ ) 136 { 137 DUNE_THROW( AlbertaIOError, 138 "Grid file '" << macroGridFileName 139 << "' is not in ALBERTA macro triangulation format." ); 140 } 141 142 setup(); 143 hIndexSet_.create(); 144 145 calcExtras(); 146 147 std::cout << typeName() << " created from macro grid file '" 148 << macroGridFileName << "'." << std::endl; 149 } 150 151 152 template< int dim, int dimworld > setup()153 inline void AlbertaGrid< dim, dimworld >::setup () 154 { 155 dofNumbering_.create( mesh_ ); 156 157 levelProvider_.create( dofNumbering_ ); 158 159 #if DUNE_ALBERTA_CACHE_COORDINATES 160 coordCache_.create( dofNumbering_ ); 161 #endif 162 } 163 164 165 template< int dim, int dimworld > removeMesh()166 inline void AlbertaGrid< dim, dimworld >::removeMesh () 167 { 168 for( size_t i = 0; i < levelIndexVec_.size(); ++i ) 169 { 170 if( levelIndexVec_[ i ] != 0 ) 171 delete levelIndexVec_[ i ]; 172 levelIndexVec_[ i ] = 0; 173 } 174 175 if( leafIndexSet_ != 0 ) 176 delete leafIndexSet_; 177 leafIndexSet_ = 0; 178 179 // release dof vectors 180 hIndexSet_.release(); 181 levelProvider_.release(); 182 #if DUNE_ALBERTA_CACHE_COORDINATES 183 coordCache_.release(); 184 #endif 185 dofNumbering_.release(); 186 187 sizeCache_.reset(); 188 189 mesh_.release(); 190 } 191 192 193 template< int dim, int dimworld > ~AlbertaGrid()194 inline AlbertaGrid< dim, dimworld >::~AlbertaGrid () 195 { 196 removeMesh(); 197 } 198 199 200 template< int dim, int dimworld > 201 template< int codim, PartitionIteratorType pitype > 202 inline typename AlbertaGrid< dim, dimworld >::Traits 203 ::template Codim< codim >::template Partition<pitype>::LevelIterator lbegin(int level) const204 AlbertaGrid< dim, dimworld >::lbegin ( int level ) const 205 { 206 typedef AlbertaGridLevelIterator< codim, pitype, const This > LevelIteratorImp; 207 assert( level >= 0 ); 208 209 if( level > maxlevel_ ) 210 return lend< codim, pitype >( level ); 211 MarkerVector &markerVector = levelMarkerVector_[ level ]; 212 213 if( (codim > 0) && !markerVector.up2Date() ) 214 markerVector.template markSubEntities< 1 >( lbegin< 0 >( level ), lend< 0 >( level ) ); 215 216 return LevelIteratorImp( *this, &markerVector, level ); 217 } 218 219 220 template< int dim, int dimworld > 221 template< int codim, PartitionIteratorType pitype > 222 inline typename AlbertaGrid< dim, dimworld >::Traits 223 ::template Codim< codim >::template Partition< pitype >::LevelIterator lend(int level) const224 AlbertaGrid < dim, dimworld >::lend ( int level ) const 225 { 226 typedef AlbertaGridLevelIterator< codim, pitype, const This > LevelIteratorImp; 227 assert( level >= 0 ); 228 229 return LevelIteratorImp( *this, level ); 230 } 231 232 233 template< int dim, int dimworld > 234 template< int codim > 235 inline typename AlbertaGrid< dim, dimworld >::Traits 236 ::template Codim< codim >::LevelIterator lbegin(int level) const237 AlbertaGrid < dim, dimworld >::lbegin ( int level ) const 238 { 239 return lbegin< codim, All_Partition >( level ); 240 } 241 242 243 template< int dim, int dimworld > 244 template< int codim > 245 inline typename AlbertaGrid< dim, dimworld >::Traits 246 ::template Codim< codim >::LevelIterator lend(int level) const247 AlbertaGrid< dim, dimworld >::lend ( int level ) const 248 { 249 return lend< codim, All_Partition >( level ); 250 } 251 252 253 template< int dim, int dimworld > 254 template< int codim, PartitionIteratorType pitype > 255 inline typename AlbertaGrid< dim, dimworld >::Traits 256 ::template Codim< codim >::template Partition< pitype >::LeafIterator leafbegin() const257 AlbertaGrid< dim, dimworld >::leafbegin () const 258 { 259 typedef AlbertaGridLeafIterator< codim, pitype, const This > LeafIteratorImp; 260 261 MarkerVector &markerVector = leafMarkerVector_; 262 const int firstMarkedCodim = 2; 263 if( (codim >= firstMarkedCodim) && !markerVector.up2Date() ) 264 markerVector.template markSubEntities< firstMarkedCodim >( leafbegin< 0 >(), leafend< 0 >() ); 265 266 return LeafIteratorImp( *this, &markerVector, maxlevel_ ); 267 } 268 269 270 template< int dim, int dimworld > 271 template< int codim, PartitionIteratorType pitype > 272 inline typename AlbertaGrid< dim, dimworld >::Traits 273 ::template Codim< codim >::template Partition< pitype >::LeafIterator leafend() const274 AlbertaGrid< dim, dimworld >::leafend () const 275 { 276 typedef AlbertaGridLeafIterator< codim, pitype, const This > LeafIteratorImp; 277 return LeafIteratorImp( *this, maxlevel_ ); 278 } 279 280 281 template< int dim, int dimworld > 282 template< int codim > 283 inline typename AlbertaGrid<dim,dimworld>::Traits 284 ::template Codim< codim >::LeafIterator leafbegin() const285 AlbertaGrid< dim, dimworld >::leafbegin () const 286 { 287 return leafbegin< codim, All_Partition >(); 288 } 289 290 291 template< int dim, int dimworld > 292 template< int codim > 293 inline typename AlbertaGrid< dim, dimworld >::Traits 294 ::template Codim< codim >::LeafIterator leafend() const295 AlbertaGrid < dim, dimworld >::leafend () const 296 { 297 return leafend< codim, All_Partition >(); 298 } 299 300 301 template< int dim, int dimworld > globalRefine(int refCount)302 inline void AlbertaGrid< dim, dimworld >::globalRefine ( int refCount ) 303 { 304 typedef typename Traits::template Codim< 0 >::LeafIterator LeafIterator; 305 306 // only MAXL levels allowed 307 assert( (refCount >= 0) && (refCount + maxlevel_ < MAXL) ); 308 309 for( int i = 0; i < refCount; ++i ) 310 { 311 // mark all interior elements 312 const LeafIterator endit = leafend< 0 >(); 313 for( LeafIterator it = leafbegin< 0 >(); it != endit; ++it ) 314 mark( 1, *it ); 315 316 preAdapt(); 317 adapt(); 318 postAdapt(); 319 } 320 } 321 322 323 template< int dim, int dimworld > 324 template< class DataHandle > 325 inline void AlbertaGrid< dim, dimworld > globalRefine(int refCount,AdaptDataHandleInterface<This,DataHandle> & handle)326 ::globalRefine ( int refCount, AdaptDataHandleInterface< This, DataHandle > &handle ) 327 { 328 typedef typename Traits::template Codim< 0 >::LeafIterator LeafIterator; 329 330 // only MAXL levels allowed 331 assert( (refCount >= 0) && (refCount + maxlevel_ < MAXL) ); 332 333 for( int i = 0; i < refCount; ++i ) 334 { 335 // mark all interior elements 336 const LeafIterator endit = leafend< 0 >(); 337 for( LeafIterator it = leafbegin< 0 >(); it != endit; ++it ) 338 mark( 1, *it ); 339 340 adapt( handle ); 341 } 342 } 343 344 345 template< int dim, int dimworld > preAdapt()346 inline bool AlbertaGrid< dim, dimworld >::preAdapt () 347 { 348 adaptationState_.preAdapt(); 349 return adaptationState_.coarsen(); 350 } 351 352 353 template < int dim, int dimworld > postAdapt()354 inline void AlbertaGrid < dim, dimworld >::postAdapt () 355 { 356 #ifndef NDEBUG 357 if( leafIndexSet_ != 0 ) 358 { 359 bool consistent = true; 360 for( int codim = 0; codim <= dimension; ++codim ) 361 { 362 if( int(leafIndexSet_->size( codim )) == mesh_.size( codim ) ) 363 continue; 364 std::cerr << "Incorrect size of leaf index set for codimension " 365 << codim << "!" << std::endl; 366 std::cerr << "DUNE index set reports: " << leafIndexSet_->size( codim ) 367 << std::endl; 368 std::cerr << "ALBERTA mesh reports: " << mesh_.size( codim ) << std::endl; 369 consistent = false; 370 } 371 if( !consistent ) 372 abort(); 373 } 374 #endif 375 376 levelProvider_.markAllOld(); 377 adaptationState_.postAdapt(); 378 } 379 380 381 template< int dim, int dimworld > 382 inline bool AlbertaGrid< dim, dimworld > mark(int refCount,const typename Traits::template Codim<0>::Entity & e)383 ::mark( int refCount, const typename Traits::template Codim< 0 >::Entity &e ) 384 { 385 // if not leaf entity, leave method 386 if( !e.isLeaf() ) 387 return false; 388 389 // don't mark macro elements for coarsening 390 if( refCount < -e.level() ) 391 return false; 392 393 // take back previous marking 394 adaptationState_.unmark( getMark( e ) ); 395 396 // set new marking 397 adaptationState_.mark( refCount ); 398 e.impl().elementInfo().setMark( refCount ); 399 400 return true; 401 } 402 403 404 template< int dim, int dimworld > 405 inline int AlbertaGrid< dim, dimworld > getMark(const typename Traits::template Codim<0>::Entity & e) const406 ::getMark( const typename Traits::template Codim< 0 >::Entity &e ) const 407 { 408 return e.impl().elementInfo().getMark(); 409 } 410 411 412 template< int dim, int dimworld > adapt()413 inline bool AlbertaGrid< dim, dimworld >::adapt () 414 { 415 // this is already done in postAdapt 416 //levelProvider_.markAllOld(); 417 418 // adapt mesh 419 hIndexSet_.preAdapt(); 420 const bool refined = mesh_.refine(); 421 const bool coarsened = (adaptationState_.coarsen() ? mesh_.coarsen() : false); 422 adaptationState_.adapt(); 423 hIndexSet_.postAdapt(); 424 425 if( refined || coarsened ) 426 calcExtras(); 427 428 // return true if elements were created 429 return refined; 430 } 431 432 433 template< int dim, int dimworld > 434 template< class DataHandle > 435 inline bool AlbertaGrid < dim, dimworld > adapt(AdaptDataHandleInterface<This,DataHandle> & handle)436 ::adapt ( AdaptDataHandleInterface< This, DataHandle > &handle ) 437 { 438 preAdapt(); 439 440 typedef Alberta::AdaptRestrictProlongHandler 441 < This, AdaptDataHandleInterface< This, DataHandle > > 442 DataHandler; 443 DataHandler dataHandler( *this, handle ); 444 445 typedef AdaptationCallback< DataHandler > Callback; 446 typename Callback::DofVectorPointer callbackVector; 447 callbackVector.create( dofNumbering_.emptyDofSpace(), "Adaptation Callback" ); 448 callbackVector.template setupInterpolation< Callback >(); 449 callbackVector.template setupRestriction< Callback >(); 450 if( Callback::DofVectorPointer::supportsAdaptationData ) 451 callbackVector.setAdaptationData( &dataHandler ); 452 else 453 Alberta::adaptationDataHandler_ = &dataHandler; 454 455 bool refined = adapt(); 456 457 if( !Callback::DofVectorPointer::supportsAdaptationData ) 458 Alberta::adaptationDataHandler_ = 0; 459 callbackVector.release(); 460 461 postAdapt(); 462 return refined; 463 } 464 465 466 template< int dim, int dimworld > 467 inline const Alberta::GlobalVector & 468 AlbertaGrid< dim, dimworld > getCoord(const ElementInfo & elementInfo,int vertex) const469 ::getCoord ( const ElementInfo &elementInfo, int vertex ) const 470 { 471 assert( (vertex >= 0) && (vertex <= dim) ); 472 #if DUNE_ALBERTA_CACHE_COORDINATES 473 return coordCache_( elementInfo, vertex ); 474 #else 475 return elementInfo.coordinate( vertex ); 476 #endif 477 } 478 479 480 template < int dim, int dimworld > maxLevel() const481 inline int AlbertaGrid < dim, dimworld >::maxLevel() const 482 { 483 return maxlevel_; 484 } 485 486 487 template< int dim, int dimworld > size(int level,int codim) const488 inline int AlbertaGrid< dim, dimworld >::size ( int level, int codim ) const 489 { 490 return ((level >= 0) && (level <= maxlevel_) ? sizeCache_.size( level, codim ) : 0); 491 } 492 493 494 template< int dim, int dimworld > size(int level,GeometryType type) const495 inline int AlbertaGrid< dim, dimworld >::size ( int level, GeometryType type ) const 496 { 497 return ((level >= 0) && (level <= maxlevel_) ? sizeCache_.size( level, type ) : 0); 498 } 499 500 501 template< int dim, int dimworld > size(int codim) const502 inline int AlbertaGrid< dim, dimworld >::size ( int codim ) const 503 { 504 assert( sizeCache_.size( codim ) == mesh_.size( codim ) ); 505 return mesh_.size( codim ); 506 } 507 508 509 template< int dim, int dimworld > size(GeometryType type) const510 inline int AlbertaGrid< dim, dimworld >::size ( GeometryType type ) const 511 { 512 return sizeCache_.size( type ); 513 } 514 515 516 template < int dim, int dimworld > 517 inline const typename AlbertaGrid < dim, dimworld > :: Traits :: LevelIndexSet & levelIndexSet(int level) const518 AlbertaGrid < dim, dimworld > :: levelIndexSet (int level) const 519 { 520 // assert that given level is in range 521 assert( (level >= 0) && (level < (int)levelIndexVec_.size()) ); 522 523 if( levelIndexVec_[ level ] == 0 ) 524 { 525 levelIndexVec_[ level ] = new typename GridFamily::LevelIndexSetImp ( dofNumbering_ ); 526 levelIndexVec_[ level ]->update( lbegin< 0 >( level ), lend< 0 >( level ) ); 527 } 528 return *(levelIndexVec_[ level ]); 529 } 530 531 template < int dim, int dimworld > 532 inline const typename AlbertaGrid < dim, dimworld > :: Traits :: LeafIndexSet & leafIndexSet() const533 AlbertaGrid < dim, dimworld > :: leafIndexSet () const 534 { 535 if( leafIndexSet_ == 0 ) 536 { 537 leafIndexSet_ = new typename GridFamily::LeafIndexSetImp( dofNumbering_ ); 538 leafIndexSet_->update( leafbegin< 0 >(), leafend< 0 >() ); 539 } 540 return *leafIndexSet_; 541 } 542 543 544 template < int dim, int dimworld > calcExtras()545 inline void AlbertaGrid < dim, dimworld >::calcExtras () 546 { 547 // determine new maxlevel 548 maxlevel_ = levelProvider_.maxLevel(); 549 assert( (maxlevel_ >= 0) && (maxlevel_ < MAXL) ); 550 551 // unset up2Dat status, if lbegin is called then this status is updated 552 for( int l = 0; l < MAXL; ++l ) 553 levelMarkerVector_[ l ].clear(); 554 555 // unset up2Dat status, if leafbegin is called then this status is updated 556 leafMarkerVector_.clear(); 557 558 sizeCache_.reset(); 559 560 // update index sets (if they exist) 561 if( leafIndexSet_ != 0 ) 562 leafIndexSet_->update( leafbegin< 0 >(), leafend< 0 >() ); 563 for( unsigned int level = 0; level < levelIndexVec_.size(); ++level ) 564 { 565 if( levelIndexVec_[ level ] ) 566 levelIndexVec_[ level ]->update( lbegin< 0 >( level ), lend< 0 >( level ) ); 567 } 568 } 569 570 571 template< int dim, int dimworld > 572 inline bool AlbertaGrid< dim, dimworld > writeGrid(const std::string & filename,ctype time) const573 ::writeGrid ( const std::string &filename, ctype time ) const 574 { 575 if( filename.size() <= 0 ) 576 DUNE_THROW( AlbertaIOError, "No filename given to writeGridXdr." ); 577 return (mesh_.write( filename, time ) && hIndexSet_.write( filename )); 578 } 579 580 581 template< int dim, int dimworld > 582 inline bool AlbertaGrid< dim, dimworld > readGrid(const std::string & filename,ctype & time)583 ::readGrid ( const std::string &filename, ctype &time ) 584 { 585 //removeMesh(); 586 587 if( filename.size() <= 0 ) 588 return false; 589 590 numBoundarySegments_ = mesh_.read( filename, time ); 591 if( !mesh_ ) 592 DUNE_THROW( AlbertaIOError, "Could not read grid file: " << filename << "." ); 593 594 setup(); 595 hIndexSet_.read( filename ); 596 597 // calc maxlevel and indexOnLevel and so on 598 calcExtras(); 599 600 return true; 601 } 602 603 604 // AlbertaGrid::AdaptationCallback 605 // ------------------------------- 606 607 template< int dim, int dimworld > 608 template< class DataHandler > 609 struct AlbertaGrid< dim, dimworld >::AdaptationCallback 610 { 611 static const int dimension = dim; 612 613 typedef Alberta::DofVectorPointer< Alberta::GlobalVector > DofVectorPointer; 614 typedef Alberta::Patch< dimension > Patch; 615 getDataHandlerDune::AlbertaGrid::AdaptationCallback616 static DataHandler &getDataHandler ( const DofVectorPointer &dofVector ) 617 { 618 DataHandler *dataHandler; 619 if( DofVectorPointer::supportsAdaptationData ) 620 dataHandler = dofVector.getAdaptationData< DataHandler >(); 621 else 622 dataHandler = (DataHandler *)Alberta::adaptationDataHandler_; 623 assert( dataHandler != 0 ); 624 return *dataHandler; 625 } 626 interpolateVectorDune::AlbertaGrid::AdaptationCallback627 static void interpolateVector ( const DofVectorPointer &dofVector, 628 const Patch &patch ) 629 { 630 DataHandler &dataHandler = getDataHandler( dofVector ); 631 for( int i = 0; i < patch.count(); ++i ) 632 dataHandler.prolongLocal( patch, i ); 633 } 634 restrictVectorDune::AlbertaGrid::AdaptationCallback635 static void restrictVector ( const DofVectorPointer &dofVector, 636 const Patch &patch ) 637 { 638 DataHandler &dataHandler = getDataHandler( dofVector ); 639 for( int i = 0; i < patch.count(); ++i ) 640 dataHandler.restrictLocal( patch, i ); 641 } 642 }; 643 644 } // namespace Dune 645 646 #endif 647