dune-grid  2.2.0
sgrid.hh
Go to the documentation of this file.
00001 #ifndef DUNE_SGRID_HH
00002 #define DUNE_SGRID_HH
00003 
00004 #include <limits>
00005 #include <vector>
00006 #include <stack>
00007 
00008 #include <dune/common/fvector.hh>
00009 #include <dune/common/fmatrix.hh>
00010 #include <dune/common/bigunsignedint.hh>
00011 #include <dune/common/collectivecommunication.hh>
00012 #include <dune/common/reservedvector.hh>
00013 #include <dune/geometry/genericgeometry/topologytypes.hh>
00014 #include <dune/grid/common/capabilities.hh>
00015 #include <dune/grid/common/grid.hh>
00016 #include <dune/grid/sgrid/numbering.hh>
00017 #include <dune/grid/common/indexidset.hh>
00018 
00024 namespace Dune {
00025 
00026 //************************************************************************
00030   typedef double sgrid_ctype; 
00031 
00032   // globally define the persistent index type
00033   const int sgrid_dim_bits = 24;   // bits for encoding each dimension
00034   const int sgrid_level_bits = 6;  // bits for encoding level number
00035   const int sgrid_codim_bits = 4;  // bits for encoding codimension
00036 
00037 //************************************************************************
00038 // forward declaration of templates
00039 
00040 template<int dim, int dimworld, class GridImp> class SGeometry;
00041 template<int codim, int dim, class GridImp> class SEntity;
00042 template<int codim, class GridImp> class SEntityPointer;
00043 template<int codim, class GridImp> class SEntitySeed;
00044 template<int codim, PartitionIteratorType, class GridImp> class SLevelIterator;
00045 template<int dim, int dimworld, class ctype> class SGrid;
00046 template<class GridImp> class SIntersection;
00047 template<class GridImp> class SIntersectionIterator;
00048 template<class GridImp> class SHierarchicIterator;
00049 
00050 namespace FacadeOptions
00051 {
00052   
00053   template<int dim, int dimworld, class ctype, int mydim, int cdim>
00054   struct StoreGeometryReference<mydim, cdim,
00055                                 SGrid<dim,dimworld,ctype>, SGeometry>
00056   {
00057     static const bool v = false;
00058   };
00059   
00060   template<int dim, int dimworld, class ctype, int mydim, int cdim>
00061   struct StoreGeometryReference<mydim, cdim,
00062                                 const SGrid<dim,dimworld,ctype>, SGeometry>
00063   {
00064     static const bool v = false;
00065   };
00066   
00067 }
00068 
00069 //************************************************************************
00100 template<int mydim, int cdim, class GridImp> 
00101 class SGeometry
00102 : public GeometryDefaultImplementation<mydim,cdim,GridImp,SGeometry>
00103 {
00104 public:
00106         typedef typename GridImp::ctype ctype;
00107 
00109         GeometryType type () const
00110         {
00111           static const GeometryType cubeType(GeometryType::cube,mydim);
00112           return cubeType;
00113         }
00114 
00116         bool affine() const { return true ; }
00117 
00119         int corners () const
00120         {
00121           return 1<<mydim;
00122         }
00123 
00125         FieldVector< ctype, cdim > corner ( const int i ) const
00126         {
00127           return c[i];
00128         }
00129 
00131         FieldVector<ctype, cdim > center ( ) const
00132         {
00133           return centroid;
00134         }
00135 
00137         FieldVector<ctype, cdim> global (const FieldVector<ctype, mydim>& local) const;
00138 
00140         FieldVector<ctype, mydim> local (const FieldVector<ctype, cdim>& global) const;
00141 
00161         ctype integrationElement (const FieldVector<ctype, mydim>& local) const
00162         {
00163           return volume();
00164         }
00165 
00167         ctype volume() const;
00168         
00169         const FieldMatrix<ctype, mydim, cdim > &jacobianTransposed ( const FieldVector< ctype, mydim > &local ) const;
00170         const FieldMatrix<ctype,cdim,mydim>& jacobianInverseTransposed (const FieldVector<ctype, mydim>& local) const;
00171 
00173         void print (std::ostream& ss, int indent) const;
00174 
00179         void make (FieldMatrix<ctype,mydim+1,cdim>& __As);
00180 
00182         SGeometry () : builtinverse(false) {};
00183 
00184 private:
00185         FieldVector<ctype, cdim> s;          
00186         FieldVector<ctype, cdim> centroid;   
00187         FieldMatrix<ctype,mydim,cdim> A;     
00188         array<FieldVector<ctype, cdim>, 1<<mydim> c; 
00189         mutable FieldMatrix<ctype,cdim,mydim> Jinv;       
00190         mutable bool builtinverse;
00191 };
00192 
00194 template<int cdim, class GridImp> 
00195 class SGeometry<0,cdim,GridImp>
00196 : public GeometryDefaultImplementation<0,cdim,GridImp,SGeometry>
00197 {
00198 public:
00200         typedef typename GridImp::ctype ctype;
00201 
00203         GeometryType type () const
00204         {
00205           static const GeometryType cubeType(GeometryType::cube,0);
00206           return cubeType;
00207         }
00208 
00210         bool affine() const { return true ; }
00211 
00213         int corners () const
00214         {
00215           return 1;
00216         }
00217 
00219         FieldVector<ctype, cdim > corner ( const int i ) const
00220         {
00221           return s;
00222         }
00223 
00225         FieldVector<ctype, cdim > center ( ) const
00226         {
00227           return s;
00228         }
00229 
00231         void print (std::ostream& ss, int indent) const;
00232 
00234         void make (FieldMatrix<ctype,1,cdim>& __As);
00235 
00237         FieldVector<ctype, cdim> global (const FieldVector<ctype, 0>& local) const { return corner(0); }
00238 
00240         FieldVector<ctype, 0> local (const FieldVector<ctype, cdim>& global) const { return FieldVector<ctype,0> (0.0); } 
00241 
00263 
00264         SGeometry () {};
00265 
00271         ctype volume() const 
00272         {
00273           return 1;
00274         }
00275 
00281         ctype integrationElement(const FieldVector<ctype, 0>& local) const {
00282           return volume();
00283         }
00284         
00285         const FieldMatrix<ctype, 0, cdim > &jacobianTransposed ( const FieldVector<ctype, 0 > &local ) const
00286         {
00287           static const FieldMatrix<ctype, 0, cdim > dummy ( ctype( 0 ) );
00288           return dummy;
00289         }
00290 
00291         const FieldMatrix<ctype,cdim,0>& jacobianInverseTransposed (const FieldVector<ctype, 0>& local) const
00292         {
00293           static const FieldMatrix<ctype,cdim,0> dummy( ctype(0) );
00294           return dummy;
00295         }
00296 
00297 protected:
00298         FieldVector<ctype, cdim> s;         
00299 };
00300 
00301 template <int mydim, int cdim, class GridImp>
00302 inline std::ostream& operator<< (std::ostream& s, SGeometry<mydim,cdim,GridImp>& e)
00303 {
00304         e.print(s,0);
00305         return s;
00306 }
00307 
00308 //************************************************************************
00313 template<int codim, int dim, class GridImp, template<int,int,class> class EntityImp>
00314 class SEntityBase :
00315   public EntityDefaultImplementation<codim,dim,GridImp,EntityImp>
00316 {
00317   friend class SEntityPointer<codim,GridImp>;
00318   friend class SIntersectionIterator<GridImp>;
00319   enum { dimworld = GridImp::dimensionworld };
00320 
00321   typedef typename GridImp::Traits::template Codim< codim >::GeometryImpl GeometryImpl;
00322 
00323 public:
00324   typedef typename GridImp::ctype ctype;
00325   typedef typename GridImp::template Codim<codim>::Geometry Geometry;
00326   typedef typename GridImp::PersistentIndexType PersistentIndexType;
00327 
00329   int level () const
00330   {
00331     return l;
00332   }
00333 
00335   int globalIndex() const;
00336   
00340   SEntitySeed<codim, GridImp> seed () const {
00341     return SEntitySeed<codim, GridImp>(l, index);
00342   }
00343   
00345   GeometryType type () const
00346   {
00347     static const GeometryType cubeType(GeometryType::cube,dim-codim);
00348     return cubeType;
00349   }
00350 
00352   Geometry geometry () const
00353   {
00354       if (!builtgeometry) makegeometry();
00355       
00356       // return result
00357       return Geometry( geo );
00358   }
00359   
00360   PartitionType partitionType () const { return InteriorEntity; }
00361 
00363   SEntityBase (GridImp* _grid, int _l, int _index) :
00364     grid(_grid),
00365     l(_l),
00366     index(_index),
00367     z(grid->z(l,index,codim)),
00368     builtgeometry(false) {}
00369 
00371   SEntityBase () :
00372     builtgeometry(false) // mark geometry as not built
00373   {}
00374 
00376   SEntityBase ( const SEntityBase& other ) : 
00377     grid(other.grid), 
00378     l(other.l),
00379     index(other.index), 
00380     z(other.z), 
00381     geo(), // do not copy geometry 
00382     builtgeometry(false) // mark geometry as not built
00383   {}
00384     
00386   void make (GridImp* _grid, int _l, int _id);
00387 
00389   void make (int _l, int _id);
00390 
00392   void makegeometry () const;
00393   
00395   PersistentIndexType persistentIndex () const 
00396   {
00397     return grid->persistentIndex(l, codim, z);
00398   }
00399 
00401   int compressedIndex () const
00402   {
00403     return index;
00404   }
00405 
00407   int compressedLeafIndex () const
00408   {
00409     // codim != dim -> there are no copies of entities
00410     // maxlevel -> ids are fine
00411     if (codim<dim || l==grid->maxLevel())
00412       return compressedIndex();
00413     
00414     // this is a vertex which is not on the finest level
00415     // move coordinates up to maxlevel (multiply by 2 for each level
00416     array<int,dim> coord;
00417     for (int k=0; k<dim; k++)
00418       coord[k] = z[k]*(1<<(grid->maxLevel()-l));
00419     
00420     // compute number with respect to maxLevel
00421     return grid->n(grid->maxLevel(),coord);
00422   }
00423 
00424 protected:
00425   // this is how we implement our elements
00426   GridImp* grid;         
00427   int l;                 
00428   int index;             
00429   array<int,dim> z;      
00430   mutable GeometryImpl geo;  
00431   mutable bool builtgeometry;   
00432 };
00433 
00434 
00440 template<int codim, int dim, class GridImp> 
00441 class SEntity : public SEntityBase<codim,dim,GridImp,SEntity>
00442 {
00443   typedef Dune::SEntityBase<codim,dim,GridImp,Dune::SEntity> SEntityBase;
00444   friend class SEntityPointer<codim,GridImp>;
00445   friend class SIntersectionIterator<GridImp>;
00446 public:
00448   SEntity (GridImp* _grid, int _l, int _id) :
00449     SEntityBase(_grid,_l,_id) {};
00450 };
00451 
00478 template<int dim, class GridImp>
00479 class SEntity<0,dim,GridImp> : public SEntityBase<0,dim,GridImp,SEntity>
00480 {
00481   enum { dimworld = GridImp::dimensionworld };
00482   typedef Dune::SEntityBase<0,dim,GridImp,Dune::SEntity> SEntityBase;
00483   using SEntityBase::grid;
00484   using SEntityBase::l;
00485   using SEntityBase::index;
00486   using SEntityBase::z;
00487 
00488   typedef typename GridImp::Traits::template Codim< 0 >::GeometryImpl GeometryImpl;
00489   typedef typename GridImp::Traits::template Codim< 0 >::LocalGeometryImpl LocalGeometryImpl;
00490 
00491   friend class SEntityPointer<0,GridImp>;
00492   friend class SIntersectionIterator<GridImp>;
00493 
00494 public:
00495   typedef typename GridImp::ctype ctype;
00496   typedef typename GridImp::template Codim<0>::Geometry Geometry;
00497   typedef typename GridImp::template Codim<0>::LocalGeometry LocalGeometry;
00498   template <int cd>
00499   struct Codim
00500   {
00501     typedef typename GridImp::template Codim<cd>::EntityPointer EntityPointer;
00502   };
00503   typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
00504   typedef typename GridImp::LeafIntersectionIterator IntersectionIterator;
00505   typedef typename GridImp::HierarchicIterator HierarchicIterator;
00506   typedef typename GridImp::PersistentIndexType PersistentIndexType;
00507 
00509   friend class SHierarchicIterator<GridImp>;
00510   
00515   template<int cc> int count () const; 
00516 
00521   template<int cc> typename Codim<cc>::EntityPointer subEntity (int i) const;
00522 
00524   int subCompressedIndex (int codim, int i) const
00525   {
00526         if (codim==0) return this->compressedIndex();
00527     // compute subIndex
00528     return (this->grid)->n(this->l, this->grid->subz(this->z,i,codim));
00529   }
00530 
00534   int subCompressedLeafIndex (int codim, int i) const
00535   {
00536         if (codim==0) return this->compressedLeafIndex();
00537 
00538     assert(this->l == this->grid->maxLevel());
00539     // compute subIndex
00540     return (this->grid)->n(this->l, this->grid->subz(this->z,i,codim));
00541   }
00542 
00544   PersistentIndexType subPersistentIndex (int codim, int i) const
00545   {
00546         if (codim==0) return this->persistentIndex();
00547     // compute subId
00548     return this->grid->persistentIndex(this->l, codim, this->grid->subz(this->z,i,codim));
00549   }
00550   
00558   IntersectionIterator ibegin () const;
00559   IntersectionIterator ileafbegin () const;
00560   IntersectionIterator ilevelbegin () const;
00562   IntersectionIterator iend () const;
00563   IntersectionIterator ileafend () const;
00564   IntersectionIterator ilevelend () const;
00565 
00571   EntityPointer father () const;
00572 
00574   bool hasFather () const
00575   {
00576     return (this->level()>0);
00577   }
00578 
00580   bool isLeaf () const
00581     {
00582       return ( this->grid->maxLevel() == this->level() );
00583     }
00584 
00596   LocalGeometry geometryInFather () const;
00597 
00604   HierarchicIterator hbegin (int maxLevel) const;
00605 
00607   HierarchicIterator hend (int maxLevel) const;
00608 
00609   // members specific to SEntity
00611   SEntity (GridImp* _grid, int _l, int _index) : 
00612     SEntityBase(_grid,_l,_index),
00613     built_father(false)
00614     {}
00615 
00616   SEntity (const SEntity& other ) :
00617     SEntityBase(other.grid, other.l, other.index ),
00618     built_father(false)
00619     {}
00620 
00622   void make (GridImp* _grid, int _l, int _id)
00623     {
00624       SEntityBase::make(_grid,_l,_id);
00625       built_father = false;
00626     }
00627 
00629   void make (int _l, int _id)
00630     {
00631       SEntityBase::make(_l,_id);
00632       built_father = false;
00633     }
00634   
00635 private:
00636 
00637   SEntity();
00638   
00639   mutable bool built_father;
00640   mutable int father_index;
00641   mutable LocalGeometryImpl in_father_local;
00642   void make_father() const;
00643 };
00644 
00645 
00646 //************************************************************************
00654 struct SHierarchicStackElem {
00655   int l;
00656   int index;
00657   SHierarchicStackElem () : l(-1), index(-1) {}
00658   SHierarchicStackElem (int _l, int _index) {l=_l; index=_index;}
00659   bool operator== (const SHierarchicStackElem& s) const {return !operator!=(s);}
00660   bool operator!= (const SHierarchicStackElem& s) const {return l!=s.l || index!=s.index;}
00661 };
00662 
00663 template<class GridImp>
00664 class SHierarchicIterator :
00665   public Dune::SEntityPointer <0,GridImp>
00666 {
00667   friend class SHierarchicIterator<const GridImp>;
00668   enum { dim = GridImp::dimension };
00669   enum { dimworld = GridImp::dimensionworld };
00670   typedef Dune::SEntityPointer<0,GridImp> SEntityPointer;
00671   using SEntityPointer::realEntity;
00672   using SEntityPointer::grid;
00673   using SEntityPointer::l;
00674   using SEntityPointer::index;
00675 public:
00676   typedef typename GridImp::template Codim<0>::Entity Entity;
00677   typedef typename GridImp::ctype ctype;
00678 
00680   void increment();
00681 
00688   SHierarchicIterator (GridImp* _grid,
00689                        const Dune::SEntity<0,GridImp::dimension,GridImp>& _e,
00690                        int _maxLevel, bool makeend) :
00691     SEntityPointer(_grid,_e.level(),_e.compressedIndex())
00692     {
00693       // without sons, we are done
00694       // (the end iterator is equal to the calling iterator)
00695       if (makeend) return;
00696       
00697       // remember element where begin has been called
00698       orig_l = this->entity().level();
00699       orig_index = _grid->getRealImplementation(this->entity()).compressedIndex();
00700       
00701       // push original element on stack
00702       SHierarchicStackElem originalElement(orig_l, orig_index);
00703       stack.push(originalElement);
00704       
00705       // compute maxLevel
00706       maxLevel = std::min(_maxLevel,this->grid->maxLevel());
00707       
00708       // ok, push all the sons as well
00709       push_sons(orig_l,orig_index);
00710       
00711       // and pop the first son
00712       increment();
00713     }
00714 
00715 private:
00716   int maxLevel;                
00717   int orig_l, orig_index;         
00718 
00720   std::stack<SHierarchicStackElem, Dune::ReservedVector<SHierarchicStackElem,GridImp::MAXL> > stack;
00721   
00722   void push_sons (int level, int fatherid); 
00723 };
00724 
00725 //************************************************************************
00732 template<class GridImp>
00733 class SIntersectionIterator
00734 {
00735   enum { dim=GridImp::dimension };
00736   enum { dimworld=GridImp::dimensionworld };
00737 
00738   typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl;
00739   typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl;
00740 
00741 public:
00742   typedef typename GridImp::template Codim<0>::Entity Entity;
00743   typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
00744   typedef typename GridImp::template Codim<1>::Geometry Geometry;
00745   typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
00746   typedef Dune::SIntersection<GridImp> IntersectionImp;
00747   typedef Dune::Intersection<const GridImp, Dune::SIntersection> Intersection;
00749   enum { dimension=dim };
00751   enum { dimensionworld=dimworld };
00753   typedef typename GridImp::ctype ctype;
00754 
00756   bool equals(const SIntersectionIterator<GridImp>& i) const;
00758   void increment();
00759 
00761   const Intersection & dereference() const
00762   {
00763     return intersection;
00764   }
00765 
00768   EntityPointer inside() const;
00769 
00772   EntityPointer outside() const;
00773   
00775   bool boundary () const;
00776 
00778   bool conforming () const;
00779 
00780   int boundaryId () const {
00781     if (boundary()) return count + 1;
00782     return 0;
00783   };
00784 
00785   int boundarySegmentIndex () const {
00786     if (boundary())
00787       return grid->boundarySegmentIndex(self.level(), count, zred);
00788     return -1;
00789   };
00790   
00792   bool neighbor () const;
00793 
00795   FieldVector<ctype, GridImp::dimensionworld> outerNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
00796     {
00797       return unitOuterNormal(local);
00798     }
00800   FieldVector<ctype, GridImp::dimensionworld> unitOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
00801     {
00802       return centerUnitOuterNormal();
00803     }
00805   FieldVector<ctype, GridImp::dimensionworld> centerUnitOuterNormal () const
00806     {
00807       // while we are at it, compute normal direction
00808       FieldVector<ctype, dimworld> normal(0.0);
00809       if (count%2)
00810         normal[count/2] =  1.0; // odd
00811       else
00812         normal[count/2] = -1.0; // even
00813       
00814       return normal; 
00815     }
00817   FieldVector<ctype, GridImp::dimensionworld> integrationOuterNormal (const FieldVector<ctype, GridImp::dimension-1>& local) const
00818     {
00819         FieldVector<ctype, dimworld> n = unitOuterNormal(local);
00820         n *= geometry().integrationElement(local);
00821         return n;
00822     }
00823 
00827   LocalGeometry geometryInInside () const;
00831   LocalGeometry geometryInOutside () const;
00835   Geometry geometry () const;
00836 
00838   GeometryType type () const
00839   {
00840     static const GeometryType cubeType(GeometryType::cube,dim-1);
00841     return cubeType;
00842   }
00843 
00845   int indexInInside () const;
00847   int indexInOutside () const;
00848 
00850   SIntersectionIterator (GridImp* _grid, const SEntity<0,dim,GridImp >* _self, int _count) :
00851     self(*_self), ne(self), grid(_grid),
00852     partition(_grid->partition(grid->getRealImplementation(ne).l,_self->z)),
00853     zred(_grid->compress(grid->getRealImplementation(ne).l,_self->z)),
00854     intersection(IntersectionImp(*this))
00855   {
00856     // make neighbor
00857     make(_count);
00858   }
00859 
00860   SIntersectionIterator (const SIntersectionIterator & other) :
00861     self(other.self), ne(other.ne), grid(other.grid),
00862     partition(other.partition), zred(other.zred),
00863     count(other.count), valid_count(other.valid_count),
00864     valid_nb(other.valid_nb), is_on_boundary(other.is_on_boundary),
00865     built_intersections(false),
00866     intersection(IntersectionImp(*this))
00867   {
00868   }
00869   
00871   SIntersectionIterator& operator = (const SIntersectionIterator& other)
00872     {
00873       /* We can't assign the grid */
00874       assert(grid == other.grid);
00875       
00876       /* Assign data from other */
00877       self = other.self;
00878       ne = other.ne;
00879       partition = other.partition;
00880       zred = other.zred;
00881       count = other.count;
00882       valid_count = other.valid_count;
00883       valid_nb = other.valid_nb;
00884       is_on_boundary = other.is_on_boundary;
00885 
00886       /* mark cached data as invalid */
00887       built_intersections = false;
00888 
00889       return *this;
00890     }
00891   
00892 private:
00893   void make (int _count) const;           
00894   void makeintersections () const;        
00895   EntityPointer self;                     
00896   mutable EntityPointer ne;               
00897   const GridImp * grid;                   
00898   int partition;                          
00899   array<int,dim> zred;                    
00900   mutable int count;                      
00901   mutable bool valid_count;               
00902   mutable bool valid_nb;                  
00903   mutable bool is_on_boundary;            
00904   mutable bool built_intersections;       
00905   mutable LocalGeometryImpl is_self_local;    
00906   mutable GeometryImpl is_global;             
00907   mutable LocalGeometryImpl is_nb_local;      
00908   Intersection intersection;
00909 };
00910 
00911 template<class GridImp>
00912 class SIntersection
00913 {
00914   enum { dim=GridImp::dimension };
00915   enum { dimworld=GridImp::dimensionworld };
00916 public:
00917   typedef typename GridImp::template Codim<0>::Entity Entity;
00918   typedef typename GridImp::template Codim<0>::EntityPointer EntityPointer;
00919   typedef typename GridImp::template Codim<1>::Geometry Geometry;
00920   typedef typename Geometry::LocalCoordinate LocalCoordinate;
00921   typedef typename Geometry::GlobalCoordinate GlobalCoordinate;
00922   typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
00923   typedef Dune::Intersection<const GridImp, Dune::SIntersectionIterator> Intersection;
00925   enum { dimension=dim };
00927   enum { dimensionworld=dimworld };
00929   typedef typename GridImp::ctype ctype;
00930 
00931   bool boundary () const
00932   {
00933     return is.boundary();
00934   }
00935 
00937   int boundaryId () const
00938   {
00939     return is.boundaryId();
00940   }
00941 
00943   size_t boundarySegmentIndex () const
00944   {
00945     return is.boundarySegmentIndex();
00946   }
00947 
00949   bool neighbor () const 
00950   {
00951     return is.neighbor();
00952   }
00953 
00955   EntityPointer inside() const
00956   {
00957     return is.inside();
00958   }
00959 
00961   EntityPointer outside() const
00962   {
00963     return is.outside();
00964   }
00965   
00967   bool conforming () const 
00968   {
00969     return is.conforming();
00970   }
00971 
00973   LocalGeometry geometryInInside () const
00974   {
00975     return is.geometryInInside();
00976   }
00977 
00979   LocalGeometry geometryInOutside () const
00980   {
00981     return is.geometryInOutside();
00982   }
00983 
00985   Geometry geometry () const
00986   {
00987     return is.geometry();
00988   }
00989 
00991   GeometryType type () const
00992   {
00993     return is.type();
00994   }
00995 
00997   int indexInInside () const
00998   {
00999     return is.indexInInside();
01000   }
01001 
01003   int indexInOutside () const
01004   {
01005     return is.indexInOutside();
01006   }
01007 
01009   GlobalCoordinate outerNormal (const LocalCoordinate& local) const
01010   {
01011     return is.outerNormal(local);
01012   }
01013 
01015   GlobalCoordinate integrationOuterNormal (const LocalCoordinate& local) const
01016   {
01017     return is.integrationOuterNormal(local);
01018   }
01019 
01021   GlobalCoordinate unitOuterNormal (const LocalCoordinate& local) const
01022   {
01023     return is.unitOuterNormal(local);
01024   }
01025 
01027   GlobalCoordinate centerUnitOuterNormal () const
01028   {
01029     return is.centerUnitOuterNormal();
01030   }
01031 
01033   SIntersection (const SIntersectionIterator<GridImp> & is_) : is(is_) {}
01034   
01035 private:
01036 #ifndef DOXYGEN // doxygen can't handle this recursive usage
01037   const SIntersectionIterator<GridImp> & is;
01038 #endif
01039 };
01040 
01041 //************************************************************************
01042 
01046 template <class T>
01047 class AutoPtrStack : public std::stack<T*>
01048 {
01049 public:
01050     ~AutoPtrStack()
01051     {
01052         while(! this->empty())
01053         {
01054             T* e = this->top();
01055             delete e;
01056             this->pop();
01057         }
01058     }
01059 };
01060 
01063 template<int codim, class GridImp>
01064 class SEntityPointer
01065 {
01066   enum { dim = GridImp::dimension };
01067   friend class SIntersectionIterator<GridImp>;
01068 public:
01069   typedef SEntityPointer<codim,GridImp> EntityPointerImp;
01070   typedef typename GridImp::template Codim<codim>::Entity Entity;
01072   enum { codimension = codim };
01073   
01075   bool equals(const SEntityPointer<codim,GridImp>& i) const;
01077   Entity& dereference() const;
01079   int level () const;
01080 
01082   SEntityPointer (GridImp * _grid, int _l, int _index) :
01083     grid(_grid), l(_l), index(_index), 
01084     e(0) 
01085   {
01086   }
01087 
01089   SEntityPointer (const SEntity<codim,dim,GridImp> & _e) :
01090     grid(_e.grid), l(_e.l), index(_e.index),
01091     e(0)
01092   {
01093   }
01094 
01096   SEntityPointer (const SEntityPointer<codim,GridImp>& other) :
01097     grid(other.grid), l(other.l), index(other.index),
01098     e( 0 )
01099   {
01100   }
01101 
01103   ~SEntityPointer()
01104   {
01105     if( e )
01106       enStack().push( e );
01107 #ifndef NDEBUG 
01108     index = -1;
01109 #endif
01110   }
01111 
01113   SEntityPointer& operator = (const SEntityPointer& other)
01114   {
01115     grid = other.grid;
01116     l = other.l;
01117     index = other.index;
01118 
01119     // free current entity
01120     if( e )
01121       enStack().push( e );
01122     e = 0;
01123 
01124     return *this;
01125   }
01126 
01127 protected:
01128   inline SEntity<codim,dim,GridImp>& realEntity() const
01129   {
01130     return grid->getRealImplementation(entity());
01131   }
01132   
01133   inline Entity& entity() const
01134   {
01135     if( ! e )
01136     {
01137       e = getEntity( grid, l, index );
01138     }
01139     return *e;
01140   }
01141 
01142   typedef AutoPtrStack< Entity > EntityStackType;
01143   static inline EntityStackType& enStack()
01144   {
01145     static EntityStackType eStack;
01146     return eStack;
01147   }
01148 
01149   inline Entity* getEntity(GridImp* _grid, int _l, int _id ) const
01150   {
01151     // get stack reference 
01152     EntityStackType& enSt = enStack();
01153 
01154     if( enSt.empty() )
01155     {
01156       return (new Entity(SEntity<codim,dim,GridImp>(_grid, _l, _id)));
01157     }
01158     else
01159     {
01160       Entity* e = enSt.top();
01161       enSt.pop();
01162       grid->getRealImplementation(*e).make(_grid, _l,_id);
01163       return e;
01164     }
01165   }
01166 
01167   GridImp* grid;                 
01168   int l;                         
01169   mutable int index;             
01170   mutable Entity* e;             
01171 };
01172 
01175 template<int codim, class GridImp>
01176 class SEntitySeed
01177 {
01178   enum { dim = GridImp::dimension };
01179 public:
01180   enum { codimension = codim };
01181   
01183   SEntitySeed (int l, int index) :
01184     _l(l), _index(index)
01185   {}
01186 
01187   int level () const { return this->_l; }
01188   int index () const { return this->_index; }
01189 
01190 private:
01191   int _l;                         
01192   int _index;                     
01193 };
01194 
01195 //************************************************************************
01196 
01197 
01200 template<int codim, PartitionIteratorType pitype, class GridImp>
01201 class SLevelIterator :
01202   public Dune::SEntityPointer <codim,GridImp>
01203 {
01204   friend class SLevelIterator<codim, pitype,const GridImp>;
01205   enum { dim = GridImp::dimension };
01206   typedef Dune::SEntityPointer<codim,GridImp> SEntityPointer;
01207   using SEntityPointer::realEntity;
01208   using SEntityPointer::l;
01209   using SEntityPointer::index;
01210 public:
01211   typedef typename GridImp::template Codim<codim>::Entity Entity;
01212 
01214   void increment();
01215 
01217   SLevelIterator (GridImp * _grid, int _l, int _id) :
01218     SEntityPointer(_grid,_l,_id) {}
01219 };
01220 
01221 
01222 //========================================================================
01227 //========================================================================
01228 
01229 template<class GridImp>
01230 class SGridLevelIndexSet : public IndexSet<GridImp,SGridLevelIndexSet<GridImp> >
01231 {
01232   typedef SGridLevelIndexSet< GridImp > This;
01233   typedef IndexSet< GridImp, This > Base;
01234 
01235   enum { dim = GridImp::dimension };
01236 
01237 public:
01238 
01240   SGridLevelIndexSet ( const GridImp &g, int l )
01241   : grid( g ),
01242     level( l )
01243   {
01244     // TODO move list of geometrytypes to grid, can be computed static (singleton)
01245     // contains a single element type;
01246         for (int codim=0; codim<=GridImp::dimension; codim++)
01247           mytypes[codim].push_back(GeometryType(GeometryType::cube,GridImp::dimension-codim));
01248   }
01249 
01251   template<int cd>
01252   int index (const typename GridImp::Traits::template Codim<cd>::Entity& e) const 
01253   {
01254       return grid.getRealImplementation(e).compressedIndex(); 
01255   }
01256 
01257   template< int cc >
01258   int subIndex ( const typename GridImp::Traits::template Codim< cc >::Entity &e,
01259                  int i, unsigned int codim ) const
01260   {
01261       if( cc == 0 )
01262           return grid.getRealImplementation(e).subCompressedIndex(codim, i);
01263       else
01264           DUNE_THROW( NotImplemented, "subIndex for higher codimension entity not implemented for SGrid." );
01265   }
01266 
01267   // return true if the given entity is contained in \f$E\f$.
01268   template< class EntityType >
01269   bool contains ( const EntityType &e ) const
01270   {
01271     return (e.level() == level);
01272   }
01273 
01275   int size (GeometryType type) const
01276   {
01277     return grid.size( level, type );
01278   }
01279   
01281   int size (int codim) const
01282   {
01283     return grid.size( level, codim );
01284   }
01285 
01287   const std::vector<GeometryType>& geomTypes (int codim) const
01288   {
01289         return mytypes[codim];
01290   }
01291 
01292 private:
01293   const GridImp& grid;
01294   int level;
01295   std::vector<GeometryType> mytypes[GridImp::dimension+1];
01296 };
01297 
01298 // Leaf Index Set
01299 
01300 template<class GridImp>
01301 class SGridLeafIndexSet : public IndexSet<GridImp,SGridLeafIndexSet<GridImp> >
01302 {
01303   typedef SGridLeafIndexSet< GridImp > This;
01304   typedef IndexSet< GridImp, This > Base;
01305 
01306   enum { dim = GridImp::dimension };
01307 
01308 public:
01309 
01311   explicit SGridLeafIndexSet ( const GridImp &g )
01312   : grid( g )
01313   {
01314     // TODO move list of geometrytypes to grid, can be computed static (singleton)
01315     // contains a single element type;
01316         for (int codim=0; codim<=dim; codim++)
01317           mytypes[codim].push_back(GeometryType(GeometryType::cube,dim-codim));
01318   }
01319 
01321   /*
01322     We use the remove_const to extract the Type from the mutable class,
01323     because the const class is not instantiated yet.
01324   */
01325   template<int cd>
01326   int index (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const 
01327   {
01328         return grid.getRealImplementation(e).compressedLeafIndex(); 
01329   }
01330 
01331   template< int cc >
01332   int subIndex ( const typename GridImp::Traits::template Codim< cc >::Entity &e,
01333                  int i, unsigned int codim ) const
01334   {
01335       if( cc == 0 )
01336           return grid.getRealImplementation(e).subCompressedIndex(codim, i);
01337       else
01338           DUNE_THROW( NotImplemented, "subIndex for higher codimension entity not implemented for SGrid." );
01339   }
01340 
01342   int size (GeometryType type) const
01343   {
01344     return grid.size( grid.maxLevel(), type );
01345   }
01346 
01348   int size (int codim) const
01349   {
01350     return grid.size( grid.maxLevel(), codim );
01351   }
01352 
01353   // return true if the given entity is contained in \f$E\f$.
01354   template< class EntityType >
01355   bool contains ( const EntityType &e ) const
01356   {
01357     return (e.level() == grid.maxLevel());
01358   }
01359 
01361   const std::vector<GeometryType>& geomTypes (int codim) const
01362   {
01363         return mytypes[codim];
01364   }
01365 
01366 private:
01367   const GridImp& grid;
01368   std::vector<GeometryType> mytypes[dim+1];
01369 };
01370 
01371 
01372 //========================================================================
01377 //========================================================================
01378 
01379 template<class GridImp>
01380 class SGridGlobalIdSet :
01381   public IdSet<GridImp,SGridGlobalIdSet<GridImp>, typename remove_const<GridImp>::type::PersistentIndexType>
01382   /*
01383     We used the remove_const to extract the Type from the mutable class,
01384     because the const class is not instantiated yet.
01385   */
01386 {
01387   typedef SGridGlobalIdSet< GridImp > This;
01388 
01389 public:
01390   
01393   using IdSet<GridImp,SGridGlobalIdSet<GridImp>, typename remove_const<GridImp>::type::PersistentIndexType>::subId;
01394 
01396   /*
01397     We use the remove_const to extract the Type from the mutable class,
01398     because the const class is not instantiated yet.
01399   */
01400   typedef typename remove_const<GridImp>::type::PersistentIndexType IdType;
01401 
01403   explicit SGridGlobalIdSet ( const GridImp &g )
01404   : grid( g )
01405   {}
01406 
01408   /*
01409     We use the remove_const to extract the Type from the mutable class,
01410     because the const class is not instantiated yet.
01411   */
01412   template<int cd>
01413   IdType id (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const 
01414   {
01415           return grid.getRealImplementation(e).persistentIndex();
01416   }
01417 
01419   /*
01420     We use the remove_const to extract the Type from the mutable class,
01421     because the const class is not instantiated yet.
01422   */
01423   IdType subId ( const typename remove_const< GridImp >::type::Traits::template Codim< 0 >::Entity &e,
01424                  int i, unsigned int codim ) const
01425   {
01426       return grid.getRealImplementation(e).subPersistentIndex(codim, i);
01427   }
01428 
01429 private:
01430   const GridImp& grid;
01431 };
01432 
01433 
01434 template<int dim, int dimworld, class ctype>
01435 struct SGridFamily
01436 {
01437   typedef GridTraits<dim,dimworld,Dune::SGrid<dim,dimworld,ctype>,
01438                      SGeometry,SEntity,
01439                      SEntityPointer,SLevelIterator,
01440                      SIntersection, // leaf intersection
01441                      SIntersection, // level intersection
01442                      SIntersectionIterator, // leaf  intersection iter 
01443                      SIntersectionIterator, // level intersection iter
01444                      SHierarchicIterator,
01445                      SLevelIterator,
01446                                          SGridLevelIndexSet<const SGrid<dim,dimworld,ctype> >,
01447                                          SGridLeafIndexSet<const SGrid<dim,dimworld,ctype> >,
01448                                          SGridGlobalIdSet<const SGrid<dim,dimworld,ctype> >,
01449                                          bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>,
01450                                          SGridGlobalIdSet<const SGrid<dim,dimworld,ctype> >,
01451                                          bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits>, 
01452                                          CollectiveCommunication<Dune::SGrid<dim,dimworld,ctype> >,
01453                      DefaultLevelGridViewTraits, DefaultLeafGridViewTraits,
01454                      SEntitySeed>
01455   Traits;
01456 };
01457 
01458 
01459 //************************************************************************
01513 template<int dim, int dimworld, typename _ctype = sgrid_ctype>
01514 class SGrid : public GridDefaultImplementation <dim,dimworld,_ctype,SGridFamily<dim,dimworld,_ctype> >
01515 {
01516 public:
01517   typedef SGridFamily<dim,dimworld,_ctype> GridFamily;
01518   typedef bigunsignedint<dim*sgrid_dim_bits+sgrid_level_bits+sgrid_codim_bits> PersistentIndexType;
01519 
01520   // need for friend declarations in entity
01521   typedef SGridLevelIndexSet<SGrid<dim,dimworld> > LevelIndexSetType;
01522   typedef SGridLeafIndexSet<SGrid<dim,dimworld> > LeafIndexSetType;
01523   typedef SGridGlobalIdSet<SGrid<dim,dimworld> > GlobalIdSetType;
01524   
01525   typedef typename SGridFamily<dim,dimworld,_ctype>::Traits Traits;
01526 
01528   enum { MAXL=32 };
01529 
01531   typedef _ctype ctype;
01532 
01533   // constructors
01534 
01542   SGrid (const int * const N_, const ctype * const H_);
01543   
01551   SGrid (const int * const N_, const ctype * const L_, const ctype * const H_);
01552 
01562   SGrid (FieldVector<int,dim> N_, FieldVector<ctype,dim> L_, FieldVector<ctype,dim> H_);
01563 
01565   SGrid ();
01566 
01568   ~SGrid ();
01569 
01572   int maxLevel() const;
01573 
01575   template<int cd, PartitionIteratorType pitype>
01576   typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lbegin (int level) const;
01577 
01579   template<int cd, PartitionIteratorType pitype>
01580   typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator lend (int level) const;
01581 
01583   template<int cd>
01584   typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lbegin (int level) const
01585     {
01586       return lbegin<cd,All_Partition>(level);
01587     }
01588 
01590   template<int cd>
01591   typename Traits::template Codim<cd>::template Partition<All_Partition>::LevelIterator lend (int level) const
01592     {
01593       return lend<cd,All_Partition>(level);
01594     }
01595 
01597   template<int cd, PartitionIteratorType pitype>
01598   typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafbegin () const;
01599 
01601   template<int cd, PartitionIteratorType pitype>
01602   typename Traits::template Codim<cd>::template Partition<pitype>::LeafIterator leafend () const;
01603 
01605   template<int cd>
01606   typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafbegin () const
01607     {
01608       return leafbegin<cd,All_Partition>();
01609     }
01610 
01612   template<int cd>
01613   typename Traits::template Codim<cd>::template Partition<All_Partition>::LeafIterator leafend () const
01614     {
01615       return leafend<cd,All_Partition>();
01616     }
01617 
01618   // \brief obtain EntityPointer from EntitySeed. */
01619   template <typename Seed>
01620   typename Traits::template Codim<Seed::codimension>::EntityPointer
01621   entityPointer(const Seed& seed) const
01622   {
01623     enum { codim = Seed::codimension };
01624     return SEntityPointer<codim,const SGrid<dim,dimworld> >(this,seed.level(),seed.index());
01625   }
01626 
01640   template<class T, template<class> class P, int codim>
01641   void communicate (T& t, InterfaceType iftype, CommunicationDirection dir, int level)
01642   {
01643           // SGrid is sequential and has no periodic boundaries, so do nothing ...
01644           return;
01645   }
01646 
01648   int size (int level, int codim) const;
01649 
01651   int size (int codim) const
01652   {
01653         return size(maxLevel(),codim);
01654   }
01655 
01657   int size (int level, GeometryType type) const
01658   {
01659       return (type.isCube()) ? size(level,dim-type.dim()) : 0;
01660   }
01661 
01663   int size (GeometryType type) const
01664   {
01665         return size(maxLevel(),type);
01666   }
01667 
01669   size_t numBoundarySegments () const
01670   {
01671     return boundarysize;
01672   }
01673 
01675   int global_size (int codim) const;
01676 
01678   int overlapSize (int level, int codim)
01679   {
01680         return 0;
01681   }
01682 
01684   int ghostSize (int level, int codim)
01685   {
01686         return 0;
01687   }
01688 
01689   // these are all members specific to sgrid
01690 
01692   void globalRefine (int refCount);
01693 
01695     const array<int, dim>& dims(int level) const {
01696         return N[level];
01697     }
01698 
01700     const FieldVector<ctype, dimworld>& lowerLeft() const {
01701         return low;
01702     }
01703 
01705     FieldVector<ctype, dimworld> upperRight() const {
01706         return H;
01707     }
01708 
01710   bool adapt ()    
01711   {
01712         globalRefine(1);
01713         return true; 
01714   }
01715   
01716   // The new index sets from DDM 11.07.2005
01717   const typename Traits::GlobalIdSet& globalIdSet() const
01718   {
01719         return *theglobalidset;
01720   }
01721   
01722   const typename Traits::LocalIdSet& localIdSet() const
01723   {
01724         return *theglobalidset;
01725   }
01726 
01727   const typename Traits::LevelIndexSet& levelIndexSet(int level) const
01728   {
01729         assert(level>=0 && level<=maxLevel());
01730         return *(indexsets[level]);
01731   }
01732 
01733   const typename Traits::LeafIndexSet& leafIndexSet() const
01734   {
01735         return *theleafindexset;
01736   }
01737 
01742   template<class DataHandle>
01743   void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int level) const
01744   {
01745   }
01746 
01747   template<class DataHandle>
01748   void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
01749   {
01750   }
01751 
01752   const CollectiveCommunication<SGrid>& comm () const
01753   {
01754         return ccobj;
01755   }
01756 
01758   int overlapSize (int level, int codim) const
01759   {
01760         return 0;
01761   }
01762 
01764   int overlapSize (int codim) const
01765   {
01766         return 0;
01767   }
01768 
01770   int ghostSize (int level, int codim) const
01771   {
01772         return 0;
01773   }
01774 
01776   int ghostSize (int codim) const
01777   {
01778         return 0;
01779   }
01780 
01781   /*
01782     @}
01783    */
01784 
01785 private:
01786   /*
01787     Make associated classes friends to grant access to the real entity
01788   */
01789   friend class Dune::SGridLevelIndexSet<Dune::SGrid<dim,dimworld> >;
01790   friend class Dune::SGridLeafIndexSet<Dune::SGrid<dim,dimworld> >;
01791   friend class Dune::SGridGlobalIdSet<Dune::SGrid<dim,dimworld> >;
01792   friend class Dune::SIntersectionIterator<Dune::SGrid<dim,dimworld> >;
01793   friend class Dune::SHierarchicIterator<Dune::SGrid<dim,dimworld> >;
01794   friend class Dune::SEntity<0,dim,Dune::SGrid<dim,dimworld> >;
01795   
01796   friend class Dune::SGridLevelIndexSet<const Dune::SGrid<dim,dimworld> >;
01797   friend class Dune::SGridLeafIndexSet<const Dune::SGrid<dim,dimworld> >;
01798   friend class Dune::SGridGlobalIdSet<const Dune::SGrid<dim,dimworld> >;
01799   friend class Dune::SIntersectionIterator<const Dune::SGrid<dim,dimworld> >;
01800   friend class Dune::SHierarchicIterator<const Dune::SGrid<dim,dimworld> >;
01801   friend class Dune::SEntity<0,dim,const Dune::SGrid<dim,dimworld> >;
01802 
01803   template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
01804   friend class Dune::SEntityBase;
01805   
01806   template<int codim_, class GridImp_>
01807   friend class Dune::SEntityPointer;
01808   
01809   template<int codim_, int dim_, class GridImp_, template<int,int,class> class EntityImp_>
01810   friend class Entity;
01811 
01813   FieldVector<ctype, dimworld> pos (int level, array<int,dim>& z) const;
01814  
01816   int calc_codim (int level, const array<int,dim>& z) const;
01817 
01819   int n (int level, const array<int,dim>& z) const;
01820 
01822   array<int,dim> z (int level, int i, int codim) const;
01823         
01825   array<int,dim> subz (const array<int,dim> & z, int i, int codim) const;
01826   
01828   array<int,dim> compress (int level, const array<int,dim>& z) const; 
01829 
01831   array<int,dim> expand (int level, const array<int,dim>& r, int b) const; 
01832 
01836   int partition (int level, const array<int,dim>& z) const; 
01837 
01839   bool exists (int level, const array<int,dim>& zred) const;
01840 
01841   // compute boundary segment index for a given zentity and a face
01842   int boundarySegmentIndex (int l, int face, const array<int,dim> & zentity) const
01843   {
01844     array<int,dim-1> zface;
01845     int dir = face/2;
01846     int side = face%2;
01847     // compute z inside the global face
01848     for (int i=0; i<dir; i++) zface[i] = zentity[i]/(1<<l);
01849     for (int i=dir+1; i<dim; i++) zface[i-1] = zentity[i]/(1<<l);
01850     zface = boundarymapper[dir].expand(zface, 0);
01851     // compute index in the face
01852     int index = boundarymapper[dir].n(zface);
01853     // compute offset
01854     for (int i=0; i<dir; i++)
01855       index += 2*boundarymapper[i].elements(0);
01856     index += side*boundarymapper[dir].elements(0);
01857     return index;
01858   }
01859   
01860   // compute persistent index for a given zentity
01861   PersistentIndexType persistentIndex (int l, int codim, const array<int,dim> & zentity) const
01862   {
01863           if (codim!=dim)
01864                 {
01865                   // encode codim, this would actually not be necessary
01866           // because z is unique in codim
01867                   PersistentIndexType id(codim);
01868 
01869                   // encode level
01870                   id = id << sgrid_level_bits;
01871                   id = id+PersistentIndexType(l);
01872         
01873                   // encode coordinates
01874                   for (int i=dim-1; i>=0; i--)
01875                         {
01876                           id = id << sgrid_dim_bits;
01877                           id = id+PersistentIndexType(zentity[i]);
01878                         }
01879                   
01880                   return id;
01881                 }
01882           else
01883                 {
01884                   // determine min number of trailing zeroes
01885                   // consider that z is on the doubled grid !
01886                   int trailing = 1000;
01887                   for (int i=0; i<dim; i++)
01888                         {
01889                           // count trailing zeros
01890                           int zeros = 0;
01891                           for (int j=0; j<l; j++)
01892                                 if (zentity[i]&(1<<(j+1)))
01893                                   break;
01894                                 else
01895                                   zeros++;
01896                           trailing = std::min(trailing,zeros);
01897                         }
01898 
01899                   // determine the level of this vertex
01900                   int level = l-trailing;
01901 
01902                   // encode codim
01903                   PersistentIndexType id(dim);
01904 
01905                   // encode level
01906                   id = id << sgrid_level_bits;
01907                   id = id+PersistentIndexType(level);
01908         
01909                   // encode coordinates
01910                   for (int i=dim-1; i>=0; i--)
01911                         {
01912                           id = id << sgrid_dim_bits;
01913                           id = id+PersistentIndexType(zentity[i]>>trailing);
01914                         }
01915         
01916                   return id;
01917                 }
01918   }
01919     
01920   // disable copy and assign
01921   SGrid(const SGrid &) {};
01922   SGrid & operator = (const SGrid &) { return *this; };
01923   // generate SGrid 
01924   void makeSGrid (const array<int,dim>& N_, const FieldVector<ctype, dim>& L_, const FieldVector<ctype, dim>& H_);
01925 
01926   /*
01927     internal data
01928    */
01929   CollectiveCommunication<SGrid> ccobj;
01930 
01931   ReservedVector<SGridLevelIndexSet<const SGrid<dim,dimworld> >*, MAXL> indexsets;
01932   SGridLeafIndexSet<const SGrid<dim,dimworld> > *theleafindexset;
01933   SGridGlobalIdSet<const SGrid<dim,dimworld> > *theglobalidset;
01934 
01935   int L;                          // number of levels in hierarchic mesh 0<=level<L
01936   FieldVector<ctype, dim> low;    // lower left corner of the grid
01937   FieldVector<ctype, dim> H;      // length of cube per direction
01938   array<int,dim> *N;              // number of elements per direction
01939   FieldVector<ctype, dim> *h;     // mesh size per direction
01940   mutable CubeMapper<dim> *mapper;// a mapper for each level
01941 
01942   // boundary segement index set
01943   array<CubeMapper<dim-1>, dim> boundarymapper;  // a mapper for each coarse grid face
01944   int boundarysize;
01945 
01946   // faster implementation of subIndex 
01947   mutable array <int,dim> zrefStatic;     // for subIndex of SEntity
01948   mutable array <int,dim> zentityStatic;  // for subIndex of SEntity
01949 };
01950 
01951 namespace Capabilities
01952 {
01953 
01965   template<int dim, int dimw>
01966   struct hasSingleGeometryType< SGrid<dim,dimw> >
01967   {
01968     static const bool v = true;
01969     static const unsigned int topologyId = GenericGeometry :: CubeTopology< dim > :: type :: id ;
01970   };
01971   
01975   template<int dim, int dimw>
01976   struct isCartesian< SGrid<dim,dimw> >
01977   {
01978     static const bool v = true;
01979   };
01980   
01984   template<int dim, int dimw, int cdim>
01985   struct hasEntity< SGrid<dim,dimw>, cdim>
01986   {
01987     static const bool v = true;
01988   };
01989   
01993   template<int dim, int dimw>
01994   struct isLevelwiseConforming< SGrid<dim,dimw> >
01995   {
01996     static const bool v = true;
01997   };
01998 
02002   template<int dim, int dimw>
02003   struct isLeafwiseConforming< SGrid<dim,dimw> >
02004   {
02005     static const bool v = true;
02006   };
02007 
02008 } // end namespace Capabilities
02009 
02010 } // end namespace Dune
02011 
02012 #include"sgrid/sgrid.cc"
02013 
02014 #endif