dune-grid  2.2.0
albertagrid/gridfactory.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_ALBERTA_GRIDFACTORY_HH
00005 #define DUNE_ALBERTA_GRIDFACTORY_HH
00006 
00012 #include <algorithm>
00013 #include <limits>
00014 #include <map>
00015 
00016 #include <dune/common/array.hh>
00017 
00018 #include <dune/geometry/referenceelements.hh>
00019 
00020 #include <dune/grid/common/gridfactory.hh>
00021 
00022 #include <dune/grid/utility/grapedataioformattypes.hh>
00023 
00024 #include <dune/grid/albertagrid/agrid.hh>
00025 
00026 #if HAVE_ALBERTA
00027 
00028 namespace Dune
00029 {
00030 
00048   template< int dim, int dimworld >
00049   class GridFactory< AlbertaGrid< dim, dimworld > >
00050   : public GridFactoryInterface< AlbertaGrid< dim, dimworld > >
00051   {
00052     typedef GridFactory< AlbertaGrid< dim, dimworld > > This;
00053 
00054   public:
00056     typedef AlbertaGrid< dim, dimworld > Grid;
00057 
00059     typedef typename Grid::ctype ctype;
00060 
00062     static const int dimension = Grid::dimension;
00064     static const int dimensionworld = Grid::dimensionworld;
00065 
00067     typedef FieldVector< ctype, dimensionworld > WorldVector;
00069     typedef FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix;
00070 
00071     typedef DuneBoundaryProjection< dimensionworld > DuneProjection;
00072     typedef Dune::shared_ptr< const DuneProjection > DuneProjectionPtr;
00073     typedef Dune::BoundarySegment< dimension, dimensionworld > BoundarySegment;
00074 
00075     template< int codim >
00076     struct Codim
00077     {
00078       typedef typename Grid::template Codim< codim >::Entity Entity;
00079     };
00080 
00081   private:
00082     typedef Dune::BoundarySegmentWrapper< dimension, dimensionworld > BoundarySegmentWrapper;
00083 
00084     static const int numVertices
00085       = Alberta::NumSubEntities< dimension, dimension >::value;
00086 
00087     typedef Alberta::MacroElement< dimension > MacroElement;
00088     typedef Alberta::ElementInfo< dimension > ElementInfo;
00089     typedef Alberta::MacroData< dimension > MacroData;
00090     typedef Alberta::NumberingMap< dimension, Alberta::Dune2AlbertaNumbering > NumberingMap;
00091     typedef Alberta::DuneBoundaryProjection< dimension > Projection;
00092 
00093     typedef array< unsigned int, dimension > FaceId;
00094     typedef std::map< FaceId, size_t > BoundaryMap;
00095 
00096     class ProjectionFactory;
00097 
00098   public:
00100     static const bool supportsBoundaryIds = true;
00102     static const bool supportPeriodicity = MacroData::supportPeriodicity;
00103 
00105     GridFactory ()
00106     : globalProjection_( (const DuneProjection *)0 )
00107     {
00108       macroData_.create();
00109     }
00110 
00111     virtual ~GridFactory ();
00112 
00117     virtual void insertVertex ( const WorldVector &pos )
00118     {
00119       macroData_.insertVertex( pos );
00120     }
00121 
00127     virtual void insertElement ( const GeometryType &type,
00128                                  const std::vector< unsigned int > &vertices )
00129     {
00130       if( (int)type.dim() != dimension )
00131         DUNE_THROW( AlbertaError, "Inserting element of wrong dimension: " << type.dim() );
00132       if( !type.isSimplex() )
00133         DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
00134 
00135       if( vertices.size() != (size_t)numVertices )
00136         DUNE_THROW( AlbertaError, "Wrong number of vertices passed: " << vertices.size() << "." );
00137 
00138       int array[ numVertices ];
00139       for( int i = 0; i < numVertices; ++i )
00140         array[ i ] = vertices[ numberingMap_.alberta2dune( dimension, i ) ];
00141       macroData_.insertElement( array );
00142     }
00143 
00152     virtual void insertBoundary ( int element, int face, int id )
00153     {
00154       if( (id <= 0) || (id > 127) )
00155         DUNE_THROW( AlbertaError, "Invalid boundary id: " << id << "." );
00156       macroData_.boundaryId( element, numberingMap_.dune2alberta( 1, face ) ) = id;
00157     }
00158 
00167     virtual void
00168     insertBoundaryProjection ( const GeometryType &type,
00169                                const std::vector< unsigned int > &vertices,
00170                                const DuneProjection *projection )
00171     {
00172       if( (int)type.dim() != dimension-1 )
00173         DUNE_THROW( AlbertaError, "Inserting boundary face of wrong dimension: " << type.dim() );
00174       if( !type.isSimplex() )
00175         DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
00176 
00177       FaceId faceId;
00178       if( vertices.size() != faceId.size() )
00179         DUNE_THROW( AlbertaError, "Wrong number of face vertices passed: " << vertices.size() << "." );
00180       for( size_t i = 0; i < faceId.size(); ++i )
00181         faceId[ i ] = vertices[ i ];
00182       std::sort( faceId.begin(), faceId.end() );
00183 
00184       typedef std::pair< typename BoundaryMap::iterator, bool > InsertResult;
00185       const InsertResult result = boundaryMap_.insert( std::make_pair( faceId, boundaryProjections_.size() ) );
00186       if( !result.second )
00187         DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
00188       boundaryProjections_.push_back( DuneProjectionPtr( projection ) );
00189     }
00190 
00191 
00198     virtual void insertBoundaryProjection ( const DuneProjection *projection )
00199     {
00200       if( globalProjection_ )
00201         DUNE_THROW( GridError, "Only one global boundary projection can be attached to a grid." );
00202       globalProjection_ = DuneProjectionPtr( projection );
00203     }
00204 
00210     virtual void
00211     insertBoundarySegment ( const std::vector< unsigned int >& vertices )
00212     {
00213       typedef typename GenericGeometry::SimplexTopology< dimension-1 >::type Topology;
00214       insertBoundaryProjection( GeometryType( Topology() ), vertices, 0 );
00215     }
00216 
00222     virtual void
00223     insertBoundarySegment ( const std::vector< unsigned int > &vertices,
00224                             const shared_ptr< BoundarySegment > &boundarySegment )
00225     {
00226       const GenericReferenceElement< ctype, dimension-1 > &refSimplex
00227         = GenericReferenceElements< ctype, dimension-1 >::simplex();
00228 
00229       if( !boundarySegment )
00230         DUNE_THROW( GridError, "Trying to insert null as a boundary segment." );
00231       if( (int)vertices.size() != refSimplex.size( dimension-1 ) )
00232         DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
00233 
00234       std::vector< WorldVector > coords( refSimplex.size( dimension-1 ) );
00235       for( int i = 0; i < dimension; ++i )
00236       {
00237         Alberta::GlobalVector &x = macroData_.vertex( vertices[ i ] );
00238         for( int j = 0; j < dimensionworld; ++j )
00239           coords[ i ][ j ] = x[ j ];
00240         if( ((*boundarySegment)( refSimplex.position( i, dimension-1 ) ) - coords[ i ]).two_norm() > 1e-6 )
00241           DUNE_THROW( GridError, "Boundary segment does not interpolate the corners." );
00242       }
00243 
00244       const GeometryType gt = refSimplex.type( 0, 0 );
00245       const DuneProjection *prj = new BoundarySegmentWrapper( gt, coords, boundarySegment );
00246       insertBoundaryProjection( gt, vertices, prj );
00247     }
00248 
00262     void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
00263 
00272     void markLongestEdge ()
00273     {
00274       macroData_.markLongestEdge();
00275     }
00276 
00289     Grid *createGrid ()
00290     {
00291       macroData_.finalize();
00292       if( macroData_.elementCount() == 0 )
00293         DUNE_THROW( GridError, "Cannot create empty AlbertaGrid." );
00294       if( dimension < 3 )
00295         macroData_.setOrientation( Alberta::Real( 1 ) );
00296       assert( macroData_.checkNeighbors() );
00297       macroData_.checkCycles();
00298       ProjectionFactory projectionFactory( *this );
00299       return new Grid( macroData_, projectionFactory );
00300     }
00301 
00306     static void destroyGrid ( Grid *grid )
00307     {
00308       delete grid;
00309     }
00310 
00319     template< GrapeIOFileFormatType type >
00320     bool write ( const std::string &filename )
00321     {
00322       dune_static_assert( type != pgm, "AlbertaGridFactory: writing pgm format is not supported." );
00323       macroData_.finalize();
00324       if( dimension < 3 )
00325         macroData_.setOrientation( Alberta::Real( 1 ) );
00326       assert( macroData_.checkNeighbors() );
00327       return macroData_.write( filename, (type == xdr) );
00328     }
00329 
00338     virtual bool write ( const std::string &filename )
00339     {
00340       return write< ascii >( filename );
00341     }
00342 
00343     virtual unsigned int
00344     insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
00345     {
00346       return insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
00347     }
00348 
00349     virtual unsigned int
00350     insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
00351     {
00352       const int elIndex = insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
00353       const typename MacroData::ElementId &elementId = macroData_.element( elIndex );
00354       return elementId[ Grid::getRealImplementation( entity ).subEntity() ];
00355     }
00356 
00357     virtual unsigned int
00358     insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
00359     {
00360       const Grid &grid = Grid::getRealImplementation( intersection ).grid();
00361       const ElementInfo &elementInfo = Grid::getRealImplementation( intersection ).elementInfo();
00362       const int face = grid.generic2alberta( 1, intersection.indexInInside() );
00363       return insertionIndex( elementInfo, face );
00364     }
00365 
00366     virtual bool
00367     wasInserted ( const typename Grid::LeafIntersection &intersection ) const
00368     {
00369       return (insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
00370     }
00371 
00372   private:
00373     unsigned int insertionIndex ( const ElementInfo &elementInfo ) const;
00374     unsigned int insertionIndex ( const ElementInfo &elementInfo, const int face ) const;
00375 
00376     FaceId faceId ( const ElementInfo &elementInfo, const int face ) const;
00377 
00378     MacroData macroData_;
00379     NumberingMap numberingMap_;
00380     DuneProjectionPtr globalProjection_;
00381     BoundaryMap boundaryMap_;
00382     std::vector< DuneProjectionPtr > boundaryProjections_;
00383   };
00384 
00385 
00386   template< int dim, int dimworld >
00387   GridFactory< AlbertaGrid< dim, dimworld > >::~GridFactory ()
00388   {
00389     macroData_.release();
00390   }
00391 
00392 
00393   template< int dim, int dimworld >
00394   inline void
00395   GridFactory< AlbertaGrid< dim, dimworld > >
00396     ::insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift )
00397   {
00398     // make sure the matrix is orthogonal
00399     for( int i = 0; i < dimworld; ++i )
00400       for( int j = 0; j < dimworld; ++j )
00401       {
00402         const ctype delta = (i == j ? ctype( 1 ) : ctype( 0 ));
00403         const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
00404 
00405         if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
00406         {
00407           DUNE_THROW( AlbertaError,
00408                       "Matrix of face transformation is not orthogonal." );
00409         }
00410       }
00411 
00412     // copy matrix
00413     Alberta::GlobalMatrix M;
00414     for( int i = 0; i < dimworld; ++i )
00415       for( int j = 0; j < dimworld; ++j )
00416         M[ i ][ j ] = matrix[ i ][ j ];
00417 
00418     // copy shift
00419     Alberta::GlobalVector t;
00420     for( int i = 0; i < dimworld; ++i )
00421       t[ i ] = shift[ i ];
00422 
00423     // insert into ALBERTA macro data
00424     macroData_.insertWallTrafo( M, t );
00425   }
00426 
00427 
00428   template< int dim, int dimworld >
00429   inline unsigned int
00430   GridFactory< AlbertaGrid< dim, dimworld > >
00431     ::insertionIndex ( const ElementInfo &elementInfo ) const
00432   {
00433     const MacroElement &macroElement = elementInfo.macroElement();
00434     const unsigned int index = macroElement.index;
00435 
00436 #ifndef NDEBUG
00437     const typename MacroData::ElementId &elementId = macroData_.element( index );
00438     for( int i = 0; i <= dimension; ++i )
00439     {
00440       const Alberta::GlobalVector &x = macroData_.vertex( elementId[ i ] );
00441       const Alberta::GlobalVector &y = macroElement.coordinate( i );
00442       for( int j = 0; j < dimensionworld; ++j )
00443       {
00444         if( x[ j ] != y[ j ] )
00445           DUNE_THROW( GridError, "Vertex in macro element does not coincide with same vertex in macro data structure." );
00446       }
00447     }
00448 #endif // #ifndef NDEBUG
00449 
00450     return index;
00451   }
00452 
00453 
00454   template< int dim, int dimworld >
00455   inline unsigned int
00456   GridFactory< AlbertaGrid< dim, dimworld > >
00457     ::insertionIndex ( const ElementInfo &elementInfo, const int face ) const
00458   {
00459     typedef typename BoundaryMap::const_iterator Iterator;
00460     const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
00461     if( it != boundaryMap_.end() )
00462       return it->second;
00463     else
00464       return std::numeric_limits< unsigned int >::max();
00465   }
00466 
00467 
00468   template< int dim, int dimworld >
00469   inline typename GridFactory< AlbertaGrid< dim, dimworld > >::FaceId
00470   GridFactory< AlbertaGrid< dim, dimworld > >
00471     ::faceId ( const ElementInfo &elementInfo, const int face ) const
00472   {
00473     const unsigned int index = insertionIndex( elementInfo );
00474     const typename MacroData::ElementId &elementId = macroData_.element( index );
00475 
00476     FaceId faceId;
00477     for( size_t i = 0; i < faceId.size(); ++i )
00478     {
00479       const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
00480       faceId[ i ] = elementId[ k ];
00481     }
00482     std::sort( faceId.begin(), faceId.end() );
00483     return faceId;
00484   }
00485 
00486 
00487 
00488   // GridFactory::ProjectionFactory
00489   // ------------------------------
00490 
00491   template< int dim, int dimworld >
00492   class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
00493   : public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
00494   {
00495     typedef ProjectionFactory This;
00496     typedef Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory > Base;
00497 
00498     typedef typename Dune::GridFactory< AlbertaGrid< dim, dimworld > > Factory;
00499 
00500   public:
00501     typedef typename Base::Projection Projection;
00502     typedef typename Base::ElementInfo ElementInfo;
00503 
00504     typedef typename Projection::Projection DuneProjection;
00505 
00506     ProjectionFactory( const GridFactory &gridFactory )
00507     : gridFactory_( gridFactory )
00508     {}
00509 
00510     bool hasProjection ( const ElementInfo &elementInfo, const int face ) const
00511     {
00512       if( gridFactory().globalProjection_ )
00513         return true;
00514 
00515       const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
00516       if( index < std::numeric_limits< unsigned int >::max() )
00517         return bool( gridFactory().boundaryProjections_[ index ] );
00518       else
00519         return false;
00520     }
00521 
00522     bool hasProjection ( const ElementInfo &elementInfo ) const
00523     {
00524       return bool( gridFactory().globalProjection_ );
00525     }
00526 
00527     Projection projection ( const ElementInfo &elementInfo, const int face ) const
00528     {
00529       const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
00530       if( index < std::numeric_limits< unsigned int >::max() )
00531       {
00532         const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
00533         if( projection )
00534           return Projection( projection );
00535       }
00536 
00537       assert( gridFactory().globalProjection_ );
00538       return Projection( gridFactory().globalProjection_ );
00539     };
00540 
00541     Projection projection ( const ElementInfo &elementInfo ) const
00542     {
00543       assert( gridFactory().globalProjection_ );
00544       return Projection( gridFactory().globalProjection_ );
00545     };
00546 
00547     const GridFactory &gridFactory () const
00548     {
00549       return gridFactory_;
00550     }
00551 
00552   private:
00553     const GridFactory &gridFactory_;
00554   };
00555 
00556 }
00557 
00558 #endif // #if HAVE_ALBERTA
00559 
00560 #endif // #ifndef DUNE_ALBERTA_GRIDFACTORY_HH