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_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 ¯oElement = 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