dune-grid
2.2.0
|
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 > ¶meters ( 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 > ¶ms ( 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