dune-grid 3.0-git
albertagrid.cc
Go to the documentation of this file.
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
17namespace Dune
18{
19
20 namespace Alberta
21 {
23 }
24
25
26 template< int dim, int dimworld >
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 >
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 >
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 >
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 >
121 ::AlbertaGrid ( const std::string &macroGridFileName )
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 >
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 >
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 >
195 {
196 removeMesh();
197 }
198
199
200 template< int dim, int dimworld >
201 template< int codim, PartitionIteratorType pitype >
203 ::template Codim< codim >::template Partition<pitype>::LevelIterator
205 {
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 >
223 ::template Codim< codim >::template Partition< pitype >::LevelIterator
224 AlbertaGrid < dim, dimworld >::lend ( int level ) const
225 {
227 assert( level >= 0 );
228
229 return LevelIteratorImp( *this, level );
230 }
231
232
233 template< int dim, int dimworld >
234 template< int codim >
236 ::template Codim< codim >::LevelIterator
237 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 >
246 ::template Codim< codim >::LevelIterator
248 {
249 return lend< codim, All_Partition >( level );
250 }
251
252
253 template< int dim, int dimworld >
254 template< int codim, PartitionIteratorType pitype >
256 ::template Codim< codim >::template Partition< pitype >::LeafIterator
258 {
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 >
273 ::template Codim< codim >::template Partition< pitype >::LeafIterator
275 {
277 return LeafIteratorImp( *this, maxlevel_ );
278 }
279
280
281 template< int dim, int dimworld >
282 template< int codim >
284 ::template Codim< codim >::LeafIterator
286 {
287 return leafbegin< codim, All_Partition >();
288 }
289
290
291 template< int dim, int dimworld >
292 template< int codim >
294 ::template Codim< codim >::LeafIterator
295 AlbertaGrid < dim, dimworld >::leafend () const
296 {
297 return leafend< codim, All_Partition >();
298 }
299
300
301 template< int dim, int dimworld >
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();
320 }
321
322
323 template< int dim, int dimworld >
324 template< class DataHandle >
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 >
347 {
348 adaptationState_.preAdapt();
349 return adaptationState_.coarsen();
350 }
351
352
353 template < int dim, int dimworld >
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( 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
381 template< int dim, int dimworld >
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 getRealImplementation( e ).elementInfo().setMark( refCount );
399
400 return true;
401 }
402
403
404 template< int dim, int dimworld >
406 ::getMark( const typename Traits::template Codim< 0 >::Entity &e ) const
407 {
408 return getRealImplementation( e ).elementInfo().getMark();
409 }
410
411
412 template< int dim, int dimworld >
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();
425 if( refined || coarsened )
426 calcExtras();
427
428 // return true if elements were created
429 return refined;
430 }
432
433 template< int dim, int dimworld >
434 template< class DataHandle >
435 inline bool AlbertaGrid < dim, dimworld >
437 {
438 preAdapt();
439
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;
464
465
466 template< int dim, int dimworld >
467 inline const Alberta::GlobalVector &
469 ::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 >
481 inline int AlbertaGrid < dim, dimworld >::maxLevel() const
482 {
483 return maxlevel_;
484 }
485
486
487 template< int dim, int dimworld >
488 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 >
495 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 >
502 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 >
510 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 &
518 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 &
533 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 >
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 >
573 ::writeGridXdr ( const std::string &filename, ctype time ) const
574 {
575 return writeGrid( filename, time );
576 }
577
578
579 template< int dim, int dimworld >
581 ::readGridXdr ( const std::string &filename, ctype &time )
582 {
583 return readGrid( filename, time );
584 }
585
586
587 template< int dim, int dimworld >
589 ::writeGrid ( const std::string &filename, ctype time ) const
590 {
591 if( filename.size() <= 0 )
592 DUNE_THROW( AlbertaIOError, "No filename given to writeGridXdr." );
593 return (mesh_.write( filename, time ) && hIndexSet_.write( filename ));
594 }
595
596
597 template< int dim, int dimworld >
599 ::readGrid ( const std::string &filename, ctype &time )
600 {
601 //removeMesh();
602
603 if( filename.size() <= 0 )
604 return false;
605
606 numBoundarySegments_ = mesh_.read( filename, time );
607 if( !mesh_ )
608 DUNE_THROW( AlbertaIOError, "Could not read grid file: " << filename << "." );
609
610 setup();
611 hIndexSet_.read( filename );
612
613 // calc maxlevel and indexOnLevel and so on
614 calcExtras();
615
616 return true;
617 }
618
619
620 // AlbertaGrid::AdaptationCallback
621 // -------------------------------
622
623 template< int dim, int dimworld >
624 template< class DataHandler >
625 struct AlbertaGrid< dim, dimworld >::AdaptationCallback
626 {
627 static const int dimension = dim;
628
630 typedef Alberta::Patch< dimension > Patch;
631
632 static DataHandler &getDataHandler ( const DofVectorPointer &dofVector )
633 {
634 DataHandler *dataHandler;
635 if( DofVectorPointer::supportsAdaptationData )
636 dataHandler = dofVector.getAdaptationData< DataHandler >();
637 else
638 dataHandler = (DataHandler *)Alberta::adaptationDataHandler_;
639 assert( dataHandler != 0 );
640 return *dataHandler;
641 }
642
643 static void interpolateVector ( const DofVectorPointer &dofVector,
644 const Patch &patch )
645 {
646 DataHandler &dataHandler = getDataHandler( dofVector );
647 for( int i = 0; i < patch.count(); ++i )
648 dataHandler.prolongLocal( patch, i );
649 }
650
651 static void restrictVector ( const DofVectorPointer &dofVector,
652 const Patch &patch )
653 {
654 DataHandler &dataHandler = getDataHandler( dofVector );
655 for( int i = 0; i < patch.count(); ++i )
656 dataHandler.restrictLocal( patch, i );
657 }
658 };
659
660} // namespace Dune
661
662#endif
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 > &macroData)
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 misc.hh:29
Definition misc.hh:33
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