3#ifndef DUNE_ALBERTAGRID_CC
4#define DUNE_ALBERTAGRID_CC
26 template<
int dim,
int dimworld >
31 "Template Parameter dimworld does not match "
32 "ALBERTA's DIM_OF_WORLD setting.");
39 template<
int dim,
int dimworld >
40 inline AlbertaGrid < dim, dimworld >::AlbertaGrid ()
43 numBoundarySegments_( 0 ),
44 hIndexSet_( dofNumbering_ ),
46 levelIndexVec_( (size_t)MAXL, 0 ),
49 leafMarkerVector_( dofNumbering_ ),
50 levelMarkerVector_( (size_t)MAXL,
MarkerVector( dofNumbering_ ) )
52 checkAlbertaDimensions< dim, dimworld>();
56 template<
int dim,
int dimworld >
57 template<
class Proj,
class Impl >
63 numBoundarySegments_( 0 ),
64 hIndexSet_( dofNumbering_ ),
66 levelIndexVec_( (size_t)MAXL, 0 ),
69 leafMarkerVector_( dofNumbering_ ),
70 levelMarkerVector_( (size_t)MAXL,
MarkerVector( dofNumbering_ ) )
72 checkAlbertaDimensions< dim, dimworld >();
74 numBoundarySegments_ = mesh_.
create( macroData, projectionFactory );
76 DUNE_THROW(
AlbertaError,
"Invalid macro data structure." );
85 template<
int dim,
int dimworld >
91 numBoundarySegments_( 0 ),
92 hIndexSet_( dofNumbering_ ),
94 levelIndexVec_( (size_t)MAXL, 0 ),
97 leafMarkerVector_( dofNumbering_ ),
98 levelMarkerVector_( (size_t)MAXL,
MarkerVector( dofNumbering_ ) )
100 checkAlbertaDimensions< dim, dimworld >();
102 if( projection != 0 )
105 numBoundarySegments_ = mesh_.
create( macroData, projectionFactory );
108 numBoundarySegments_ = mesh_.
create( macroData );
110 DUNE_THROW(
AlbertaError,
"Invalid macro data structure." );
119 template <
int dim,
int dimworld >
124 hIndexSet_( dofNumbering_ ),
125 idSet_( hIndexSet_ ),
126 levelIndexVec_( (size_t)MAXL, 0 ),
129 leafMarkerVector_( dofNumbering_ ),
130 levelMarkerVector_( (size_t)MAXL,
MarkerVector( dofNumbering_ ) )
132 checkAlbertaDimensions< dim, dimworld >();
134 numBoundarySegments_ = mesh_.
create( macroGridFileName );
138 "Grid file '" << macroGridFileName
139 <<
"' is not in ALBERTA macro triangulation format." );
147 std::cout <<
typeName() <<
" created from macro grid file '"
148 << macroGridFileName <<
"'." << std::endl;
152 template<
int dim,
int dimworld >
155 dofNumbering_.create( mesh_ );
157 levelProvider_.create( dofNumbering_ );
159#if DUNE_ALBERTA_CACHE_COORDINATES
160 coordCache_.create( dofNumbering_ );
165 template<
int dim,
int dimworld >
166 inline void AlbertaGrid< dim, dimworld >::removeMesh ()
168 for(
size_t i = 0; i < levelIndexVec_.size(); ++i )
170 if( levelIndexVec_[ i ] != 0 )
171 delete levelIndexVec_[ i ];
172 levelIndexVec_[ i ] = 0;
175 if( leafIndexSet_ != 0 )
176 delete leafIndexSet_;
180 hIndexSet_.release();
181 levelProvider_.release();
182#if DUNE_ALBERTA_CACHE_COORDINATES
183 coordCache_.release();
185 dofNumbering_.release();
193 template<
int dim,
int dimworld >
200 template<
int dim,
int dimworld >
201 template<
int codim, PartitionIteratorType pitype >
203 ::template Codim< codim >::template Partition<pitype>::LevelIterator
207 assert( level >= 0 );
209 if( level > maxlevel_ )
210 return lend< codim, pitype >( level );
211 MarkerVector &markerVector = levelMarkerVector_[ level ];
213 if( (codim > 0) && !markerVector.
up2Date() )
214 markerVector.template markSubEntities< 1 >( lbegin< 0 >( level ), lend< 0 >( level ) );
216 return LevelIteratorImp( *
this, &markerVector, level );
220 template<
int dim,
int dimworld >
221 template<
int codim, PartitionIteratorType pitype >
223 ::template Codim< codim >::template Partition< pitype >::LevelIterator
224 AlbertaGrid < dim, dimworld >::lend (
int level )
const
227 assert( level >= 0 );
229 return LevelIteratorImp( *
this, level );
233 template<
int dim,
int dimworld >
234 template<
int codim >
236 ::template Codim< codim >::LevelIterator
237 AlbertaGrid < dim, dimworld >::lbegin (
int level )
const
239 return lbegin< codim, All_Partition >( level );
243 template<
int dim,
int dimworld >
244 template<
int codim >
246 ::template Codim< codim >::LevelIterator
249 return lend< codim, All_Partition >( level );
253 template<
int dim,
int dimworld >
254 template<
int codim, PartitionIteratorType pitype >
256 ::template Codim< codim >::template Partition< pitype >::LeafIterator
262 const int firstMarkedCodim = 2;
263 if( (codim >= firstMarkedCodim) && !markerVector.
up2Date() )
264 markerVector.template markSubEntities< firstMarkedCodim >( leafbegin< 0 >(), leafend< 0 >() );
266 return LeafIteratorImp( *
this, &markerVector, maxlevel_ );
270 template<
int dim,
int dimworld >
271 template<
int codim, PartitionIteratorType pitype >
273 ::template Codim< codim >::template Partition< pitype >::LeafIterator
277 return LeafIteratorImp( *
this, maxlevel_ );
281 template<
int dim,
int dimworld >
282 template<
int codim >
284 ::template Codim< codim >::LeafIterator
287 return leafbegin< codim, All_Partition >();
291 template<
int dim,
int dimworld >
292 template<
int codim >
294 ::template Codim< codim >::LeafIterator
295 AlbertaGrid < dim, dimworld >::leafend ()
const
297 return leafend< codim, All_Partition >();
301 template<
int dim,
int dimworld >
307 assert( (refCount >= 0) && (refCount + maxlevel_ < MAXL) );
309 for(
int i = 0; i < refCount; ++i )
312 const LeafIterator endit = leafend< 0 >();
313 for( LeafIterator it = leafbegin< 0 >(); it != endit; ++it )
323 template<
int dim,
int dimworld >
324 template<
class DataHandle >
331 assert( (refCount >= 0) && (refCount + maxlevel_ < MAXL) );
333 for(
int i = 0; i < refCount; ++i )
336 const LeafIterator endit = leafend< 0 >();
337 for( LeafIterator it = leafbegin< 0 >(); it != endit; ++it )
345 template<
int dim,
int dimworld >
349 return adaptationState_.coarsen();
353 template <
int dim,
int dimworld >
354 inline void AlbertaGrid < dim, dimworld >::postAdapt ()
357 if( leafIndexSet_ != 0 )
359 bool consistent =
true;
360 for(
int codim = 0; codim <= dimension; ++codim )
362 if( leafIndexSet_->size( codim ) == mesh_.size( codim ) )
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 )
368 std::cerr <<
"ALBERTA mesh reports: " << mesh_.size( codim ) << std::endl;
376 levelProvider_.markAllOld();
377 adaptationState_.postAdapt();
381 template<
int dim,
int dimworld >
390 if( refCount < -e.level() )
394 adaptationState_.unmark( getMark( e ) );
397 adaptationState_.
mark( refCount );
398 getRealImplementation( e ).elementInfo().setMark( refCount );
404 template<
int dim,
int dimworld >
408 return getRealImplementation( e ).elementInfo().
getMark();
412 template<
int dim,
int dimworld >
420 const bool refined = mesh_.refine();
421 const bool coarsened = (adaptationState_.coarsen() ? mesh_.coarsen() :
false);
422 adaptationState_.adapt();
423 hIndexSet_.postAdapt();
425 if( refined || coarsened )
433 template<
int dim,
int dimworld >
434 template<
class DataHandle >
435 inline bool AlbertaGrid < dim, dimworld >
443 DataHandler dataHandler( *
this, handle );
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 );
453 Alberta::adaptationDataHandler_ = &dataHandler;
455 bool refined = adapt();
457 if( !Callback::DofVectorPointer::supportsAdaptationData )
458 Alberta::adaptationDataHandler_ = 0;
459 callbackVector.release();
466 template<
int dim,
int dimworld >
471 assert( (vertex >= 0) && (vertex <= dim) );
472#if DUNE_ALBERTA_CACHE_COORDINATES
473 return coordCache_( elementInfo, vertex );
480 template <
int dim,
int dimworld >
481 inline int AlbertaGrid < dim, dimworld >::maxLevel()
const
487 template<
int dim,
int dimworld >
490 return ((level >= 0) && (level <= maxlevel_) ? sizeCache_.size( level, codim ) : 0);
494 template<
int dim,
int dimworld >
497 return ((level >= 0) && (level <= maxlevel_) ? sizeCache_.size( level, type ) : 0);
501 template<
int dim,
int dimworld >
504 assert( sizeCache_.size( codim ) == mesh_.size( codim ) );
505 return mesh_.size( codim );
509 template<
int dim,
int dimworld >
512 return sizeCache_.
size( type );
516 template <
int dim,
int dimworld >
517 inline const typename AlbertaGrid < dim, dimworld > :: Traits :: LevelIndexSet &
518 AlbertaGrid < dim, dimworld > :: levelIndexSet (
int level)
const
521 assert( (level >= 0) && (level < (
int)levelIndexVec_.size()) );
523 if( levelIndexVec_[ level ] == 0 )
526 levelIndexVec_[ level ]->update( lbegin< 0 >( level ), lend< 0 >( level ) );
528 return *(levelIndexVec_[ level ]);
531 template <
int dim,
int dimworld >
532 inline const typename AlbertaGrid < dim, dimworld > :: Traits :: LeafIndexSet &
533 AlbertaGrid < dim, dimworld > :: leafIndexSet ()
const
535 if( leafIndexSet_ == 0 )
538 leafIndexSet_->update( leafbegin< 0 >(), leafend< 0 >() );
540 return *leafIndexSet_;
544 template <
int dim,
int dimworld >
545 inline void AlbertaGrid < dim, dimworld >::calcExtras ()
548 maxlevel_ = levelProvider_.maxLevel();
549 assert( (maxlevel_ >= 0) && (maxlevel_ < MAXL) );
552 for(
int l = 0; l < MAXL; ++l )
553 levelMarkerVector_[ l ].clear();
556 leafMarkerVector_.clear();
561 if( leafIndexSet_ != 0 )
562 leafIndexSet_->update( leafbegin< 0 >(), leafend< 0 >() );
563 for(
unsigned int level = 0; level < levelIndexVec_.size(); ++level )
565 if( levelIndexVec_[ level ] )
566 levelIndexVec_[ level ]->update( lbegin< 0 >( level ), lend< 0 >( level ) );
571 template<
int dim,
int dimworld >
572 inline bool AlbertaGrid< dim, dimworld >
573 ::writeGridXdr (
const std::string &filename,
ctype time )
const
575 return writeGrid( filename, time );
579 template<
int dim,
int dimworld >
583 return readGrid( filename, time );
587 template<
int dim,
int dimworld >
591 if( filename.size() <= 0 )
592 DUNE_THROW(
AlbertaIOError,
"No filename given to writeGridXdr." );
593 return (mesh_.write( filename, time ) && hIndexSet_.write( filename ));
597 template<
int dim,
int dimworld >
603 if( filename.size() <= 0 )
606 numBoundarySegments_ = mesh_.read( filename, time );
608 DUNE_THROW(
AlbertaIOError,
"Could not read grid file: " << filename <<
"." );
611 hIndexSet_.read( filename );
623 template<
int dim,
int dimworld >
624 template<
class DataHandler >
625 struct AlbertaGrid< dim, dimworld >::AdaptationCallback
627 static const int dimension = dim;
632 static DataHandler &getDataHandler (
const DofVectorPointer &dofVector )
634 DataHandler *dataHandler;
635 if( DofVectorPointer::supportsAdaptationData )
636 dataHandler = dofVector.getAdaptationData< DataHandler >();
638 dataHandler = (DataHandler *)Alberta::adaptationDataHandler_;
639 assert( dataHandler != 0 );
643 static void interpolateVector (
const DofVectorPointer &dofVector,
646 DataHandler &dataHandler = getDataHandler( dofVector );
647 for(
int i = 0; i < patch.count(); ++i )
648 dataHandler.prolongLocal( patch, i );
651 static void restrictVector (
const DofVectorPointer &dofVector,
654 DataHandler &dataHandler = getDataHandler( dofVector );
655 for(
int i = 0; i < patch.count(); ++i )
656 dataHandler.restrictLocal( patch, i );
Include standard header files.
Definition agrid.hh:60
static void checkAlbertaDimensions()
Definition albertagrid.cc:27
static const int dimWorld
Definition misc.hh:43
static void * adaptationDataHandler_
Definition albertagrid.cc:22
ALBERTA REAL_D GlobalVector
Definition misc.hh:47
[ provides Dune::Grid ]
Definition agrid.hh:140
static std::string typeName()
Definition agrid.hh:443
GridFamily::ctype ctype
Definition agrid.hh:175
int size(int level, int codim) const
Number of grid entities per level and codim because lbegin and lend are none const,...
Definition albertagrid.cc:488
int getMark(const typename Traits::template Codim< 0 >::Entity &e) const
returns adaptation mark for given entity
Definition albertagrid.cc:406
bool preAdapt()
returns true, if a least one element is marked for coarsening
Definition albertagrid.cc:346
bool mark(int refCount, const typename Traits::template Codim< 0 >::Entity &e)
Marks an entity to be refined/coarsened in a subsequent adapt.
Definition albertagrid.cc:383
Definition albertagrid/datahandle.hh:25
unsigned int create(const MacroData< dim > ¯oData)
Definition meshpointer.hh:297
Definition dofvector.hh:177
const GlobalVector & coordinate(int vertex) const
Definition elementinfo.hh:683
void create()
Definition indexsets.cc:146
Definition albertagrid/indexsets.hh:334
Definition leafiterator.hh:21
Definition leveliterator.hh:21
Definition albertagrid/gridfamily.hh:97
Definition macrodata.hh:28
Definition albertagrid/projection.hh:78
Definition albertagrid/projection.hh:161
Definition refinement.hh:38
marker assigning subentities to one element containing them
Definition treeiterator.hh:30
bool up2Date() const
return true if marking is up to date
Definition treeiterator.hh:90
Interface class for the Grid's adapt method where the parameter is a AdaptDataHandleInterface.
Definition adaptcallback.hh:31
Interface class for vertex projection at the boundary.
Definition boundaryprojection.hh:24
A Traits struct that collects all associated types of one implementation.
Definition common/grid.hh:415