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> &macroData,
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> &macroData,
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 &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 >
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