4#ifndef DUNE_ALBERTA_GRIDFACTORY_HH
5#define DUNE_ALBERTA_GRIDFACTORY_HH
18#include <dune/geometry/referenceelements.hh>
22#include <dune/grid/utility/grapedataioformattypes.hh>
48 template<
int dim,
int dimworld >
59 typedef typename Grid::ctype
ctype;
78 typedef typename Grid::template Codim< codim >::Entity
Entity;
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;
94 typedef std::map< FaceId, size_t > BoundaryMap;
96 class ProjectionFactory;
100 static const bool supportsBoundaryIds =
true;
102 static const bool supportPeriodicity = MacroData::supportPeriodicity;
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 if( (
id <= 0) || (
id > 127) )
155 DUNE_THROW( AlbertaError,
"Invalid boundary id: " <<
id <<
"." );
156 macroData_.boundaryId( element, numberingMap_.dune2alberta( 1, face ) ) = id;
169 const std::vector< unsigned int > &vertices,
173 DUNE_THROW( AlbertaError,
"Inserting boundary face of wrong dimension: " << type.dim() );
174 if( !type.isSimplex() )
175 DUNE_THROW( AlbertaError,
"Alberta supports only simplices." );
178 if( vertices.size() != faceId.
size() )
179 DUNE_THROW( AlbertaError,
"Wrong number of face vertices passed: " << vertices.size() <<
"." );
180 for(
size_t i = 0; i < faceId.
size(); ++i )
181 faceId[ i ] = vertices[ i ];
182 std::sort( faceId.begin(), faceId.end() );
184 typedef std::pair< typename BoundaryMap::iterator, bool > InsertResult;
185 const InsertResult result = boundaryMap_.insert( std::make_pair( faceId, boundaryProjections_.
size() ) );
200 if( globalProjection_ )
213 typedef typename GenericGeometry::SimplexTopology<
dimension-1 >::type Topology;
214 insertBoundaryProjection(
GeometryType( Topology() ), vertices, 0 );
229 if( !boundarySegment )
231 if( (
int)vertices.size() != refSimplex.size(
dimension-1 ) )
232 DUNE_THROW(
GridError,
"Wrong number of face vertices passed: " << vertices.size() <<
"." );
234 std::vector< WorldVector > coords( refSimplex.size(
dimension-1 ) );
237 Alberta::GlobalVector &x = macroData_.vertex( vertices[ i ] );
238 for(
int j = 0; j < dimensionworld; ++j )
239 coords[ i ][ j ] = x[ j ];
240 if( ((*boundarySegment)( refSimplex.position( i,
dimension-1 ) ) - coords[ i ]).two_norm() > 1e-6 )
246 insertBoundaryProjection(
gt, vertices, prj );
262 void insertFaceTransformation (
const WorldMatrix &matrix,
const WorldVector &shift );
274 macroData_.markLongestEdge();
291 macroData_.finalize();
292 if( macroData_.elementCount() == 0 )
295 macroData_.setOrientation( Alberta::Real( 1 ) );
296 assert( macroData_.checkNeighbors() );
297 macroData_.checkCycles();
298 ProjectionFactory projectionFactory( *
this );
299 return new Grid( macroData_, projectionFactory );
319 template< GrapeIOFileFormatType type >
320 bool write (
const std::string &filename )
322 dune_static_assert( type != pgm,
"AlbertaGridFactory: writing pgm format is not supported." );
323 macroData_.finalize();
325 macroData_.setOrientation( Alberta::Real( 1 ) );
326 assert( macroData_.checkNeighbors() );
327 return macroData_.write( filename, (type ==
xdr) );
338 virtual bool write (
const std::string &filename )
340 return write< ascii >( filename );
346 return insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
352 const int elIndex =
insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
353 const typename MacroData::ElementId &elementId = macroData_.element( elIndex );
354 return elementId[ Grid::getRealImplementation( entity ).subEntity() ];
360 const Grid &grid = Grid::getRealImplementation( intersection ).grid();
361 const ElementInfo &elementInfo = Grid::getRealImplementation( intersection ).elementInfo();
362 const int face = grid.generic2alberta( 1, intersection.indexInInside() );
369 return (
insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
373 unsigned int insertionIndex (
const ElementInfo &elementInfo )
const;
374 unsigned int insertionIndex (
const ElementInfo &elementInfo,
const int face )
const;
376 FaceId faceId (
const ElementInfo &elementInfo,
const int face )
const;
378 MacroData macroData_;
379 NumberingMap numberingMap_;
380 DuneProjectionPtr globalProjection_;
381 BoundaryMap boundaryMap_;
382 std::vector< DuneProjectionPtr > boundaryProjections_;
386 template<
int dim,
int dimworld >
387 GridFactory< AlbertaGrid< dim, dimworld > >::~GridFactory ()
389 macroData_.release();
393 template<
int dim,
int dimworld >
395 GridFactory< AlbertaGrid< dim, dimworld > >
399 for(
int i = 0; i < dimworld; ++i )
400 for(
int j = 0; j < dimworld; ++j )
403 const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
405 if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
408 "Matrix of face transformation is not orthogonal." );
413 Alberta::GlobalMatrix M;
414 for(
int i = 0; i < dimworld; ++i )
415 for(
int j = 0; j < dimworld; ++j )
416 M[ i ][ j ] = matrix[ i ][ j ];
419 Alberta::GlobalVector t;
420 for(
int i = 0; i < dimworld; ++i )
424 macroData_.insertWallTrafo( M, t );
428 template<
int dim,
int dimworld >
431 ::insertionIndex (
const ElementInfo &elementInfo )
const
433 const MacroElement ¯oElement = elementInfo.macroElement();
434 const unsigned int index = macroElement.index;
437 const typename MacroData::ElementId &elementId = macroData_.element( index );
438 for(
int i = 0; i <= dimension; ++i )
440 const Alberta::GlobalVector &x = macroData_.vertex( elementId[ i ] );
441 const Alberta::GlobalVector &y = macroElement.coordinate( i );
442 for(
int j = 0; j < dimensionworld; ++j )
444 if( x[ j ] != y[ j ] )
445 DUNE_THROW(
GridError,
"Vertex in macro element does not coincide with same vertex in macro data structure." );
454 template<
int dim,
int dimworld >
456 GridFactory< AlbertaGrid< dim, dimworld > >
457 ::insertionIndex (
const ElementInfo &elementInfo,
const int face )
const
459 typedef typename BoundaryMap::const_iterator Iterator;
460 const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
461 if( it != boundaryMap_.end() )
464 return std::numeric_limits< unsigned int >::max();
468 template<
int dim,
int dimworld >
469 inline typename GridFactory< AlbertaGrid< dim, dimworld > >::FaceId
470 GridFactory< AlbertaGrid< dim, dimworld > >
471 ::faceId (
const ElementInfo &elementInfo,
const int face )
const
473 const unsigned int index = insertionIndex( elementInfo );
474 const typename MacroData::ElementId &elementId = macroData_.element( index );
477 for(
size_t i = 0; i < faceId.size(); ++i )
479 const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
480 faceId[ i ] = elementId[ k ];
482 std::sort( faceId.begin(), faceId.end() );
491 template<
int dim,
int dimworld >
492 class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
493 :
public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
495 typedef ProjectionFactory This;
496 typedef Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory > Base;
501 typedef typename Base::Projection Projection;
502 typedef typename Base::ElementInfo ElementInfo;
504 typedef typename Projection::Projection DuneProjection;
506 ProjectionFactory(
const GridFactory &gridFactory )
507 : gridFactory_( gridFactory )
510 bool hasProjection (
const ElementInfo &elementInfo,
const int face )
const
512 if( gridFactory().globalProjection_ )
515 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
516 if( index < std::numeric_limits< unsigned int >::max() )
517 return bool( gridFactory().boundaryProjections_[ index ] );
522 bool hasProjection (
const ElementInfo &elementInfo )
const
524 return bool( gridFactory().globalProjection_ );
527 Projection projection (
const ElementInfo &elementInfo,
const int face )
const
529 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
530 if( index < std::numeric_limits< unsigned int >::max() )
532 const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
534 return Projection( projection );
537 assert( gridFactory().globalProjection_ );
538 return Projection( gridFactory().globalProjection_ );
541 Projection projection (
const ElementInfo &elementInfo )
const
543 assert( gridFactory().globalProjection_ );
544 return Projection( gridFactory().globalProjection_ );
provides the AlbertaGrid class
Fallback implementation of the std::array class (a static array)
[ provides Dune::Grid ]
Definition: agrid.hh:140
Definition: boundaryprojection.hh:67
Wrapper class for entities.
Definition: entity.hh:57
A dense n x m matrix.
Definition: fmatrix.hh:67
vector space out of a tensor product of fields.
Definition: fvector.hh:92
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:25
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:18
Provide a generic factory class for unstructured grids.
Definition: gridfactory.hh:74
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition: gridfactory.hh:179
static const int dimension
dimension of the grid
Definition: gridfactory.hh:78
virtual bool wasInserted(const typename GridType::LeafIntersection &intersection) const
determine whether an intersection was inserted
Definition: gridfactory.hh:245
specialization of the generic GridFactory for AlbertaGrid
Definition: gridfactory.hh:51
FieldVector< ctype, dimensionworld > WorldVector
type of vector for world coordinates
Definition: gridfactory.hh:67
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition: gridfactory.hh:344
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices, const shared_ptr< BoundarySegment > &boundarySegment)
insert a shaped boundary segment into the macro grid
Definition: gridfactory.hh:223
AlbertaGrid< dim, dimworld > Grid
type of grid this factory is for
Definition: gridfactory.hh:56
virtual void insertBoundaryProjection(const GeometryType &type, const std::vector< unsigned int > &vertices, const DuneProjection *projection)
insert a boundary projection into the macro grid
Definition: gridfactory.hh:168
FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix
type of matrix from world coordinates to world coordinates
Definition: gridfactory.hh:69
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
insert an element into the macro grid
Definition: gridfactory.hh:127
void markLongestEdge()
mark the longest edge as refinemet edge
Definition: gridfactory.hh:272
bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition: gridfactory.hh:320
virtual bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition: gridfactory.hh:338
Grid::ctype ctype
type of (scalar) coordinates
Definition: gridfactory.hh:59
virtual void insertVertex(const WorldVector &pos)
insert a vertex into the macro grid
Definition: gridfactory.hh:117
virtual void insertBoundary(int element, int face, int id)
mark a face as boundary (and assign a boundary id)
Definition: gridfactory.hh:152
GridFactory()
Definition: gridfactory.hh:105
static void destroyGrid(Grid *grid)
destroy a grid previously obtain from this factory
Definition: gridfactory.hh:306
virtual unsigned int insertionIndex(const typename Codim< dimension >::Entity &entity) const
obtain a vertex' insertion index
Definition: gridfactory.hh:350
virtual void insertBoundaryProjection(const DuneProjection *projection)
insert a global (boundary) projection into the macro grid
Definition: gridfactory.hh:198
Grid * createGrid()
finalize grid creation and hand over the grid
Definition: gridfactory.hh:289
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment into the macro grid
Definition: gridfactory.hh:211
Provide a generic factory class for unstructured grids.
Definition: gridfactory.hh:263
GridFactory()
Default constructor.
Definition: gridfactory.hh:274
Grid abstract base class.
Definition: grid.hh:386
@ dimension
The dimension of the grid.
Definition: grid.hh:400
GridFamily::Traits::LeafIntersection LeafIntersection
A type that is a model of Dune::Intersection, an intersections of two codimension 1 of two codimensio...
Definition: grid.hh:484
@ dimensionworld
The dimension of the world the grid lives in.
Definition: grid.hh:406
This class provides access to geometric and topological properties of a reference element....
Definition: referenceelements.hh:58
Simple fixed size array class. This replaces std::array, if that is not available.
Definition: array.hh:40
size_type size() const
Return array size.
Definition: array.hh:71
Provide a generic factory class for unstructured grids.
#define dune_static_assert(COND, MSG)
Helper template so that compilation fails if condition is not true.
Definition: static_assert.hh:79
#define DUNE_THROW(E, m)
Definition: exceptions.hh:244
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:132
Dune namespace.
Definition: alignment.hh:14
@ xdr
Definition: grapedataioformattypes.hh:16
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:30
Interface class for vertex projection at the boundary.
Definition: boundaryprojection.hh:24
static const ReferenceElement< ctype, dim > & simplex()
get simplex reference elements
Definition: referenceelements.hh:574