dune-grid  2.3.1
agrid.hh
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_IMP_HH
4 #define DUNE_ALBERTAGRID_IMP_HH
5 
11 #if HAVE_ALBERTA || DOXYGEN
12 
13 #include <iostream>
14 #include <fstream>
15 
16 #include <algorithm>
17 #include <cassert>
18 #include <vector>
19 
20 // Dune includes
21 #include <dune/common/fvector.hh>
22 #include <dune/common/fmatrix.hh>
23 #include <dune/common/stdstreams.hh>
24 #include <dune/common/parallel/collectivecommunication.hh>
25 
26 #include <dune/grid/common/grid.hh>
30 
31 //- Local includes
32 // some cpp defines and include of alberta.h
33 #include "albertaheader.hh"
34 
35 // grape data io
37 
41 
49 
50 #include "indexsets.hh"
51 #include "geometry.hh"
52 #include "entity.hh"
53 #include "entitypointer.hh"
54 #include "hierarchiciterator.hh"
55 #include "treeiterator.hh"
56 #include "leveliterator.hh"
57 #include "leafiterator.hh"
58 
59 namespace Dune
60 {
61 
62  // External Forward Declarations
63  // -----------------------------
64 
65  template< class Grid >
66  struct DGFGridFactory;
67 
68 
69 
70  // AlbertaGrid
71  // -----------
72 
136  template< int dim, int dimworld = Alberta::dimWorld >
139  < dim, dimworld, Alberta::Real, AlbertaGridFamily< dim, dimworld > >
140  {
144  Base;
145 
146  template< int, int, class > friend class AlbertaGridEntity;
147  template< int, class > friend class AlbertaGridEntityPointer;
148  template< class, PartitionIteratorType > friend class AlbertaLevelGridView;
149  template< class, PartitionIteratorType > friend class AlbertaLeafGridView;
150 
151  friend class GridFactory< This >;
152  friend struct DGFGridFactory< This >;
153 
155 
156  friend class AlbertaGridIntersectionBase< const This >;
157  friend class AlbertaGridLeafIntersection< const This >;
158 
159  friend class AlbertaMarkerVector< dim, dimworld >;
160 #if (__GNUC__ < 4) && !(defined __ICC)
161  // add additional friend decls for gcc 3.4
162  friend struct AlbertaMarkerVector< dim, dimworld >::MarkSubEntities<true>;
163  friend struct AlbertaMarkerVector< dim, dimworld >::MarkSubEntities<false>;
164 #endif
165  friend class AlbertaGridIndexSet< dim, dimworld >;
166  friend class AlbertaGridHierarchicIndexSet< dim, dimworld >;
167 
168  template< class, class >
170 
171  public:
173  typedef AlbertaGridFamily< dim, dimworld > GridFamily;
174 
175  typedef typename GridFamily::ctype ctype;
176 
177  static const int dimension = GridFamily::dimension;
179 
180  // the Traits
181  typedef typename AlbertaGridFamily< dim, dimworld >::Traits Traits;
182 
187 
190 
194  typedef typename Traits::LocalIdSet LocalIdSet;
195 
198 
199  private:
201  typedef typename Traits::template Codim<0>::LeafIterator LeafIterator;
202 
204  typedef AlbertaGridIdSet<dim,dimworld> IdSetImp;
205 
207  struct AdaptationState
208  {
209  enum Phase { ComputationPhase, PreAdaptationPhase, PostAdaptationPhase };
210 
211  private:
212  Phase phase_;
213  int coarsenMarked_;
214  int refineMarked_;
215 
216  public:
217  AdaptationState ()
218  : phase_( ComputationPhase ),
219  coarsenMarked_( 0 ),
220  refineMarked_( 0 )
221  {}
222 
223  void mark ( int count )
224  {
225  if( count < 0 )
226  ++coarsenMarked_;
227  if( count > 0 )
228  refineMarked_ += (2 << count);
229  }
230 
231  void unmark ( int count )
232  {
233  if( count < 0 )
234  --coarsenMarked_;
235  if( count > 0 )
236  refineMarked_ -= (2 << count);
237  }
238 
239  bool coarsen () const
240  {
241  return (coarsenMarked_ > 0);
242  }
243 
244  int refineMarked () const
245  {
246  return refineMarked_;
247  }
248 
249  void preAdapt ()
250  {
251  if( phase_ != ComputationPhase )
252  error( "preAdapt may only be called in computation phase." );
253  phase_ = PreAdaptationPhase;
254  }
255 
256  void adapt ()
257  {
258  if( phase_ != PreAdaptationPhase )
259  error( "adapt may only be called in preadapdation phase." );
260  phase_ = PostAdaptationPhase;
261  }
262 
263  void postAdapt ()
264  {
265  if( phase_ != PostAdaptationPhase )
266  error( "postAdapt may only be called in postadaptation phase." );
267  phase_ = ComputationPhase;
268 
269  coarsenMarked_ = 0;
270  refineMarked_ = 0;
271  }
272 
273  private:
274  void error ( const std::string &message )
275  {
276  DUNE_THROW( InvalidStateException, message );
277  }
278  };
279 
280  template< class DataHandler >
281  struct AdaptationCallback;
282 
283  // max number of allowed levels is 64
284  static const int MAXL = 64;
285 
286  typedef Alberta::ElementInfo< dimension > ElementInfo;
287  typedef Alberta::MeshPointer< dimension > MeshPointer;
288  typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
289  typedef AlbertaGridLevelProvider< dimension > LevelProvider;
290 
291  // forbid copying and assignment
292  AlbertaGrid ( const This & );
293  This &operator= ( const This & );
294 
295  public:
297  AlbertaGrid ();
298 
304  AlbertaGrid ( const Alberta::MacroData< dimension > &macroData,
305  const Dune::shared_ptr< DuneBoundaryProjection< dimensionworld > > &projection
306  = Dune::shared_ptr< DuneBoundaryProjection< dimensionworld > >() );
307 
308  template< class Proj, class Impl >
309  AlbertaGrid ( const Alberta::MacroData< dimension > &macroData,
310  const Alberta::ProjectionFactoryInterface< Proj, Impl > &projectionFactory );
311 
316  AlbertaGrid ( const std::string &macroGridFileName );
317 
319  ~AlbertaGrid ();
320 
323  int maxLevel () const;
324 
326  template<int cd, PartitionIteratorType pitype>
327  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator
328  lbegin (int level) const;
329 
331  template<int cd, PartitionIteratorType pitype>
332  typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator
333  lend (int level) const;
334 
336  template< int codim >
337  typename Traits::template Codim< codim >::LevelIterator
338  lbegin ( int level ) const;
339 
341  template< int codim >
342  typename Traits::template Codim< codim >::LevelIterator
343  lend ( int level ) const;
344 
346  template< int codim, PartitionIteratorType pitype >
347  typename Traits
348  ::template Codim< codim >::template Partition< pitype >::LeafIterator
349  leafbegin () const;
350 
352  template< int codim, PartitionIteratorType pitype >
353  typename Traits
354  ::template Codim< codim >::template Partition< pitype >::LeafIterator
355  leafend () const;
356 
358  template< int codim >
359  typename Traits::template Codim< codim >::LeafIterator
360  leafbegin () const;
361 
363  template< int codim >
364  typename Traits::template Codim< codim >::LeafIterator
365  leafend () const;
366 
371  int size (int level, int codim) const;
372 
374  int size (int level, GeometryType type) const;
375 
377  int size (int codim) const;
378 
380  int size (GeometryType type) const;
381 
383  size_t numBoundarySegments () const
384  {
385  return numBoundarySegments_;
386  }
387 
389  template< PartitionIteratorType pitype >
390  typename Traits::template Partition< pitype >::LevelGridView
391  levelView ( int level ) const
392  {
393  typedef typename Traits::template Partition< pitype >::LevelGridView View;
394  typedef typename View::GridViewImp ViewImp;
395  return View( ViewImp( *this, level ) );
396  }
397 
399  template< PartitionIteratorType pitype >
400  typename Traits::template Partition< pitype >::LeafGridView leafView () const
401  {
402  typedef typename Traits::template Partition< pitype >::LeafGridView View;
403  typedef typename View::GridViewImp ViewImp;
404  return View( ViewImp( *this ) );
405  }
406 
408  typename Traits::template Partition< All_Partition >::LevelGridView
409  levelView ( int level ) const
410  {
411  typedef typename Traits::template Partition< All_Partition >::LevelGridView View;
412  typedef typename View::GridViewImp ViewImp;
413  return View( ViewImp( *this, level ) );
414  }
415 
417  typename Traits::template Partition< All_Partition >::LeafGridView leafView () const
418  {
419  typedef typename Traits::template Partition< All_Partition >::LeafGridView View;
420  typedef typename View::GridViewImp ViewImp;
421  return View( ViewImp( *this ) );
422  }
423 
424  public:
425  //***************************************************************
426  // Interface for Adaptation
427  //***************************************************************
428  using Base::getMark;
429  using Base::mark;
430 
432  int getMark ( const typename Traits::template Codim< 0 >::Entity &e ) const;
433 
435  bool mark ( int refCount, const typename Traits::template Codim< 0 >::Entity &e );
436 
438  void globalRefine ( int refCount );
439 
440  template< class DataHandle >
441  void globalRefine ( int refCount, AdaptDataHandleInterface< This, DataHandle > &handle );
442 
444  bool adapt ();
445 
447  template< class DataHandle >
449 
451  bool preAdapt ();
452 
454  void postAdapt();
455 
459  {
460  return comm_;
461  }
462 
463  static std::string typeName ()
464  {
465  std::ostringstream s;
466  s << "AlbertaGrid< " << dim << ", " << dimworld << " >";
467  return s.str();
468  }
469 
471  template< class EntitySeed >
472  typename Traits::template Codim< EntitySeed::codimension >::EntityPointer
473  entityPointer ( const EntitySeed &seed ) const
474  {
475  typedef typename Traits::template Codim< EntitySeed::codimension >::EntityPointerImpl EntityPointerImpl;
476  return EntityPointerImpl( *this, this->getRealImplementation(seed).elementInfo( meshPointer() ), this->getRealImplementation(seed).subEntity() );
477  }
478 
479  //**********************************************************
480  // End of Interface Methods
481  //**********************************************************
483  template< GrapeIOFileFormatType ftype >
484  bool writeGrid( const std::string &filename, ctype time ) const;
485 
487  template< GrapeIOFileFormatType ftype >
488  bool readGrid( const std::string &filename, ctype &time );
489 
490  // return hierarchic index set
491  const HierarchicIndexSet & hierarchicIndexSet () const { return hIndexSet_; }
492 
494  const typename Traits :: LevelIndexSet & levelIndexSet (int level) const;
495 
497  const typename Traits :: LeafIndexSet & leafIndexSet () const;
498 
500  const GlobalIdSet &globalIdSet () const
501  {
502  return idSet_;
503  }
504 
506  const LocalIdSet &localIdSet () const
507  {
508  return idSet_;
509  }
510 
511  // access to mesh pointer, needed by some methods
512  ALBERTA MESH* getMesh () const
513  {
514  return mesh_;
515  };
516 
517  const MeshPointer &meshPointer () const
518  {
519  return mesh_;
520  }
521 
522  const DofNumbering &dofNumbering () const
523  {
524  return dofNumbering_;
525  }
526 
527  const LevelProvider &levelProvider () const
528  {
529  return levelProvider_;
530  }
531 
532  int dune2alberta ( int codim, int i ) const
533  {
534  return numberingMap_.dune2alberta( codim, i );
535  }
536 
537  int alberta2dune ( int codim, int i ) const
538  {
539  return numberingMap_.alberta2dune( codim, i );
540  }
541 
542  int generic2alberta ( int codim, int i ) const
543  {
544  return genericNumberingMap_.dune2alberta( codim, i );
545  }
546 
547  int alberta2generic ( int codim, int i ) const
548  {
549  return genericNumberingMap_.alberta2dune( codim, i );
550  }
551 
552  // write ALBERTA mesh file
553  bool writeGridXdr ( const std::string &filename, ctype time ) const;
554 
556  bool readGridXdr ( const std::string &filename, ctype &time );
557 
558  private:
560 
561  typedef std::vector<int> ArrayType;
562 
563  void setup ();
564 
565  // make the calculation of indexOnLevel and so on.
566  // extra method because of Reihenfolge
567  void calcExtras();
568 
569  private:
570  // delete mesh and all vectors
571  void removeMesh();
572 
573  //***********************************************************************
574  // MemoryManagement for Entitys and Geometrys
575  //**********************************************************************
577  EntityObject;
578 
579  public:
581 
582  template< int codim >
583  static int
584  getTwist ( const typename Traits::template Codim< codim >::Entity &entity )
585  {
586  return getRealImplementation( entity ).twist();
587  }
588 
589  template< int codim >
590  static int
591  getTwist ( const typename Traits::template Codim< 0 >::Entity &entity, int subEntity )
592  {
593  return getRealImplementation( entity ).template twist< codim >( subEntity );
594  }
595 
596  static int
597  getTwistInInside ( const typename Traits::LeafIntersection &intersection )
598  {
599  return getRealImplementation( intersection ).twistInInside();
600  }
601 
602  static int
603  getTwistInOutside ( const typename Traits::LeafIntersection &intersection )
604  {
605  return getRealImplementation( intersection ).twistInOutside();
606  }
607 
609  getRealIntersection ( const typename Traits::LeafIntersection &intersection ) const
610  {
611  return getRealImplementation( intersection );
612  }
613 
614  public:
615  // read global element number from elNumbers_
616  const Alberta::GlobalVector &
617  getCoord ( const ElementInfo &elementInfo, int vertex ) const;
618 
619  private:
620  // pointer to an Albert Mesh, which contains the data
621  MeshPointer mesh_;
622 
623  // collective communication
625 
626  // maximum level of the mesh
627  int maxlevel_;
628 
629  // number of boundary segments within the macro grid
630  size_t numBoundarySegments_;
631 
632  // map between ALBERTA and DUNE numbering
635 
636  DofNumbering dofNumbering_;
637 
638  LevelProvider levelProvider_;
639 
640  // hierarchical numbering of AlbertaGrid, unique per codim
641  HierarchicIndexSet hIndexSet_;
642 
643  // the id set of this grid
644  IdSetImp idSet_;
645 
646  // the level index set, is generated from the HierarchicIndexSet
647  // is generated, when accessed
648  mutable std::vector< typename GridFamily::LevelIndexSetImp * > levelIndexVec_;
649 
650  // the leaf index set, is generated from the HierarchicIndexSet
651  // is generated, when accessed
652  mutable typename GridFamily::LeafIndexSetImp* leafIndexSet_;
653 
654  SizeCache< This > sizeCache_;
655 
656  typedef AlbertaMarkerVector< dim, dimworld > MarkerVector;
657 
658  // needed for VertexIterator, mark on which element a vertex is treated
659  mutable MarkerVector leafMarkerVector_;
660 
661  // needed for VertexIterator, mark on which element a vertex is treated
662  mutable std::vector< MarkerVector > levelMarkerVector_;
663 
664 #if DUNE_ALBERTA_CACHE_COORDINATES
666 #endif
667 
668  // current state of adaptation
669  AdaptationState adaptationState_;
670  };
671 
672 } // namespace Dune
673 
674 #include "albertagrid.cc"
675 
676 // undef all dangerous defines
677 #undef DIM
678 #undef DIM_OF_WORLD
679 
680 #ifdef _ABS_NOT_DEFINED_
681 #undef ABS
682 #endif
683 
684 #ifdef _MIN_NOT_DEFINED_
685 #undef MIN
686 #endif
687 
688 #ifdef _MAX_NOT_DEFINED_
689 #undef MAX
690 #endif
691 
692 #if DUNE_ALBERTA_VERSION >= 0x300
693 #ifdef obstack_chunk_alloc
694 #undef obstack_chunk_alloc
695 #endif
696 #ifdef obstack_chunk_free
697 #undef obstack_chunk_free
698 #endif
700 #else
702 #endif
703 
704 // We use MEM_ALLOC, so undefine it here.
705 #undef MEM_ALLOC
706 
707 // We use MEM_REALLOC, so undefine it here.
708 #undef MEM_REALLOC
709 
710 // We use MEM_CALLOC, so undefine it here.
711 #undef MEM_CALLOC
712 
713 // We use MEM_FREE, so undefine it here.
714 #undef MEM_FREE
715 
716 // Macro ERROR may be defined by alberta_util.h. If so, undefine it.
717 #ifdef ERROR
718 #undef ERROR
719 #endif // #ifdef ERROR
720 
721 // Macro ERROR_EXIT may be defined by alberta_util.h. If so, undefine it.
722 #ifdef ERROR_EXIT
723 #undef ERROR_EXIT
724 #endif // #ifdef ERROR_EXIT
725 
726 // Macro WARNING may be defined by alberta_util.h. If so, undefine it.
727 #ifdef WARNING
728 #undef WARNING
729 #endif // #ifdef WARNING
730 
731 // Macro TEST may be defined by alberta_util.h. If so, undefine it.
732 #ifdef TEST
733 #undef TEST
734 #endif // #ifdef TEST
735 
736 // Macro TEST_EXIT may be defined by alberta_util.h. If so, undefine it.
737 #ifdef TEST_EXIT
738 #undef TEST_EXIT
739 #endif // #ifdef TEST_EXIT
740 
741 // Macro DEBUG_TEST may be defined by alberta_util.h. If so, undefine it.
742 #ifdef DEBUG_TEST
743 #undef DEBUG_TEST
744 #endif // #ifdef DEBUG_TEST
745 
746 // Macro DEBUG_TEST_EXIT may be defined by alberta_util.h. If so, undefine it.
747 #ifdef DEBUG_TEST_EXIT
748 #undef DEBUG_TEST_EXIT
749 #endif // #ifdef DEBUG_TEST_EXIT
750 
751 // Macro INFO may be defined by alberta_util.h. If so, undefine it.
752 #ifdef INFO
753 #undef INFO
754 #endif // #ifdef INFO
755 
756 // Macro PRINT_INFO may be defined by alberta_util.h. If so, undefine it.
757 #ifdef PRINT_INFO
758 #undef PRINT_INFO
759 #endif // #ifdef PRINT_INFO
760 
761 // Macro PRINT_INT_VEC may be defined by alberta_util.h. If so, undefine it.
762 #ifdef PRINT_INT_VEC
763 #undef PRINT_INT_VEC
764 #endif // #ifdef PRINT_INT_VEC
765 
766 // Macro PRINT_REAL_VEC may be defined by alberta_util.h. If so, undefine it.
767 #ifdef PRINT_REAL_VEC
768 #undef PRINT_REAL_VEC
769 #endif // #ifdef PRINT_REAL_VEC
770 
771 // Macro WAIT may be defined by alberta_util.h. If so, undefine it.
772 #ifdef WAIT
773 #undef WAIT
774 #endif // #ifdef WAIT
775 
776 // Macro WAIT_REALLY may be defined by alberta_util.h. If so, undefine it.
777 #ifdef WAIT_REALLY
778 #undef WAIT_REALLY
779 #endif // #ifdef WAIT_REALLY
780 
781 // Macro GET_WORKSPACE may be defined by alberta_util.h. If so, undefine it.
782 #ifdef GET_WORKSPACE
783 #undef GET_WORKSPACE
784 #endif // #ifdef GET_WORKSPACE
785 
786 // Macro FREE_WORKSPACE may be defined by alberta_util.h. If so, undefine it.
787 #ifdef FREE_WORKSPACE
788 #undef FREE_WORKSPACE
789 #endif // #ifdef FREE_WORKSPACE
790 
791 // Macro MAT_ALLOC may be defined by alberta_util.h. If so, undefine it.
792 #ifdef MAT_ALLOC
793 #undef MAT_ALLOC
794 #endif // #ifdef MAT_ALLOC
795 
796 // Macro MAT_FREE may be defined by alberta_util.h. If so, undefine it.
797 #ifdef MAT_FREE
798 #undef MAT_FREE
799 #endif // #ifdef MAT_FREE
800 
801 // Macro NAME may be defined by alberta_util.h. If so, undefine it.
802 #ifdef NAME
803 #undef NAME
804 #endif // #ifdef NAME
805 
806 // Macro GET_STRUCT may be defined by alberta_util.h. If so, undefine it.
807 #ifdef GET_STRUCT
808 #undef GET_STRUCT
809 #endif // #ifdef GET_STRUCT
810 
811 // Macro ADD_PARAMETER may be defined by alberta_util.h. If so, undefine it.
812 #ifdef ADD_PARAMETER
813 #undef ADD_PARAMETER
814 #endif // #ifdef ADD_PARAMETER
815 
816 // Macro GET_PARAMETER may be defined by alberta_util.h. If so, undefine it.
817 #ifdef GET_PARAMETER
818 #undef GET_PARAMETER
819 #endif // #ifdef GET_PARAMETER
820 
821 #define _ALBERTA_H_
822 
823 #endif // HAVE_ALBERTA || DOXYGEN
824 
825 #endif