4 #ifndef DUNE_ALBERTA_GRIDFACTORY_HH
5 #define DUNE_ALBERTA_GRIDFACTORY_HH
18 #include <dune/common/to_unique_ptr.hh>
20 #include <dune/geometry/referenceelements.hh>
48 template<
int dim,
int dimworld >
49 class GridFactory< AlbertaGrid< dim, dimworld > >
50 :
public GridFactoryInterface< AlbertaGrid< dim, dimworld > >
52 typedef GridFactory< AlbertaGrid< dim, dimworld > > This;
56 typedef AlbertaGrid< dim, dimworld > Grid;
67 typedef FieldVector< ctype, dimensionworld > WorldVector;
69 typedef FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix;
71 typedef DuneBoundaryProjection< dimensionworld > DuneProjection;
72 typedef std::shared_ptr< const DuneProjection > DuneProjectionPtr;
84 static const int numVertices
85 = Alberta::NumSubEntities< dimension, dimension >::value;
87 typedef Alberta::MacroElement< dimension > MacroElement;
88 typedef Alberta::ElementInfo< dimension > ElementInfo;
89 typedef Alberta::MacroData< dimension > MacroData;
90 typedef Alberta::NumberingMap< dimension, Alberta::Dune2AlbertaNumbering > NumberingMap;
91 typedef Alberta::DuneBoundaryProjection< dimension > Projection;
93 typedef std::array< unsigned int, dimension > FaceId;
94 typedef std::map< FaceId, size_t > BoundaryMap;
96 class ProjectionFactory;
100 static const bool supportsBoundaryIds =
true;
102 static const bool supportPeriodicity = MacroData::supportPeriodicity;
106 : globalProjection_( (const DuneProjection *) 0 )
119 macroData_.insertVertex( pos );
128 const std::vector< unsigned int > &vertices )
131 DUNE_THROW( AlbertaError,
"Inserting element of wrong dimension: " << type.dim() );
132 if( !type.isSimplex() )
133 DUNE_THROW( AlbertaError,
"Alberta supports only simplices." );
135 if( vertices.size() != (
size_t)numVertices )
136 DUNE_THROW( AlbertaError,
"Wrong number of vertices passed: " << vertices.size() <<
"." );
138 int array[ numVertices ];
139 for(
int i = 0; i < numVertices; ++i )
140 array[ i ] = vertices[ numberingMap_.alberta2dune(
dimension, i ) ];
141 macroData_.insertElement( array );
154 virtual void insertBoundary (
int element,
int face,
int id )
156 if( (
id <= 0) || (
id > 127) )
157 DUNE_THROW( AlbertaError,
"Invalid boundary id: " <<
id <<
"." );
158 macroData_.boundaryId( element, numberingMap_.dune2alberta( 1, face ) ) = id;
173 const std::vector< unsigned int > &vertices,
174 const DuneProjection *projection )
177 DUNE_THROW( AlbertaError,
"Inserting boundary face of wrong dimension: " << type.dim() );
178 if( !type.isSimplex() )
179 DUNE_THROW( AlbertaError,
"Alberta supports only simplices." );
182 if( vertices.size() != faceId.size() )
183 DUNE_THROW( AlbertaError,
"Wrong number of face vertices passed: " << vertices.size() <<
"." );
184 for(
size_t i = 0; i < faceId.size(); ++i )
185 faceId[ i ] = vertices[ i ];
186 std::sort( faceId.begin(), faceId.end() );
188 typedef std::pair< typename BoundaryMap::iterator, bool > InsertResult;
189 const InsertResult result = boundaryMap_.insert( std::make_pair( faceId, boundaryProjections_.size() ) );
191 DUNE_THROW( GridError,
"Only one boundary projection can be attached to a face." );
192 boundaryProjections_.push_back( DuneProjectionPtr( projection ) );
204 virtual void insertBoundaryProjection (
const DuneProjection *projection )
206 if( globalProjection_ )
207 DUNE_THROW( GridError,
"Only one global boundary projection can be attached to a grid." );
208 globalProjection_ = DuneProjectionPtr( projection );
219 insertBoundaryProjection( GeometryTypes::simplex(
dimension-1 ), vertices, 0 );
229 const std::shared_ptr< BoundarySegment > &boundarySegment )
231 auto refSimplex = ReferenceElements< ctype, dimension-1 >::simplex();
233 if( !boundarySegment )
234 DUNE_THROW( GridError,
"Trying to insert null as a boundary segment." );
235 if( (
int)vertices.size() != refSimplex.size(
dimension-1 ) )
236 DUNE_THROW( GridError,
"Wrong number of face vertices passed: " << vertices.size() <<
"." );
238 std::vector< WorldVector > coords( refSimplex.size(
dimension-1 ) );
241 Alberta::GlobalVector &x = macroData_.vertex( vertices[ i ] );
242 for(
int j = 0; j < dimensionworld; ++j )
243 coords[ i ][ j ] = x[ j ];
244 if( ((*boundarySegment)( refSimplex.position( i,
dimension-1 ) ) - coords[ i ]).two_norm() > 1e-6 )
245 DUNE_THROW( GridError,
"Boundary segment does not interpolate the corners." );
249 const DuneProjection *prj =
new BoundarySegmentWrapper( gt, coords, boundarySegment );
250 insertBoundaryProjection( gt, vertices, prj );
266 void insertFaceTransformation (
const WorldMatrix &matrix,
const WorldVector &shift );
276 void markLongestEdge ()
278 macroData_.markLongestEdge();
295 macroData_.finalize();
296 if( macroData_.elementCount() == 0 )
297 DUNE_THROW( GridError,
"Cannot create empty AlbertaGrid." );
299 macroData_.setOrientation( Alberta::Real( 1 ) );
300 assert( macroData_.checkNeighbors() );
301 macroData_.checkCycles();
302 ProjectionFactory projectionFactory( *
this );
303 return new Grid( macroData_, projectionFactory );
310 static void destroyGrid ( Grid *grid )
321 bool write (
const std::string &filename )
323 macroData_.finalize();
325 macroData_.setOrientation( Alberta::Real( 1 ) );
326 assert( macroData_.checkNeighbors() );
327 return macroData_.write( filename,
false );
339 const int elIndex =
insertionIndex( entity.impl().elementInfo() );
340 const typename MacroData::ElementId &elementId = macroData_.element( elIndex );
341 return elementId[ entity.impl().subEntity() ];
347 const Grid &grid = intersection.impl().grid();
348 const ElementInfo &elementInfo = intersection.impl().elementInfo();
349 const int face = grid.generic2alberta( 1, intersection.indexInInside() );
356 return (
insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
360 unsigned int insertionIndex (
const ElementInfo &elementInfo )
const;
361 unsigned int insertionIndex (
const ElementInfo &elementInfo,
const int face )
const;
363 FaceId faceId (
const ElementInfo &elementInfo,
const int face )
const;
365 MacroData macroData_;
366 NumberingMap numberingMap_;
367 DuneProjectionPtr globalProjection_;
368 BoundaryMap boundaryMap_;
369 std::vector< DuneProjectionPtr > boundaryProjections_;
373 template<
int dim,
int dimworld >
374 GridFactory< AlbertaGrid< dim, dimworld > >::~GridFactory ()
376 macroData_.release();
380 template<
int dim,
int dimworld >
382 GridFactory< AlbertaGrid< dim, dimworld > >
383 ::insertFaceTransformation (
const WorldMatrix &matrix,
const WorldVector &shift )
386 for(
int i = 0; i < dimworld; ++i )
387 for(
int j = 0; j < dimworld; ++j )
389 const ctype delta = (i == j ? ctype( 1 ) : ctype( 0 ));
390 const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
392 if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
394 DUNE_THROW( AlbertaError,
395 "Matrix of face transformation is not orthogonal." );
400 Alberta::GlobalMatrix M;
401 for(
int i = 0; i < dimworld; ++i )
402 for(
int j = 0; j < dimworld; ++j )
403 M[ i ][ j ] = matrix[ i ][ j ];
406 Alberta::GlobalVector t;
407 for(
int i = 0; i < dimworld; ++i )
411 macroData_.insertWallTrafo( M, t );
415 template<
int dim,
int dimworld >
417 GridFactory< AlbertaGrid< dim, dimworld > >
418 ::insertionIndex (
const ElementInfo &elementInfo )
const
420 const MacroElement ¯oElement = elementInfo.macroElement();
421 const unsigned int index = macroElement.index;
424 const typename MacroData::ElementId &elementId = macroData_.element( index );
425 for(
int i = 0; i <= dimension; ++i )
427 const Alberta::GlobalVector &x = macroData_.vertex( elementId[ i ] );
428 const Alberta::GlobalVector &y = macroElement.coordinate( i );
429 for(
int j = 0; j < dimensionworld; ++j )
431 if( x[ j ] != y[ j ] )
432 DUNE_THROW( GridError,
"Vertex in macro element does not coincide with same vertex in macro data structure." );
435 #endif // #ifndef NDEBUG
441 template<
int dim,
int dimworld >
443 GridFactory< AlbertaGrid< dim, dimworld > >
444 ::insertionIndex (
const ElementInfo &elementInfo,
const int face )
const
446 typedef typename BoundaryMap::const_iterator Iterator;
447 const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
448 if( it != boundaryMap_.end() )
451 return std::numeric_limits< unsigned int >::max();
455 template<
int dim,
int dimworld >
456 inline typename GridFactory< AlbertaGrid< dim, dimworld > >::FaceId
457 GridFactory< AlbertaGrid< dim, dimworld > >
458 ::faceId (
const ElementInfo &elementInfo,
const int face )
const
460 const unsigned int index = insertionIndex( elementInfo );
461 const typename MacroData::ElementId &elementId = macroData_.element( index );
464 for(
size_t i = 0; i < faceId.size(); ++i )
466 const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
467 faceId[ i ] = elementId[ k ];
469 std::sort( faceId.begin(), faceId.end() );
478 template<
int dim,
int dimworld >
479 class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
480 :
public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
482 typedef ProjectionFactory This;
483 typedef Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory > Base;
488 typedef typename Base::Projection Projection;
489 typedef typename Base::ElementInfo ElementInfo;
491 typedef typename Projection::Projection DuneProjection;
493 ProjectionFactory(
const GridFactory &gridFactory )
494 : gridFactory_( gridFactory )
497 bool hasProjection (
const ElementInfo &elementInfo,
const int face )
const
499 if( gridFactory().globalProjection_ )
502 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
503 if( index < std::numeric_limits< unsigned int >::max() )
504 return bool( gridFactory().boundaryProjections_[ index ] );
509 bool hasProjection (
const ElementInfo &elementInfo )
const
511 return bool( gridFactory().globalProjection_ );
514 Projection projection (
const ElementInfo &elementInfo,
const int face )
const
516 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
517 if( index < std::numeric_limits< unsigned int >::max() )
519 const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
521 return Projection( projection );
524 assert( gridFactory().globalProjection_ );
525 return Projection( gridFactory().globalProjection_ );
528 Projection projection (
const ElementInfo &elementInfo )
const
530 assert( gridFactory().globalProjection_ );
531 return Projection( gridFactory().globalProjection_ );
545 #endif // #if HAVE_ALBERTA
547 #endif // #ifndef DUNE_ALBERTA_GRIDFACTORY_HH