dune-grid  2.2.0
gridptr.hh
Go to the documentation of this file.
00001 #ifndef DUNE_DGF_GRIDPTR_HH
00002 #define DUNE_DGF_GRIDPTR_HH
00003 
00004 #include <iostream>
00005 #include <string>
00006 #include <vector>
00007 #include <memory>
00008 #include <map>
00009 #include <assert.h>
00010 
00011 //- Dune includes 
00012 #include <dune/common/mpihelper.hh>
00013 #include <dune/grid/common/gridenums.hh>
00014 #include <dune/grid/common/datahandleif.hh>
00015 
00016 #include <dune/grid/io/file/dgfparser/dgfexception.hh>
00017 #include <dune/grid/io/file/dgfparser/entitykey.hh>
00018 #include <dune/grid/io/file/dgfparser/parser.hh>
00019 
00020 #include <dune/grid/common/intersection.hh>
00021 
00022 namespace Dune
00023 {
00024   // forward declarations
00025   // --------------------
00026   template < class G >
00027   struct DGFGridFactory;
00028 
00029   template< class GridImp, template < class > class IntersectionImp >
00030   class Intersection;
00031 
00044   template< class GridType >
00045   struct GridPtr
00046   {
00047     typedef MPIHelper::MPICommunicator MPICommunicatorType;
00048     static const int dimension = GridType::dimension;
00049 
00051     explicit GridPtr ( const std::string &filename,
00052                        MPICommunicatorType comm = MPIHelper::getCommunicator() )
00053     : gridPtr_( 0 ),
00054       elParam_(),
00055       vtxParam_(),
00056       bndParam_(),
00057       bndId_(),
00058       emptyParam_(),
00059       nofElParam_( 0 ),
00060       nofVtxParam_( 0 ),
00061       haveBndParam_( false )
00062     {
00063       DGFGridFactory< GridType > dgfFactory( filename, comm );
00064       initialize( dgfFactory );
00065     }
00066 
00068     explicit GridPtr ( std::istream &input,
00069                        MPICommunicatorType comm = MPIHelper::getCommunicator() )
00070     : gridPtr_( 0 ),
00071       elParam_(),
00072       vtxParam_(),
00073       bndParam_(),
00074       bndId_(),
00075       emptyParam_(),
00076       nofElParam_( 0 ),
00077       nofVtxParam_( 0 ),
00078       haveBndParam_( false )
00079     {
00080       DGFGridFactory< GridType > dgfFactory( input, comm );
00081       initialize( dgfFactory );
00082     }
00083 
00085     GridPtr()
00086     : gridPtr_(0),
00087       elParam_(),
00088       vtxParam_(),
00089       bndParam_(),
00090       bndId_(),
00091       emptyParam_(),
00092       nofElParam_(0),
00093       nofVtxParam_(0),
00094       haveBndParam_( false )
00095     {}
00096 
00098     explicit GridPtr( GridType *grd )
00099     : gridPtr_(grd),
00100       elParam_(),
00101       vtxParam_(),
00102       bndParam_(),
00103       bndId_(),
00104       emptyParam_(),
00105       nofElParam_(0),
00106       nofVtxParam_(0),
00107       haveBndParam_( false )
00108     {}
00109 
00111     GridPtr( const GridPtr &org )
00112     : gridPtr_(org.gridPtr_),
00113       elParam_(org.elParam_),
00114       vtxParam_(org.vtxParam_),
00115       bndParam_(org.bndParam_),
00116       bndId_(org.bndId_),
00117       emptyParam_( org.emptyParam_ ),
00118       nofElParam_(org.nofElParam_),
00119       nofVtxParam_(org.nofVtxParam_),
00120       haveBndParam_(org.haveBndParam_)
00121     {}
00122 
00124     GridPtr &operator= ( const GridPtr &org )
00125     {
00126       gridPtr_    = org.gridPtr_;
00127       elParam_    = org.elParam_;
00128       vtxParam_   = org.vtxParam_; 
00129       bndParam_   = org.bndParam_; 
00130       bndId_      = org.bndId_;
00131       emptyParam_ = org.emptyParam_;
00132 
00133       nofElParam_ = org.nofElParam_;
00134       nofVtxParam_ = org.nofVtxParam_; 
00135       haveBndParam_ = org.haveBndParam_; 
00136       return *this;
00137     }
00138     
00140     GridPtr & operator = (GridType * grd)
00141     {
00142       gridPtr_ = std::auto_ptr<GridType>(grd);
00143       elParam_.resize(0);
00144       vtxParam_.resize(0);
00145       bndParam_.resize(0);
00146       bndId_.resize(0);
00147       emptyParam_.resize(0);
00148 
00149       nofVtxParam_ = 0;
00150       nofElParam_ = 0;
00151       haveBndParam_ = false;
00152       return *this;
00153     }
00154 
00156     GridType& operator*() {
00157       return *gridPtr_;
00158     }
00159 
00161     GridType* operator->() {
00162       return gridPtr_.operator -> ();
00163     }
00164 
00166     const GridType& operator*() const {
00167       return *gridPtr_;
00168     }
00169     
00171     const GridType* operator->() const {
00172       return gridPtr_.operator -> ();
00173     }
00174     
00176     GridType* release () {
00177       return gridPtr_.release();
00178     }
00179 
00181     int nofParameters(int cdim) const {
00182       switch (cdim) {
00183       case 0: return nofElParam_; break;
00184       case GridType::dimension: return nofVtxParam_; break;
00185       }
00186       return 0;
00187     }
00188 
00190     template <class Entity>
00191     int nofParameters ( const Entity & ) const
00192     {
00193       return nofParameters( (int) Entity::codimension );
00194     }
00195 
00197     template< class GridImp, template< class > class IntersectionImp >
00198     int nofParameters ( const Intersection< GridImp, IntersectionImp > & intersection ) const
00199     {
00200       return parameters( intersection ).size();
00201     }
00202 
00204     template <class Entity>
00205     const std::vector< double > &parameters ( const Entity &entity ) const
00206     {
00207       typedef typename GridType::LevelGridView GridView;
00208       GridView gridView = gridPtr_->levelView( 0 ); 
00209       switch( (int)Entity::codimension )
00210       {
00211       case 0:
00212         if( nofElParam_ > 0 )
00213         {
00214           assert( (unsigned int)gridView.indexSet().index( entity ) < elParam_.size() );
00215           return elParam_[ gridView.indexSet().index( entity ) ];
00216         }
00217         break;
00218       case GridType::dimension:
00219         if( nofVtxParam_ > 0 )
00220         {
00221           assert( (unsigned int)gridView.indexSet().index( entity ) < vtxParam_.size() );
00222           return vtxParam_[ gridView.indexSet().index( entity ) ];
00223         }
00224         break;
00225       }
00226       return emptyParam_;
00227     }
00228 
00230     template< class GridImp, template< class > class IntersectionImp >
00231     const DGFBoundaryParameter::type & parameters ( const Intersection< GridImp, IntersectionImp > & intersection ) const
00232     {
00233       // if no parameters given return empty vecto
00234       if ( !haveBndParam_ )
00235         return DGFBoundaryParameter::defaultValue();
00236 
00237       return bndParam_[ intersection.boundarySegmentIndex() ];
00238     }
00239 
00240     void loadBalance() 
00241     {
00242       if ( gridPtr_->comm().size() == 1 )
00243         return;
00244       int params = nofElParam_ + nofVtxParam_;
00245       if ( haveBndParam_ )
00246         params += 1;
00247       if ( gridPtr_->comm().max( params ) > 0 ) 
00248       {
00249         DataHandle dh(*this);
00250         gridPtr_->loadBalance( dh.interface() );
00251         gridPtr_->communicate( dh.interface(), InteriorBorder_All_Interface,ForwardCommunication);
00252       } else 
00253       {
00254         gridPtr_->loadBalance();
00255       }
00256     }
00257 
00258   protected:
00259     void initialize ( DGFGridFactory< GridType > &dgfFactory )
00260     {
00261       gridPtr_ = std::auto_ptr< GridType >( dgfFactory.grid() );
00262 
00263       typedef typename GridType::LevelGridView GridView;
00264       GridView gridView = gridPtr_->levelView( 0 ); 
00265       const typename GridView::IndexSet &indexSet = gridView.indexSet();
00266 
00267       nofElParam_ = dgfFactory.template numParameters< 0 >();
00268       nofVtxParam_ = dgfFactory.template numParameters< dimension >();
00269       haveBndParam_ = dgfFactory.haveBoundaryParameters();
00270       
00271       if ( nofElParam_ > 0 )
00272         elParam_.resize( indexSet.size(0) );
00273       if ( nofVtxParam_ > 0 )
00274         vtxParam_.resize( indexSet.size(dimension) );
00275       bndId_.resize( indexSet.size(1) );
00276       if ( haveBndParam_ )
00277         bndParam_.resize( gridPtr_->numBoundarySegments() );
00278 
00279       const PartitionIteratorType partType = Interior_Partition;
00280       typedef typename GridView::template Codim< 0 >::template Partition< partType >::Iterator Iterator;
00281       const Iterator enditer = gridView.template end< 0, partType >();
00282       for( Iterator iter = gridView.template begin< 0, partType >(); iter != enditer; ++iter )
00283       {
00284         const typename Iterator::Entity &el = *iter;
00285         if ( nofElParam_ > 0 ) {
00286           std::swap( elParam_[ indexSet.index(el) ], dgfFactory.parameter(el) );
00287           assert( elParam_[ indexSet.index(el) ].size()  == (size_t)nofElParam_ );
00288         }
00289         if ( nofVtxParam_ > 0 )
00290         {
00291           for ( int v = 0; v < el.template count<dimension>(); ++v) 
00292           {
00293             typename GridView::IndexSet::IndexType index = indexSet.subIndex(el,v,dimension);
00294             if ( vtxParam_[ index ].empty() ) 
00295               std::swap( vtxParam_[ index ], dgfFactory.parameter(*el.template subEntity<dimension>(v) ) );
00296             assert( vtxParam_[ index ].size()  == (size_t)nofVtxParam_ );
00297           }
00298         }
00299         if ( el.hasBoundaryIntersections() )
00300         {
00301           typedef typename GridView::IntersectionIterator IntersectionIterator;
00302           const IntersectionIterator iend = gridView.iend(el);
00303           for( IntersectionIterator iiter = gridView.ibegin(el); iiter != iend; ++iiter )
00304           {
00305             const typename IntersectionIterator::Intersection &inter = *iiter;
00306             // dirty hack: check for "none" to make corner point grid work
00307             if ( inter.boundary() && !inter.type().isNone() )
00308             {
00309               const int k = indexSet.subIndex(el,inter.indexInInside(),1);
00310               bndId_[ k ] = dgfFactory.boundaryId( inter );
00311               if ( haveBndParam_ )
00312                 bndParam_[ inter.boundarySegmentIndex() ] = dgfFactory.boundaryParameter( inter );
00313             }
00314           }
00315         }
00316       }
00317     }
00318 
00319     template <class Entity>
00320     std::vector< double > &params ( const Entity &entity ) 
00321     {
00322       typedef typename GridType::LevelGridView GridView;
00323       GridView gridView = gridPtr_->levelView( 0 ); 
00324       switch( (int)Entity::codimension )
00325       {
00326       case 0:
00327         if( nofElParam_ > 0 ) {
00328           if ( gridView.indexSet().index( entity ) >= elParam_.size() )
00329             elParam_.resize( gridView.indexSet().index( entity ) );
00330           return elParam_[ gridView.indexSet().index( entity ) ];
00331         }
00332         break;
00333       case GridType::dimension:
00334         if( nofVtxParam_ > 0 ) {
00335           if ( gridView.indexSet().index( entity ) >= vtxParam_.size() )
00336             vtxParam_.resize( gridView.indexSet().index( entity ) );
00337           return vtxParam_[ gridView.indexSet().index( entity ) ];
00338         }
00339         break;
00340       }
00341       return emptyParam_;
00342     }
00343 
00344     void setNofParams( int cdim, int nofP )
00345     {
00346       switch (cdim) {
00347       case 0: nofElParam_ = nofP; break;
00348       case GridType::dimension: nofVtxParam_ = nofP; break;
00349       }
00350     }
00351 
00352     struct DataHandle : public CommDataHandleIF<DataHandle,double>
00353     {
00354       DataHandle( GridPtr& gridPtr) : 
00355         gridPtr_(gridPtr),
00356         idSet_(gridPtr->localIdSet())
00357       {
00358         typedef typename GridType::LevelGridView GridView;
00359         GridView gridView = gridPtr_->levelView( 0 ); 
00360         const typename GridView::IndexSet &indexSet = gridView.indexSet();
00361 
00362         const PartitionIteratorType partType = Interior_Partition;
00363         typedef typename GridView::template Codim< 0 >::template Partition< partType >::Iterator Iterator;
00364         const Iterator enditer = gridView.template end< 0, partType >();
00365         for( Iterator iter = gridView.template begin< 0, partType >(); iter != enditer; ++iter )
00366         {
00367           const typename Iterator::Entity &el = *iter;
00368           if ( gridPtr_.nofElParam_ > 0 ) 
00369             std::swap( gridPtr_.elParam_[ indexSet.index(el) ], elData_[ idSet_.id(el) ] );
00370           if ( gridPtr_.nofVtxParam_ > 0 )
00371           {
00372             for ( int v = 0; v < el.template count<dimension>(); ++v) 
00373             {
00374               typename GridView::IndexSet::IndexType index = indexSet.subIndex(el,v,dimension);
00375               if ( ! gridPtr_.vtxParam_[ index ].empty() ) 
00376                 std::swap( gridPtr_.vtxParam_[ index ], vtxData_[ idSet_.subId(el,v,dimension) ] );
00377             }
00378           }
00379         }
00380       }
00381 
00382       ~DataHandle() 
00383       {
00384         typedef typename GridType::LevelGridView GridView;
00385         GridView gridView = gridPtr_->levelView( 0 ); 
00386         const typename GridView::IndexSet &indexSet = gridView.indexSet();
00387 
00388         if ( gridPtr_.nofElParam_ > 0 )
00389           gridPtr_.elParam_.resize( indexSet.size(0) );
00390         if ( gridPtr_.nofVtxParam_ > 0 )
00391           gridPtr_.vtxParam_.resize( indexSet.size(dimension) );
00392 
00393         const PartitionIteratorType partType = All_Partition;
00394         typedef typename GridView::template Codim< 0 >::template Partition< partType >::Iterator Iterator;
00395         const Iterator enditer = gridView.template end< 0, partType >();
00396         for( Iterator iter = gridView.template begin< 0, partType >(); iter != enditer; ++iter )
00397         {
00398           const typename Iterator::Entity &el = *iter;
00399           if ( gridPtr_.nofElParam_ > 0 ) 
00400           {
00401             std::swap( gridPtr_.elParam_[ indexSet.index(el) ], elData_[ idSet_.id(el) ] );
00402             assert( gridPtr_.elParam_[ indexSet.index(el) ].size() == (unsigned int)gridPtr_.nofElParam_ );
00403           }
00404           if ( gridPtr_.nofVtxParam_ > 0 )
00405           {
00406             for ( int v = 0; v < el.template count<dimension>(); ++v) 
00407             {
00408               typename GridView::IndexSet::IndexType index = indexSet.subIndex(el,v,dimension);
00409               if ( gridPtr_.vtxParam_[ index ].empty() ) 
00410                 std::swap( gridPtr_.vtxParam_[ index ], vtxData_[ idSet_.subId(el,v,dimension) ] );
00411               assert( gridPtr_.vtxParam_[ index ].size() == (unsigned int)gridPtr_.nofVtxParam_ );
00412             }
00413           }
00414         }
00415       }
00416 
00417       CommDataHandleIF<DataHandle,double> &interface() 
00418       {
00419         return *this;
00420       }
00421 
00422       bool contains (int dim, int codim) const
00423       {
00424         return (codim==dim || codim==0);
00425       }
00426 
00427       bool fixedsize (int dim, int codim) const
00428       {
00429         return false;
00430       }
00431 
00432       template<class EntityType>
00433       size_t size (const EntityType& e) const
00434       { 
00435         return gridPtr_.nofParameters( (int) e.codimension);
00436       }
00437 
00438       template<class MessageBufferImp, class EntityType>
00439       void gather (MessageBufferImp& buff, const EntityType& e) const
00440       {
00441         const std::vector<double> &v = (e.codimension==0)? elData_[idSet_.id(e)]:vtxData_[idSet_.id(e)];
00442         const size_t s = v.size();
00443         for (size_t i=0;i<s;++i)
00444           buff.write( v[i] );
00445         assert( s == (size_t)gridPtr_.nofParameters(e.codimension) );
00446       }
00447 
00448       template<class MessageBufferImp, class EntityType>
00449       void scatter (MessageBufferImp& buff, const EntityType& e, size_t n)
00450       { 
00451         std::vector<double> &v = (e.codimension==0)? elData_[idSet_.id(e)]:vtxData_[idSet_.id(e)];
00452         v.resize( n );
00453         gridPtr_.setNofParams( e.codimension, n );
00454         for (size_t i=0;i<n;++i)
00455           buff.read( v[i] );
00456       }
00457 
00458       private:
00459       typedef typename GridType::LocalIdSet IdSet;
00460       GridPtr &gridPtr_;
00461       const IdSet &idSet_;
00462       mutable std::map< typename IdSet::IdType, std::vector<double> > elData_, vtxData_;
00463     };
00464 
00465     // grid auto pointer
00466     mutable std::auto_ptr<GridType> gridPtr_;
00467     // element and vertex parameters
00468     std::vector< std::vector< double > > elParam_;
00469     std::vector< std::vector< double > > vtxParam_;
00470     std::vector< DGFBoundaryParameter::type > bndParam_;
00471     std::vector< int > bndId_;
00472     std::vector < double > emptyParam_;
00473 
00474     int  nofElParam_;
00475     int  nofVtxParam_;
00476     bool haveBndParam_;
00477   }; // end of class GridPtr
00478 
00479 } // end namespace Dune
00480 
00481 #endif