dune-grid  2.2.0
sizecache.hh
Go to the documentation of this file.
00001 #ifndef DUNE_SIZECACHE_HH
00002 #define DUNE_SIZECACHE_HH
00003 
00004 #include <cassert>
00005 #include <vector>
00006 #include <set>
00007 
00008 #include <dune/common/forloop.hh>
00009 #include <dune/common/exceptions.hh>
00010 
00011 #include <dune/geometry/type.hh>
00012 #include <dune/geometry/referenceelements.hh>
00013 
00014 #include <dune/grid/common/gridenums.hh>
00015 #include <dune/grid/common/capabilities.hh>
00016 
00023 namespace Dune {
00024 
00026 template <class GridImp> 
00027 class SizeCache 
00028 {
00029   typedef SizeCache<GridImp> ThisType;
00031   enum { dim    = GridImp::dimension   };
00032 
00034   enum { nCodim = GridImp::dimension+1 };
00035 
00036   // type of grid 
00037   typedef GridImp GridType;
00038 
00039   // coordinate type 
00040   typedef typename GridType :: ctype ctype ;
00041     
00042   // stores all sizes of the levels 
00043   mutable std::vector< int > levelSizes_[nCodim]; 
00044   
00045   // stores all sizes of the levels 
00046   mutable std::vector< std::vector< int > > levelTypeSizes_[nCodim]; 
00047   
00048   // stores all sizes of leafs 
00049   mutable int leafSizes_[nCodim]; 
00050 
00051   // stores all sizes of leafs 
00052   mutable std::vector< int > leafTypeSizes_[nCodim]; 
00053 
00054   // the grid 
00055   const GridType & grid_;
00056 
00057   // count elements of set by iterating the grid 
00058   template < int codim, bool gridHasCodim > 
00059   struct CountLevelEntitiesBase
00060   {
00061     template < class SzCacheType >
00062     static void apply(const SzCacheType & sc, int level, int cd)
00063     {   
00064       if( cd == codim )
00065       { 
00066         sc.template countLevelEntities<All_Partition,codim> (level);
00067       }
00068     }
00069   };
00070 
00071   template < int codim > 
00072   struct CountLevelEntitiesBase< codim, false >
00073   {
00074     template < class SzCacheType >
00075     static void apply(const SzCacheType & sc, int level, int cd)
00076     {   
00077       if( cd == codim ) 
00078       {
00079         sc.template countLevelEntitiesNoCodim<All_Partition,codim> (level);
00080       }
00081     }
00082   };
00083 
00084   template < int codim > 
00085   struct CountLevelEntities 
00086     : public CountLevelEntitiesBase< codim, Capabilities :: hasEntity< GridType, codim > :: v >
00087   {
00088   };
00089 
00090   // count elements of set by iterating the grid 
00091   template < int codim, bool gridHasCodim >
00092   struct CountLeafEntitiesBase
00093   {
00094     template <class SzCacheType>
00095     static void apply(const SzCacheType & sc, int cd)
00096     {   
00097       if( cd == codim )
00098       { 
00099         sc.template countLeafEntities<All_Partition,codim> ();
00100       }
00101     }
00102   };
00103 
00104   // count elements of set by iterating the grid 
00105   template < int codim > 
00106   struct CountLeafEntitiesBase< codim, false >
00107   {
00108     template <class SzCacheType>
00109     static void apply(const SzCacheType & sc, int cd)
00110     {
00111       if( cd == codim ) 
00112       {
00113         sc.template countLeafEntitiesNoCodim<All_Partition,codim> ();
00114       }
00115     }
00116   };
00117 
00118   template < int codim > 
00119   struct CountLeafEntities 
00120     : public CountLeafEntitiesBase< codim, Capabilities :: hasEntity< GridType, codim > :: v >
00121   {
00122   };
00123 
00124   int gtIndex( const GeometryType& type ) const 
00125   {
00126     return type.id() >> 1 ;
00127   }
00128 
00129   int sizeCodim( const int codim ) const 
00130   {
00131     const int mydim = GridType :: dimension - codim;
00132     return ((1 << mydim) + 1) / 2;
00133   }
00134   
00135   // private copy constructor 
00136   SizeCache (const SizeCache & );
00137 public:
00139   SizeCache (const GridType & grid) : grid_( grid ) 
00140   {
00141     reset();
00142   }
00143 
00145   void reset() 
00146   {
00147     for(int codim=0; codim<nCodim; ++codim) 
00148     {
00149       leafSizes_[ codim ] = -1;
00150       leafTypeSizes_[ codim ].resize( sizeCodim( codim ), -1 );
00151     }
00152     
00153     const int numMxl = grid_.maxLevel()+1;
00154     for(int codim=0; codim<nCodim; ++codim)
00155     {  
00156       std::vector<int> & vec = levelSizes_[codim]; 
00157       vec.resize(numMxl);
00158       levelTypeSizes_[codim].resize( numMxl );
00159       for(int level = 0; level<numMxl; ++level)
00160       {
00161         vec[level] = -1;
00162         levelTypeSizes_[codim][level].resize( sizeCodim( codim ), -1 );
00163       }
00164     }
00165   }
00166 
00167   //********************************************************************
00168   // level sizes 
00169   //********************************************************************
00171   int size (int level, int codim) const
00172   {
00173     assert( codim >= 0 );
00174     assert( codim < nCodim );
00175     assert( level >= 0 );
00176     if( level >= (int) levelSizes_[codim].size() ) return 0;
00177 
00178     if( levelSizes_[codim][level] < 0)
00179       ForLoop< CountLevelEntities, 0, dim > :: apply( *this, level, codim );
00180 
00181     //  CountLevelEntities<ThisType,All_Partition,dim>::count(*this,level,codim);
00182 
00183     assert( levelSizes_[codim][level] >= 0 );
00184     return levelSizes_[codim][level];
00185   }
00186   
00188   int size (int level, GeometryType type) const
00189   {
00190     const int codim = GridType ::dimension - type.dim();
00191     if( levelSizes_[codim][level] < 0)
00192       ForLoop< CountLevelEntities, 0, dim > :: apply( *this, level, codim );
00193 
00194     assert( levelTypeSizes_[codim][level][gtIndex( type )] >= 0 );
00195     return levelTypeSizes_[codim][level][gtIndex( type )];
00196   }
00197  
00198   //********************************************************************
00199   // leaf sizes 
00200   //********************************************************************
00202   int size (int codim) const
00203   {
00204     assert( codim >= 0 );
00205     assert( codim < nCodim );
00206     if( leafSizes_[codim] < 0 )
00207       ForLoop< CountLeafEntities, 0, dim > :: apply( *this, codim );
00208 
00209     assert( leafSizes_[codim] >= 0 );
00210     return leafSizes_[codim];
00211   }; 
00212   
00214   int size ( const GeometryType type ) const
00215   {
00216     const int codim = GridType :: dimension - type.dim();
00217     if( leafSizes_[codim] < 0 )
00218       ForLoop< CountLeafEntities, 0, dim > :: apply( *this, codim );
00219 
00220     assert( leafTypeSizes_[codim][ gtIndex( type )] >= 0 );
00221     return leafTypeSizes_[codim][ gtIndex( type )];
00222   }
00223 
00224 private:
00225   template <PartitionIteratorType pitype, int codim> 
00226   void countLevelEntities(int level) const 
00227   {
00228     typedef typename GridType :: LevelGridView GridView ;
00229     typedef typename GridView :: template Codim< codim > :: template Partition<pitype>  :: Iterator Iterator ;
00230     GridView gridView = grid_.levelView( level );
00231     Iterator it  = gridView.template begin<codim,pitype> ();
00232     Iterator end = gridView.template end<codim,pitype>   ();
00233     levelSizes_[codim][level] = countElements(it,end, levelTypeSizes_[codim][level]);
00234   }
00235   
00236   template <PartitionIteratorType pitype, int codim> 
00237   void countLeafEntities() const 
00238   {
00239     // count All_Partition entities 
00240     typedef typename GridType :: LeafGridView GridView ;
00241     typedef typename GridView :: template Codim< codim > :: template Partition<pitype>  :: Iterator Iterator ;
00242     GridView gridView = grid_.leafView();
00243     Iterator it  = gridView.template begin<codim,pitype> ();
00244     Iterator end = gridView.template end<codim,pitype>   ();
00245     leafSizes_[codim] = countElements(it,end, leafTypeSizes_[codim] );
00246   }
00247   
00248   // counts entities with given type for given iterator 
00249   template <class IteratorType> 
00250   int countElements(IteratorType & it, const IteratorType & end, std::vector<int>& typeSizes) const 
00251   {
00252     int overall = 0;
00253     const size_t types = typeSizes.size();
00254     for(size_t i=0; i<types; ++i) typeSizes[i] = 0;
00255     for( ; it != end; ++it )
00256     {
00257       const GeometryType type = it->type();
00258       ++typeSizes[ gtIndex( type ) ];
00259       ++overall;
00260     }
00261 
00262     int sumtypes = 0;
00263     for(size_t i=0; i<types; ++i) sumtypes += typeSizes[i];
00264 
00265     assert( overall == sumtypes );
00266     return overall;
00267   }
00268 
00269   template <PartitionIteratorType pitype, int codim> 
00270   void countLevelEntitiesNoCodim(int level) const 
00271   {
00272     typedef typename GridType :: LevelGridView GridView ;
00273     typedef typename GridView :: template Codim< 0 > :: template Partition<pitype>  :: Iterator Iterator ;
00274     GridView gridView = grid_.levelView( level );
00275     Iterator it  = gridView.template begin< 0, pitype> ();
00276     Iterator end = gridView.template end< 0, pitype>   ();
00277     levelSizes_[codim][level] = countElementsNoCodim< codim >(it,end, levelTypeSizes_[codim][level]);
00278   }
00279   
00280   template <PartitionIteratorType pitype, int codim> 
00281   void countLeafEntitiesNoCodim() const 
00282   {
00283     // count All_Partition entities 
00284     typedef typename GridType :: LeafGridView GridView ;
00285     typedef typename GridView :: template Codim< 0 > :: template Partition<pitype>  :: Iterator Iterator ;
00286     GridView gridView = grid_.leafView();
00287     Iterator it  = gridView.template begin< 0, pitype > ();
00288     Iterator end = gridView.template end< 0, pitype >   ();
00289     leafSizes_[codim] = countElementsNoCodim< codim >(it,end, leafTypeSizes_[codim] );
00290   }
00291   
00292   // counts entities with given type for given iterator 
00293   template < int codim, class IteratorType > 
00294   int countElementsNoCodim(IteratorType & it, const IteratorType & end, std::vector<int>& typeSizes) const 
00295   {
00296     typedef typename GridType :: LocalIdSet LocalIdSet ;
00297     typedef typename LocalIdSet :: IdType  IdType ;
00298 
00299     typedef GenericReferenceElement< ctype, dim > ReferenceElementType;
00300     typedef GenericReferenceElements< ctype, dim > ReferenceElementContainerType;
00301 
00302     typedef std::set< IdType > CodimIdSetType ;
00303 
00304     typedef typename IteratorType :: Entity  ElementType ;
00305 
00306     // get id set 
00307     const LocalIdSet& idSet = grid_.localIdSet();
00308 
00309     const size_t types = typeSizes.size();
00310     for(size_t i=0; i<types; ++i) typeSizes[ i ] = 0;
00311 
00312     std::vector< CodimIdSetType > typeCount( types );
00313 
00314     // count all elements of codimension codim 
00315     for( ; it != end; ++it )
00316     {
00317       // get entity 
00318       const ElementType& element = *it ;
00319       // get reference element 
00320       const ReferenceElementType& refElem = 
00321         ReferenceElementContainerType :: general( element.type() );
00322 
00323       // count all sub entities of codimension codim 
00324       const int count = element.template count< codim > ();
00325       for( int i=0; i< count; ++ i ) 
00326       {
00327         // get geometry type 
00328         const GeometryType geomType = refElem.type( i, codim );
00329         // get id of sub entity 
00330         const IdType id = idSet.subId( element, i, codim );
00331         // insert id into set 
00332         typeCount[ gtIndex( geomType ) ].insert( id );
00333       }
00334     }
00335 
00336     // accumulate numbers 
00337     int overall = 0;
00338     for(size_t i=0; i<types; ++i) 
00339     {
00340       typeSizes[ i ] = typeCount[ i ].size();
00341       overall += typeSizes[ i ];
00342     }
00343 
00344     return overall;
00345   }
00346 };
00347 
00349 template <class GridImp> 
00350 struct SingleTypeSizeCache : public SizeCache< GridImp > 
00351 {
00352   DUNE_DEPRECATED SingleTypeSizeCache( const GridImp& grid, bool isSimplex , bool isCube, bool notWorry = false ) : SizeCache< GridImp >( grid ) {}
00353 };
00354 
00355 } // end namespace Dune 
00356 #endif