dune-grid 3.0-git
albertagrid/gridfactory.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4#ifndef DUNE_ALBERTA_GRIDFACTORY_HH
5#define DUNE_ALBERTA_GRIDFACTORY_HH
6
12#include <algorithm>
13#include <array>
14#include <limits>
15#include <map>
16#include <memory>
17
18#include <dune/geometry/referenceelements.hh>
19
21
23
24#if HAVE_ALBERTA
25
26namespace Dune
27{
28
46 template< int dim, int dimworld >
47 class GridFactory< AlbertaGrid< dim, dimworld > >
48 : public GridFactoryInterface< AlbertaGrid< dim, dimworld > >
49 {
51
52 public:
55
57 typedef typename Grid::ctype ctype;
58
60 static const int dimension = Grid::dimension;
62 static const int dimensionworld = Grid::dimensionworld;
63
65 typedef FieldVector< ctype, dimensionworld > WorldVector;
67 typedef FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix;
68
70 typedef std::shared_ptr< const DuneProjection > DuneProjectionPtr;
72
73 template< int codim >
74 struct Codim
75 {
76 typedef typename Grid::template Codim< codim >::Entity Entity;
77 };
78
79 private:
81
82 static const int numVertices
84
90
91 typedef std::array< unsigned int, dimension > FaceId;
92 typedef std::map< FaceId, size_t > BoundaryMap;
93
94 class ProjectionFactory;
95
96 public:
98 static const bool supportsBoundaryIds = true;
100 static const bool supportPeriodicity = MacroData::supportPeriodicity;
101
104 : globalProjection_( (const DuneProjection *) 0 )
105 {
106 macroData_.create();
107 }
108
109 virtual ~GridFactory ();
110
115 virtual void insertVertex ( const WorldVector &pos )
116 {
117 macroData_.insertVertex( pos );
118 }
119
125 virtual void insertElement ( const GeometryType &type,
126 const std::vector< unsigned int > &vertices )
127 {
128 if( (int)type.dim() != dimension )
129 DUNE_THROW( AlbertaError, "Inserting element of wrong dimension: " << type.dim() );
130 if( !type.isSimplex() )
131 DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
132
133 if( vertices.size() != (size_t)numVertices )
134 DUNE_THROW( AlbertaError, "Wrong number of vertices passed: " << vertices.size() << "." );
135
136 int array[ numVertices ];
137 for( int i = 0; i < numVertices; ++i )
138 array[ i ] = vertices[ numberingMap_.alberta2dune( dimension, i ) ];
139 macroData_.insertElement( array );
140 }
141
152 virtual void insertBoundary ( int element, int face, int id )
153 {
154 if( (id <= 0) || (id > 127) )
155 DUNE_THROW( AlbertaError, "Invalid boundary id: " << id << "." );
156 macroData_.boundaryId( element, numberingMap_.dune2alberta( 1, face ) ) = id;
157 }
158
169 virtual void
170 insertBoundaryProjection ( const GeometryType &type,
171 const std::vector< unsigned int > &vertices,
172 const DuneProjection *projection )
173 {
174 if( (int)type.dim() != dimension-1 )
175 DUNE_THROW( AlbertaError, "Inserting boundary face of wrong dimension: " << type.dim() );
176 if( !type.isSimplex() )
177 DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
178
179 FaceId faceId;
180 if( vertices.size() != faceId.size() )
181 DUNE_THROW( AlbertaError, "Wrong number of face vertices passed: " << vertices.size() << "." );
182 for( size_t i = 0; i < faceId.size(); ++i )
183 faceId[ i ] = vertices[ i ];
184 std::sort( faceId.begin(), faceId.end() );
185
186 typedef std::pair< typename BoundaryMap::iterator, bool > InsertResult;
187 const InsertResult result = boundaryMap_.insert( std::make_pair( faceId, boundaryProjections_.size() ) );
188 if( !result.second )
189 DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
190 boundaryProjections_.push_back( DuneProjectionPtr( projection ) );
191 }
192
193
202 virtual void insertBoundaryProjection ( const DuneProjection *projection )
203 {
204 if( globalProjection_ )
205 DUNE_THROW( GridError, "Only one global boundary projection can be attached to a grid." );
206 globalProjection_ = DuneProjectionPtr( projection );
207 }
208
214 virtual void
215 insertBoundarySegment ( const std::vector< unsigned int >& vertices )
216 {
217 typedef typename GenericGeometry::SimplexTopology< dimension-1 >::type Topology;
218 insertBoundaryProjection( GeometryType( Topology() ), vertices, 0 );
219 }
220
226 virtual void
227 insertBoundarySegment ( const std::vector< unsigned int > &vertices,
228 const std::shared_ptr< BoundarySegment > &boundarySegment )
229 {
230 const ReferenceElement< ctype, dimension-1 > &refSimplex
231 = ReferenceElements< ctype, dimension-1 >::simplex();
232
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() << "." );
237
238 std::vector< WorldVector > coords( refSimplex.size( dimension-1 ) );
239 for( int i = 0; i < dimension; ++i )
240 {
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." );
246 }
247
248 const GeometryType gt = refSimplex.type( 0, 0 );
249 const DuneProjection *prj = new BoundarySegmentWrapper( gt, coords, boundarySegment );
250 insertBoundaryProjection( gt, vertices, prj );
251 }
252
266 void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
267
277 {
278 macroData_.markLongestEdge();
279 }
280
294 {
295 macroData_.finalize();
296 if( macroData_.elementCount() == 0 )
297 DUNE_THROW( GridError, "Cannot create empty AlbertaGrid." );
298 if( dimension < 3 )
299 macroData_.setOrientation( Alberta::Real( 1 ) );
300 assert( macroData_.checkNeighbors() );
301 macroData_.checkCycles();
302 ProjectionFactory projectionFactory( *this );
303 return new Grid( macroData_, projectionFactory );
304 }
305
310 static void destroyGrid ( Grid *grid )
311 {
312 delete grid;
313 }
314
321 bool write ( const std::string &filename )
322 {
323 macroData_.finalize();
324 if( dimension < 3 )
325 macroData_.setOrientation( Alberta::Real( 1 ) );
326 assert( macroData_.checkNeighbors() );
327 return macroData_.write( filename, false );
328 }
329
330 virtual unsigned int
331 insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
332 {
333 return insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
334 }
335
336 virtual unsigned int
337 insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
338 {
339 const int elIndex = insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
340 const typename MacroData::ElementId &elementId = macroData_.element( elIndex );
341 return elementId[ Grid::getRealImplementation( entity ).subEntity() ];
342 }
343
344 virtual unsigned int
345 insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
346 {
347 const Grid &grid = Grid::getRealImplementation( intersection ).grid();
348 const ElementInfo &elementInfo = Grid::getRealImplementation( intersection ).elementInfo();
349 const int face = grid.generic2alberta( 1, intersection.indexInInside() );
350 return insertionIndex( elementInfo, face );
351 }
352
353 virtual bool
354 wasInserted ( const typename Grid::LeafIntersection &intersection ) const
355 {
356 return (insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
357 }
358
359 private:
360 unsigned int insertionIndex ( const ElementInfo &elementInfo ) const;
361 unsigned int insertionIndex ( const ElementInfo &elementInfo, const int face ) const;
362
363 FaceId faceId ( const ElementInfo &elementInfo, const int face ) const;
364
365 MacroData macroData_;
366 NumberingMap numberingMap_;
367 DuneProjectionPtr globalProjection_;
368 BoundaryMap boundaryMap_;
369 std::vector< DuneProjectionPtr > boundaryProjections_;
370 };
371
372
373 template< int dim, int dimworld >
375 {
376 macroData_.release();
377 }
378
379
380 template< int dim, int dimworld >
381 inline void
383 ::insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift )
384 {
385 // make sure the matrix is orthogonal
386 for( int i = 0; i < dimworld; ++i )
387 for( int j = 0; j < dimworld; ++j )
388 {
389 const ctype delta = (i == j ? ctype( 1 ) : ctype( 0 ));
390 const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
391
392 if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
393 {
394 DUNE_THROW( AlbertaError,
395 "Matrix of face transformation is not orthogonal." );
396 }
397 }
398
399 // copy matrix
401 for( int i = 0; i < dimworld; ++i )
402 for( int j = 0; j < dimworld; ++j )
403 M[ i ][ j ] = matrix[ i ][ j ];
404
405 // copy shift
407 for( int i = 0; i < dimworld; ++i )
408 t[ i ] = shift[ i ];
409
410 // insert into ALBERTA macro data
411 macroData_.insertWallTrafo( M, t );
412 }
413
414
415 template< int dim, int dimworld >
416 inline unsigned int
418 ::insertionIndex ( const ElementInfo &elementInfo ) const
419 {
420 const MacroElement &macroElement = elementInfo.macroElement();
421 const unsigned int index = macroElement.index;
422
423#ifndef NDEBUG
424 const typename MacroData::ElementId &elementId = macroData_.element( index );
425 for( int i = 0; i <= dimension; ++i )
426 {
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 )
430 {
431 if( x[ j ] != y[ j ] )
432 DUNE_THROW( GridError, "Vertex in macro element does not coincide with same vertex in macro data structure." );
433 }
434 }
435#endif // #ifndef NDEBUG
436
437 return index;
438 }
439
440
441 template< int dim, int dimworld >
442 inline unsigned int
443 GridFactory< AlbertaGrid< dim, dimworld > >
444 ::insertionIndex ( const ElementInfo &elementInfo, const int face ) const
445 {
446 typedef typename BoundaryMap::const_iterator Iterator;
447 const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
448 if( it != boundaryMap_.end() )
449 return it->second;
450 else
451 return std::numeric_limits< unsigned int >::max();
452 }
453
454
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
459 {
460 const unsigned int index = insertionIndex( elementInfo );
461 const typename MacroData::ElementId &elementId = macroData_.element( index );
462
463 FaceId faceId;
464 for( size_t i = 0; i < faceId.size(); ++i )
465 {
466 const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
467 faceId[ i ] = elementId[ k ];
468 }
469 std::sort( faceId.begin(), faceId.end() );
470 return faceId;
471 }
472
473
474
475 // GridFactory::ProjectionFactory
476 // ------------------------------
477
478 template< int dim, int dimworld >
479 class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
480 : public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
481 {
482 typedef ProjectionFactory This;
484
486
487 public:
488 typedef typename Base::Projection Projection;
490
492
493 ProjectionFactory( const GridFactory &gridFactory )
494 : gridFactory_( gridFactory )
495 {}
496
497 bool hasProjection ( const ElementInfo &elementInfo, const int face ) const
498 {
499 if( gridFactory().globalProjection_ )
500 return true;
501
502 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
503 if( index < std::numeric_limits< unsigned int >::max() )
504 return bool( gridFactory().boundaryProjections_[ index ] );
505 else
506 return false;
507 }
508
509 bool hasProjection ( const ElementInfo &elementInfo ) const
510 {
511 return bool( gridFactory().globalProjection_ );
512 }
513
514 Projection projection ( const ElementInfo &elementInfo, const int face ) const
515 {
516 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
517 if( index < std::numeric_limits< unsigned int >::max() )
518 {
519 const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
520 if( projection )
521 return Projection( projection );
522 }
523
524 assert( gridFactory().globalProjection_ );
525 return Projection( gridFactory().globalProjection_ );
526 };
527
528 Projection projection ( const ElementInfo &elementInfo ) const
529 {
530 assert( gridFactory().globalProjection_ );
531 return Projection( gridFactory().globalProjection_ );
532 };
533
534 const GridFactory &gridFactory () const
535 {
536 return gridFactory_;
537 }
538
539 private:
540 const GridFactory &gridFactory_;
541 };
542
543}
544
545#endif // #if HAVE_ALBERTA
546
547#endif // #ifndef DUNE_ALBERTA_GRIDFACTORY_HH
provides the AlbertaGrid class
Include standard header files.
Definition agrid.hh:60
ALBERTA REAL_DD GlobalMatrix
Definition misc.hh:48
ALBERTA REAL Real
Definition misc.hh:45
ALBERTA REAL_D GlobalVector
Definition misc.hh:47
[ provides Dune::Grid ]
Definition agrid.hh:140
int generic2alberta(int codim, int i) const
Definition agrid.hh:520
GridFamily::ctype ctype
Definition agrid.hh:175
specialization of the generic GridFactory for AlbertaGrid
Definition albertagrid/gridfactory.hh:49
DuneBoundaryProjection< dimensionworld > DuneProjection
Definition albertagrid/gridfactory.hh:69
FieldVector< ctype, dimensionworld > WorldVector
type of vector for world coordinates
Definition albertagrid/gridfactory.hh:65
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition albertagrid/gridfactory.hh:331
virtual bool wasInserted(const typename Grid::LeafIntersection &intersection) const
Definition albertagrid/gridfactory.hh:354
std::shared_ptr< const DuneProjection > DuneProjectionPtr
Definition albertagrid/gridfactory.hh:70
AlbertaGrid< dim, dimworld > Grid
type of grid this factory is for
Definition albertagrid/gridfactory.hh:54
virtual void insertBoundaryProjection(const GeometryType &type, const std::vector< unsigned int > &vertices, const DuneProjection *projection)
insert a boundary projection into the macro grid
Definition albertagrid/gridfactory.hh:170
FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix
type of matrix from world coordinates to world coordinates
Definition albertagrid/gridfactory.hh:67
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
insert an element into the macro grid
Definition albertagrid/gridfactory.hh:125
Dune::BoundarySegment< dimension, dimensionworld > BoundarySegment
Definition albertagrid/gridfactory.hh:71
void markLongestEdge()
mark the longest edge as refinemet edge
Definition albertagrid/gridfactory.hh:276
Grid::ctype ctype
type of (scalar) coordinates
Definition albertagrid/gridfactory.hh:57
virtual void insertVertex(const WorldVector &pos)
insert a vertex into the macro grid
Definition albertagrid/gridfactory.hh:115
virtual void insertBoundary(int element, int face, int id)
mark a face as boundary (and assign a boundary id)
Definition albertagrid/gridfactory.hh:152
GridFactory()
Definition albertagrid/gridfactory.hh:103
static void destroyGrid(Grid *grid)
destroy a grid previously obtain from this factory
Definition albertagrid/gridfactory.hh:310
virtual unsigned int insertionIndex(const typename Codim< dimension >::Entity &entity) const
obtain a vertex' insertion index
Definition albertagrid/gridfactory.hh:337
bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition albertagrid/gridfactory.hh:321
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices, const std::shared_ptr< BoundarySegment > &boundarySegment)
insert a shaped boundary segment into the macro grid
Definition albertagrid/gridfactory.hh:227
virtual void insertBoundaryProjection(const DuneProjection *projection)
insert a global (boundary) projection into the macro grid
Definition albertagrid/gridfactory.hh:202
Grid * createGrid()
finalize grid creation and hand over the grid
Definition albertagrid/gridfactory.hh:293
virtual unsigned int insertionIndex(const typename Grid::LeafIntersection &intersection) const
Definition albertagrid/gridfactory.hh:345
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment into the macro grid
Definition albertagrid/gridfactory.hh:215
Grid::template Codim< codim >::Entity Entity
Definition albertagrid/gridfactory.hh:76
Projection projection(const ElementInfo &elementInfo) const
Definition albertagrid/gridfactory.hh:528
bool hasProjection(const ElementInfo &elementInfo) const
Definition albertagrid/gridfactory.hh:509
Base::ElementInfo ElementInfo
Definition albertagrid/gridfactory.hh:489
Base::Projection Projection
Definition albertagrid/gridfactory.hh:488
const GridFactory & gridFactory() const
Definition albertagrid/gridfactory.hh:534
ProjectionFactory(const GridFactory &gridFactory)
Definition albertagrid/gridfactory.hh:493
Projection projection(const ElementInfo &elementInfo, const int face) const
Definition albertagrid/gridfactory.hh:514
Projection::Projection DuneProjection
Definition albertagrid/gridfactory.hh:491
bool hasProjection(const ElementInfo &elementInfo, const int face) const
Definition albertagrid/gridfactory.hh:497
Definition macroelement.hh:22
Definition misc.hh:29
Definition misc.hh:145
Definition albertagrid/projection.hh:133
Definition albertagrid/projection.hh:34
Interface class for vertex projection at the boundary.
Definition boundaryprojection.hh:24
Definition boundaryprojection.hh:67
Base class for classes implementing geometries of boundary segments.
Definition boundarysegment.hh:30
Base class for exceptions in Dune grid modules.
Definition exceptions.hh:18
GridFamily::Traits::LeafIntersection LeafIntersection
A type that is a model of Dune::Intersection, an intersections of two codimension 1 of two codimensio...
Definition common/grid.hh:460
@ dimension
The dimension of the grid.
Definition common/grid.hh:387
@ dimensionworld
The dimension of the world the grid lives in.
Definition common/grid.hh:393
Provide a generic factory class for unstructured grids.
Definition common/gridfactory.hh:74
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition common/gridfactory.hh:179
static const int dimension
dimension of the grid
Definition common/gridfactory.hh:78
Provide a generic factory class for unstructured grids.
Definition common/gridfactory.hh:263
Provide a generic factory class for unstructured grids.