dune-grid  2.2.0
alu3dgridfactory.hh
Go to the documentation of this file.
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