dune-grid
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=8 sw=2 sts=2: 00003 00004 #ifndef DUNE_ALU3DGRID_FACTORY_HH 00005 #define DUNE_ALU3DGRID_FACTORY_HH 00006 00007 #include <map> 00008 #include <vector> 00009 00010 #ifdef ENABLE_ALUGRID 00011 00012 #include <dune/common/array.hh> 00013 #include <dune/common/mpihelper.hh> 00014 00015 #include <dune/geometry/referenceelements.hh> 00016 00017 #include <dune/grid/common/gridfactory.hh> 00018 #include <dune/grid/common/boundaryprojection.hh> 00019 00020 #include <dune/grid/alugrid/common/transformation.hh> 00021 #include <dune/grid/alugrid/3d/alugrid.hh> 00022 00023 namespace Dune 00024 { 00025 00027 template< class ALUGrid > 00028 class ALU3dGridFactory 00029 : public GridFactoryInterface< ALUGrid > 00030 { 00031 typedef ALU3dGridFactory< ALUGrid > ThisType; 00032 typedef GridFactoryInterface< ALUGrid > BaseType; 00033 00034 public: 00035 typedef ALUGrid Grid; 00036 00037 typedef typename Grid::ctype ctype; 00038 00039 static const ALU3dGridElementType elementType = Grid::elementType; 00040 00041 static const unsigned int dimension = Grid::dimension; 00042 static const unsigned int dimensionworld = Grid::dimensionworld; 00043 00044 typedef typename Grid::MPICommunicatorType MPICommunicatorType; 00045 00047 typedef DuneBoundaryProjection< 3 > DuneBoundaryProjectionType; 00048 00049 template< int codim > 00050 struct Codim 00051 { 00052 typedef typename Grid::template Codim< codim >::Entity Entity; 00053 }; 00054 00055 typedef unsigned int VertexId; 00056 00057 typedef ALUGridTransformation< ctype, dimensionworld > Transformation; 00058 00060 typedef typename Transformation::WorldVector WorldVector; 00062 typedef typename Transformation::WorldMatrix WorldMatrix; 00063 00064 private: 00065 dune_static_assert( (elementType == tetra || elementType == hexa), 00066 "ALU3dGridFactory supports only grids containing " 00067 "tetrahedrons or hexahedrons exclusively." ); 00068 00069 typedef Dune::BoundarySegmentWrapper< dimension, dimensionworld > BoundarySegmentWrapperType; 00070 00071 static const unsigned int numCorners = EntityCount< elementType >::numVertices; 00072 static const unsigned int numFaces = EntityCount< elementType >::numFaces; 00073 static const unsigned int numFaceCorners = EntityCount< elementType >::numVerticesPerFace; 00074 00075 typedef ElementTopologyMapping< elementType > ElementTopologyMappingType; 00076 typedef FaceTopologyMapping< elementType > FaceTopologyMappingType; 00077 00078 typedef FieldVector< ctype, dimensionworld > VertexType; 00079 typedef std::vector< unsigned int > ElementType; 00080 typedef array< unsigned int, numFaceCorners > FaceType; 00081 00082 struct FaceLess; 00083 00084 typedef std::vector< std::pair< VertexType, size_t > > VertexVector; 00085 typedef std::vector< ElementType > ElementVector; 00086 typedef std::pair< FaceType, int > BndPair ; 00087 typedef std::map< FaceType, int > BoundaryIdMap; 00088 typedef std::vector< std::pair< BndPair, BndPair > > PeriodicBoundaryVector; 00089 typedef std::pair< unsigned int, int > SubEntity; 00090 typedef std::map< FaceType, SubEntity, FaceLess > FaceMap; 00091 00092 typedef std::map< FaceType, const DuneBoundaryProjectionType* > BoundaryProjectionMap; 00093 typedef std::vector< const DuneBoundaryProjectionType* > BoundaryProjectionVector; 00094 00095 typedef std::vector< Transformation > FaceTransformationVector; 00096 00097 // copy vertex numbers and store smalled #dimension ones 00098 void copyAndSort ( const std::vector< unsigned int > &vertices, FaceType &faceId ) const 00099 { 00100 std::vector<unsigned int> tmp( vertices ); 00101 std::sort( tmp.begin(), tmp.end() ); 00102 00103 // copy only the first dimension vertices (enough for key) 00104 for( size_t i = 0; i < faceId.size(); ++i ) faceId[ i ] = tmp[ i ]; 00105 } 00106 00107 private: 00108 // return grid object 00109 virtual Grid* createGridObj( BoundaryProjectionVector* bndProjections, const std::string& name ) const 00110 { 00111 return ( allowGridGeneration_ ) ? 00112 new Grid( communicator_, globalProjection_, bndProjections , name, realGrid_ ) : 00113 new Grid( communicator_ ); 00114 } 00115 00116 public: 00118 explicit ALU3dGridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator(), 00119 bool removeGeneratedFile = true ); 00120 00122 explicit ALU3dGridFactory ( const std::string &filename, 00123 const MPICommunicatorType &communicator = Grid::defaultCommunicator() ); 00124 00126 explicit ALU3dGridFactory ( const bool verbose, const MPICommunicatorType &communicator ); 00127 00129 virtual ~ALU3dGridFactory (); 00130 00135 virtual void insertVertex ( const VertexType &pos ); 00136 00137 // for testing parallel GridFactory 00138 VertexId insertVertex ( const VertexType &pos, const size_t globalId ); 00139 00148 virtual void 00149 insertElement ( const GeometryType &geometry, 00150 const std::vector< VertexId > &vertices ); 00151 00162 virtual void 00163 insertBoundary ( const GeometryType &geometry, 00164 const std::vector< VertexId > &faceVertices, 00165 const int id ); 00166 00173 virtual void insertBoundary ( const int element, const int face, const int id ); 00174 00175 // for testing parallel GridFactory 00176 void insertProcessBorder ( const int element, const int face ) 00177 { 00178 insertBoundary( element, face, ALU3DSPACE ProcessorBoundary_t ); 00179 } 00180 00189 virtual void 00190 insertBoundaryProjection ( const GeometryType &type, 00191 const std::vector< VertexId > &vertices, 00192 const DuneBoundaryProjectionType *projection ); 00193 00198 virtual void 00199 insertBoundarySegment ( const std::vector< VertexId >& vertices ) ; 00200 00206 virtual void 00207 insertBoundarySegment ( const std::vector< VertexId >& vertices, 00208 const shared_ptr<BoundarySegment<3,3> >& boundarySegment ) ; 00209 00214 virtual void insertBoundaryProjection ( const DuneBoundaryProjectionType& bndProjection ); 00215 00225 void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift ); 00226 00231 Grid *createGrid (); 00232 00233 Grid *createGrid ( const bool addMissingBoundaries, const std::string dgfName = "" ); 00234 00235 Grid *createGrid ( const bool addMissingBoundaries, bool temporary, const std::string dgfName = "" ); 00236 00237 virtual unsigned int 00238 insertionIndex ( const typename Codim< 0 >::Entity &entity ) const 00239 { 00240 return Grid::getRealImplementation( entity ).getIndex(); 00241 } 00242 virtual unsigned int 00243 insertionIndex ( const typename Codim< dimension >::Entity &entity ) const 00244 { 00245 return Grid::getRealImplementation( entity ).getIndex(); 00246 } 00247 virtual unsigned int 00248 insertionIndex ( const typename Grid::LeafIntersection &intersection ) const 00249 { 00250 return intersection.boundarySegmentIndex(); 00251 } 00252 virtual bool 00253 wasInserted ( const typename Grid::LeafIntersection &intersection ) const 00254 { 00255 return intersection.boundary() && 00256 ( insertionIndex(intersection) < numFacesInserted_ ); 00257 } 00258 00259 private: 00260 size_t globalId ( const VertexId &id ) const 00261 { 00262 assert( id < vertices_.size() ); 00263 return vertices_[ id ].second; 00264 } 00265 00266 const VertexType &position ( const VertexId &id ) const 00267 { 00268 assert( id < vertices_.size() ); 00269 return vertices_[ id ].first; 00270 } 00271 00272 void assertGeometryType( const GeometryType &geometry ); 00273 static void generateFace ( const ElementType &element, const int f, FaceType &face ); 00274 void generateFace ( const SubEntity &subEntity, FaceType &face ) const; 00275 void correctElementOrientation (); 00276 bool identifyFaces ( const Transformation &transformation, const FaceType &key1, const FaceType &key2, const int defaultId ); 00277 void searchPeriodicNeighbor ( FaceMap &faceMap, const typename FaceMap::iterator &pos, const int defaultId ); 00278 void reinsertBoundary ( const FaceMap &faceMap, const typename FaceMap::const_iterator &pos, const int id ); 00279 void recreateBoundaryIds ( const int defaultId = 1 ); 00280 00281 int rank_; 00282 00283 VertexVector vertices_; 00284 ElementVector elements_; 00285 BoundaryIdMap boundaryIds_; 00286 PeriodicBoundaryVector periodicBoundaries_; 00287 const DuneBoundaryProjectionType* globalProjection_ ; 00288 BoundaryProjectionMap boundaryProjections_; 00289 FaceTransformationVector faceTransformations_; 00290 unsigned int numFacesInserted_; 00291 bool realGrid_; 00292 const bool allowGridGeneration_; 00293 00294 MPICommunicatorType communicator_; 00295 }; 00296 00297 00298 00299 template< class ALUGrid > 00300 struct ALU3dGridFactory< ALUGrid >::FaceLess 00301 : public std::binary_function< FaceType, FaceType, bool > 00302 { 00303 bool operator() ( const FaceType &a, const FaceType &b ) const 00304 { 00305 for( unsigned int i = 0; i < numFaceCorners; ++i ) 00306 { 00307 if( a[ i ] != b[ i ] ) 00308 return (a[ i ] < b[ i ]); 00309 } 00310 return false; 00311 } 00312 }; 00313 00314 00315 template< class ALUGrid > 00316 inline void ALU3dGridFactory< ALUGrid > 00317 ::assertGeometryType( const GeometryType &geometry ) 00318 { 00319 if( elementType == tetra ) 00320 { 00321 if( !geometry.isSimplex() ) 00322 DUNE_THROW( GridError, "Only simplex geometries can be inserted into " 00323 "ALUSimplexGrid< 3, 3 >." ); 00324 } 00325 else 00326 { 00327 if( !geometry.isCube() ) 00328 DUNE_THROW( GridError, "Only cube geometries can be inserted into " 00329 "ALUCubeGrid< 3, 3 >." ); 00330 } 00331 } 00332 00333 00337 template<> 00338 class GridFactory< ALUSimplexGrid< 3, 3 > > 00339 : public ALU3dGridFactory< ALUSimplexGrid< 3, 3 > > 00340 { 00341 typedef GridFactory< ALUSimplexGrid< 3, 3 > > ThisType; 00342 typedef ALU3dGridFactory< ALUSimplexGrid< 3, 3 > > BaseType; 00343 00344 public: 00345 typedef BaseType::Grid Grid; 00346 00347 typedef BaseType::MPICommunicatorType MPICommunicatorType; 00348 00350 explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00351 : BaseType( communicator ) 00352 {} 00353 00355 GridFactory ( const std::string &filename, 00356 const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00357 : BaseType( filename, communicator ) 00358 {} 00359 00360 protected: 00361 template< class, class, int > friend class ALULocalGeometryStorage; 00363 GridFactory ( const bool realGrid, 00364 const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00365 : BaseType( realGrid, communicator ) 00366 {} 00367 }; 00368 00369 00370 00374 template<> 00375 class GridFactory< ALUCubeGrid< 3, 3 > > 00376 : public ALU3dGridFactory< ALUCubeGrid< 3, 3 > > 00377 { 00378 typedef GridFactory< ALUCubeGrid< 3, 3 > > ThisType; 00379 typedef ALU3dGridFactory< ALUCubeGrid< 3, 3 > > BaseType; 00380 00381 public: 00382 typedef BaseType::Grid Grid; 00383 00384 typedef BaseType::MPICommunicatorType MPICommunicatorType; 00385 00387 explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00388 : BaseType( communicator ) 00389 {} 00390 00392 GridFactory ( const std::string &filename, 00393 const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00394 : BaseType( filename, communicator ) 00395 {} 00396 00397 protected: 00398 template< class, class, int > friend class ALULocalGeometryStorage; 00400 GridFactory ( const bool realGrid, 00401 const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00402 : BaseType( realGrid, communicator ) 00403 {} 00404 }; 00405 00406 00410 template<ALUGridElementType eltype, ALUGridRefinementType refinementtype , class Comm > 00411 class GridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > 00412 : public ALU3dGridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > 00413 { 00414 typedef GridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > ThisType; 00415 typedef ALU3dGridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > BaseType; 00416 00417 public: 00418 typedef typename BaseType::Grid Grid; 00419 00420 typedef typename BaseType::MPICommunicatorType MPICommunicatorType; 00421 00423 explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00424 : BaseType( communicator ) 00425 {} 00426 00428 GridFactory ( const std::string &filename, 00429 const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00430 : BaseType( filename, communicator ) 00431 {} 00432 00433 protected: 00434 template< class, class, int > friend class ALULocalGeometryStorage; 00436 GridFactory ( const bool realGrid, 00437 const MPICommunicatorType &communicator = Grid::defaultCommunicator() ) 00438 : BaseType( realGrid, communicator ) 00439 {} 00440 }; 00441 00442 00443 00444 // Implementation of ALU3dGridFactory 00445 // ---------------------------------- 00446 00447 template< class ALUGrid > 00448 inline 00449 ALU3dGridFactory< ALUGrid > 00450 :: ALU3dGridFactory ( const MPICommunicatorType &communicator, 00451 bool removeGeneratedFile ) 00452 : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ), 00453 globalProjection_ ( 0 ), 00454 numFacesInserted_ ( 0 ), 00455 realGrid_( true ), 00456 allowGridGeneration_( rank_ == 0 ), 00457 communicator_( communicator ) 00458 {} 00459 00460 template< class ALUGrid > 00461 inline 00462 ALU3dGridFactory< ALUGrid > 00463 :: ALU3dGridFactory ( const std::string &filename, 00464 const MPICommunicatorType &communicator ) 00465 : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ), 00466 globalProjection_ ( 0 ), 00467 numFacesInserted_ ( 0 ), 00468 realGrid_( true ), 00469 allowGridGeneration_( rank_ == 0 ), 00470 communicator_( communicator ) 00471 {} 00472 00473 template< class ALUGrid > 00474 inline 00475 ALU3dGridFactory< ALUGrid > 00476 :: ALU3dGridFactory ( const bool realGrid, 00477 const MPICommunicatorType &communicator ) 00478 : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ), 00479 globalProjection_ ( 0 ), 00480 numFacesInserted_ ( 0 ), 00481 realGrid_( realGrid ), 00482 allowGridGeneration_( true ), 00483 communicator_( communicator ) 00484 {} 00485 00486 template< class ALUGrid > 00487 inline void ALU3dGridFactory< ALUGrid > :: 00488 insertBoundarySegment ( const std::vector< unsigned int >& vertices ) 00489 { 00490 FaceType faceId; 00491 copyAndSort( vertices, faceId ); 00492 00493 if( vertices.size() != numFaceCorners ) 00494 DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." ); 00495 00496 if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() ) 00497 DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." ); 00498 // DUNE_THROW( NotImplemented, "insertBoundarySegment with a single argument" ); 00499 00500 boundaryProjections_[ faceId ] = 0; 00501 00502 BndPair boundaryId; 00503 for( unsigned int i = 0; i < numFaceCorners; ++i ) 00504 { 00505 const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i ); 00506 boundaryId.first[ j ] = vertices[ i ]; 00507 } 00508 boundaryId.second = 1; 00509 boundaryIds_.insert( boundaryId ); 00510 } 00511 00512 template< class ALUGrid > 00513 inline void ALU3dGridFactory< ALUGrid > :: 00514 insertBoundarySegment ( const std::vector< unsigned int >& vertices, 00515 const shared_ptr<BoundarySegment<3,3> >& boundarySegment ) 00516 { 00517 FaceType faceId; 00518 copyAndSort( vertices, faceId ); 00519 00520 if( vertices.size() != numFaceCorners ) 00521 DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." ); 00522 00523 if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() ) 00524 DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." ); 00525 00526 const size_t numVx = vertices.size(); 00527 GeometryType type; 00528 if( numVx == 3 ) 00529 type.makeSimplex( dimension-1 ); 00530 else 00531 type.makeCube( dimension-1 ); 00532 00533 // we need double here because of the structure of BoundarySegment 00534 // and BoundarySegmentWrapper which have double as coordinate type 00535 typedef FieldVector< double, dimensionworld > CoordType; 00536 std::vector< CoordType > coords( numVx ); 00537 for( size_t i = 0; i < numVx; ++i ) 00538 { 00539 // if this assertions is thrown vertices were not inserted at first 00540 assert( vertices_.size() > vertices[ i ] ); 00541 00542 // get global coordinate and copy it 00543 const VertexType &x = position( vertices[ i ] ); 00544 for( unsigned int j = 0; j < dimensionworld; ++j ) 00545 coords[ i ][ j ] = x[ j ]; 00546 } 00547 00548 BoundarySegmentWrapperType* prj 00549 = new BoundarySegmentWrapperType( type, coords, boundarySegment ); 00550 boundaryProjections_[ faceId ] = prj; 00551 #ifndef NDEBUG 00552 // consistency check 00553 for( size_t i = 0; i < numVx; ++i ) 00554 { 00555 CoordType global = (*prj)( coords [ i ] ); 00556 if( (global - coords[ i ]).two_norm() > 1e-6 ) 00557 DUNE_THROW(GridError,"BoundarySegment does not map face vertices to face vertices."); 00558 } 00559 #endif 00560 00561 BndPair boundaryId; 00562 for( unsigned int i = 0; i < numFaceCorners; ++i ) 00563 { 00564 const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i ); 00565 boundaryId.first[ j ] = vertices[ i ]; 00566 } 00567 boundaryId.second = 1; 00568 boundaryIds_.insert( boundaryId ); 00569 } 00570 00571 00572 template< class ALUGrid > 00573 inline void ALU3dGridFactory< ALUGrid > 00574 ::generateFace ( const SubEntity &subEntity, FaceType &face ) const 00575 { 00576 generateFace( elements_[ subEntity.first ], subEntity.second, face ); 00577 } 00578 00579 } // end namespace Dune 00580 00581 #endif // #ifdef ENABLE_ALUGRID 00582 00583 #if COMPILE_ALUGRID_INLINE 00584 #include "alu3dgridfactory.cc" 00585 #endif 00586 #endif