Dune Core Modules (2.4.1)

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 <limits>
14#include <map>
15
16#include <dune/common/array.hh>
17
18#include <dune/geometry/referenceelements.hh>
19
21
22#include <dune/grid/utility/grapedataioformattypes.hh>
23
25
26#if HAVE_ALBERTA
27
28namespace Dune
29{
30
48 template< int dim, int dimworld >
49 class GridFactory< AlbertaGrid< dim, dimworld > >
50 : public GridFactoryInterface< AlbertaGrid< dim, dimworld > >
51 {
53
54 public:
57
59 typedef typename Grid::ctype ctype;
60
62 static const int dimension = Grid::dimension;
64 static const int dimensionworld = Grid::dimensionworld;
65
70
72 typedef Dune::shared_ptr< const DuneProjection > DuneProjectionPtr;
74
75 template< int codim >
76 struct Codim
77 {
78 typedef typename Grid::template Codim< codim >::Entity Entity;
79 };
80
81 private:
83
84 static const int numVertices
85 = Alberta::NumSubEntities< dimension, dimension >::value;
86
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;
92
93 typedef array< unsigned int, dimension > FaceId;
94 typedef std::map< FaceId, size_t > BoundaryMap;
95
96 class ProjectionFactory;
97
98 public:
100 static const bool supportsBoundaryIds = true;
102 static const bool supportPeriodicity = MacroData::supportPeriodicity;
103
106 : globalProjection_( (const DuneProjection *) 0 )
107 {
108 macroData_.create();
109 }
110
111 virtual ~GridFactory ();
112
117 virtual void insertVertex ( const WorldVector &pos )
118 {
119 macroData_.insertVertex( pos );
120 }
121
127 virtual void insertElement ( const GeometryType &type,
128 const std::vector< unsigned int > &vertices )
129 {
130 if( (int)type.dim() != dimension )
131 DUNE_THROW( AlbertaError, "Inserting element of wrong dimension: " << type.dim() );
132 if( !type.isSimplex() )
133 DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
134
135 if( vertices.size() != (size_t)numVertices )
136 DUNE_THROW( AlbertaError, "Wrong number of vertices passed: " << vertices.size() << "." );
137
138 int array[ numVertices ];
139 for( int i = 0; i < numVertices; ++i )
140 array[ i ] = vertices[ numberingMap_.alberta2dune( dimension, i ) ];
141 macroData_.insertElement( array );
142 }
143
154 virtual void insertBoundary ( int element, int face, int id )
155 {
156 if( (id <= 0) || (id > 127) )
157 DUNE_THROW( AlbertaError, "Invalid boundary id: " << id << "." );
158 macroData_.boundaryId( element, numberingMap_.dune2alberta( 1, face ) ) = id;
159 }
160
171 virtual void
173 const std::vector< unsigned int > &vertices,
174 const DuneProjection *projection )
175 {
176 if( (int)type.dim() != dimension-1 )
177 DUNE_THROW( AlbertaError, "Inserting boundary face of wrong dimension: " << type.dim() );
178 if( !type.isSimplex() )
179 DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
180
181 FaceId faceId;
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() );
187
188 typedef std::pair< typename BoundaryMap::iterator, bool > InsertResult;
189 const InsertResult result = boundaryMap_.insert( std::make_pair( faceId, boundaryProjections_.size() ) );
190 if( !result.second )
191 DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
192 boundaryProjections_.push_back( DuneProjectionPtr( projection ) );
193 }
194
195
204 virtual void insertBoundaryProjection ( const DuneProjection *projection )
205 {
206 if( globalProjection_ )
207 DUNE_THROW( GridError, "Only one global boundary projection can be attached to a grid." );
208 globalProjection_ = DuneProjectionPtr( projection );
209 }
210
216 virtual void
217 insertBoundarySegment ( const std::vector< unsigned int >& vertices )
218 {
219 typedef typename GenericGeometry::SimplexTopology< dimension-1 >::type Topology;
220 insertBoundaryProjection( GeometryType( Topology() ), vertices, 0 );
221 }
222
228 virtual void
229 insertBoundarySegment ( const std::vector< unsigned int > &vertices,
230 const shared_ptr< BoundarySegment > &boundarySegment )
231 {
232 const ReferenceElement< ctype, dimension-1 > &refSimplex
234
235 if( !boundarySegment )
236 DUNE_THROW( GridError, "Trying to insert null as a boundary segment." );
237 if( (int)vertices.size() != refSimplex.size( dimension-1 ) )
238 DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
239
240 std::vector< WorldVector > coords( refSimplex.size( dimension-1 ) );
241 for( int i = 0; i < dimension; ++i )
242 {
243 Alberta::GlobalVector &x = macroData_.vertex( vertices[ i ] );
244 for( int j = 0; j < dimensionworld; ++j )
245 coords[ i ][ j ] = x[ j ];
246 if( ((*boundarySegment)( refSimplex.position( i, dimension-1 ) ) - coords[ i ]).two_norm() > 1e-6 )
247 DUNE_THROW( GridError, "Boundary segment does not interpolate the corners." );
248 }
249
250 const GeometryType gt = refSimplex.type( 0, 0 );
251 const DuneProjection *prj = new BoundarySegmentWrapper( gt, coords, boundarySegment );
252 insertBoundaryProjection( gt, vertices, prj );
253 }
254
268 void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
269
279 {
280 macroData_.markLongestEdge();
281 }
282
296 {
297 macroData_.finalize();
298 if( macroData_.elementCount() == 0 )
299 DUNE_THROW( GridError, "Cannot create empty AlbertaGrid." );
300 if( dimension < 3 )
301 macroData_.setOrientation( Alberta::Real( 1 ) );
302 assert( macroData_.checkNeighbors() );
303 macroData_.checkCycles();
304 ProjectionFactory projectionFactory( *this );
305 return new Grid( macroData_, projectionFactory );
306 }
307
312 static void destroyGrid ( Grid *grid )
313 {
314 delete grid;
315 }
316
325 template< GrapeIOFileFormatType type >
326 bool write ( const std::string &filename )
327 {
328 static_assert( type != pgm, "AlbertaGridFactory: writing pgm format is not supported." );
329 macroData_.finalize();
330 if( dimension < 3 )
331 macroData_.setOrientation( Alberta::Real( 1 ) );
332 assert( macroData_.checkNeighbors() );
333 return macroData_.write( filename, (type == xdr) );
334 }
335
344 virtual bool write ( const std::string &filename )
345 {
346 return write< ascii >( filename );
347 }
348
349 virtual unsigned int
350 insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
351 {
352 return insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
353 }
354
355 virtual unsigned int
356 insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
357 {
358 const int elIndex = insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
359 const typename MacroData::ElementId &elementId = macroData_.element( elIndex );
360 return elementId[ Grid::getRealImplementation( entity ).subEntity() ];
361 }
362
363 virtual unsigned int
364 insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
365 {
366 const Grid &grid = Grid::getRealImplementation( intersection ).grid();
367 const ElementInfo &elementInfo = Grid::getRealImplementation( intersection ).elementInfo();
368 const int face = grid.generic2alberta( 1, intersection.indexInInside() );
369 return insertionIndex( elementInfo, face );
370 }
371
372 virtual bool
373 wasInserted ( const typename Grid::LeafIntersection &intersection ) const
374 {
375 return (insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
376 }
377
378 private:
379 unsigned int insertionIndex ( const ElementInfo &elementInfo ) const;
380 unsigned int insertionIndex ( const ElementInfo &elementInfo, const int face ) const;
381
382 FaceId faceId ( const ElementInfo &elementInfo, const int face ) const;
383
384 MacroData macroData_;
385 NumberingMap numberingMap_;
386 DuneProjectionPtr globalProjection_;
387 BoundaryMap boundaryMap_;
388 std::vector< DuneProjectionPtr > boundaryProjections_;
389 };
390
391
392 template< int dim, int dimworld >
393 GridFactory< AlbertaGrid< dim, dimworld > >::~GridFactory ()
394 {
395 macroData_.release();
396 }
397
398
399 template< int dim, int dimworld >
400 inline void
401 GridFactory< AlbertaGrid< dim, dimworld > >
402 ::insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift )
403 {
404 // make sure the matrix is orthogonal
405 for( int i = 0; i < dimworld; ++i )
406 for( int j = 0; j < dimworld; ++j )
407 {
408 const ctype delta = (i == j ? ctype( 1 ) : ctype( 0 ));
409 const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
410
411 if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
412 {
413 DUNE_THROW( AlbertaError,
414 "Matrix of face transformation is not orthogonal." );
415 }
416 }
417
418 // copy matrix
419 Alberta::GlobalMatrix M;
420 for( int i = 0; i < dimworld; ++i )
421 for( int j = 0; j < dimworld; ++j )
422 M[ i ][ j ] = matrix[ i ][ j ];
423
424 // copy shift
425 Alberta::GlobalVector t;
426 for( int i = 0; i < dimworld; ++i )
427 t[ i ] = shift[ i ];
428
429 // insert into ALBERTA macro data
430 macroData_.insertWallTrafo( M, t );
431 }
432
433
434 template< int dim, int dimworld >
435 inline unsigned int
437 ::insertionIndex ( const ElementInfo &elementInfo ) const
438 {
439 const MacroElement &macroElement = elementInfo.macroElement();
440 const unsigned int index = macroElement.index;
441
442#ifndef NDEBUG
443 const typename MacroData::ElementId &elementId = macroData_.element( index );
444 for( int i = 0; i <= dimension; ++i )
445 {
446 const Alberta::GlobalVector &x = macroData_.vertex( elementId[ i ] );
447 const Alberta::GlobalVector &y = macroElement.coordinate( i );
448 for( int j = 0; j < dimensionworld; ++j )
449 {
450 if( x[ j ] != y[ j ] )
451 DUNE_THROW( GridError, "Vertex in macro element does not coincide with same vertex in macro data structure." );
452 }
453 }
454#endif // #ifndef NDEBUG
455
456 return index;
457 }
458
459
460 template< int dim, int dimworld >
461 inline unsigned int
462 GridFactory< AlbertaGrid< dim, dimworld > >
463 ::insertionIndex ( const ElementInfo &elementInfo, const int face ) const
464 {
465 typedef typename BoundaryMap::const_iterator Iterator;
466 const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
467 if( it != boundaryMap_.end() )
468 return it->second;
469 else
470 return std::numeric_limits< unsigned int >::max();
471 }
472
473
474 template< int dim, int dimworld >
475 inline typename GridFactory< AlbertaGrid< dim, dimworld > >::FaceId
476 GridFactory< AlbertaGrid< dim, dimworld > >
477 ::faceId ( const ElementInfo &elementInfo, const int face ) const
478 {
479 const unsigned int index = insertionIndex( elementInfo );
480 const typename MacroData::ElementId &elementId = macroData_.element( index );
481
482 FaceId faceId;
483 for( size_t i = 0; i < faceId.size(); ++i )
484 {
485 const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
486 faceId[ i ] = elementId[ k ];
487 }
488 std::sort( faceId.begin(), faceId.end() );
489 return faceId;
490 }
491
492
493
494 // GridFactory::ProjectionFactory
495 // ------------------------------
496
497 template< int dim, int dimworld >
498 class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
499 : public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
500 {
501 typedef ProjectionFactory This;
502 typedef Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory > Base;
503
504 typedef typename Dune::GridFactory< AlbertaGrid< dim, dimworld > > Factory;
505
506 public:
507 typedef typename Base::Projection Projection;
508 typedef typename Base::ElementInfo ElementInfo;
509
510 typedef typename Projection::Projection DuneProjection;
511
512 ProjectionFactory( const GridFactory &gridFactory )
513 : gridFactory_( gridFactory )
514 {}
515
516 bool hasProjection ( const ElementInfo &elementInfo, const int face ) const
517 {
518 if( gridFactory().globalProjection_ )
519 return true;
520
521 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
522 if( index < std::numeric_limits< unsigned int >::max() )
523 return bool( gridFactory().boundaryProjections_[ index ] );
524 else
525 return false;
526 }
527
528 bool hasProjection ( const ElementInfo &elementInfo ) const
529 {
530 return bool( gridFactory().globalProjection_ );
531 }
532
533 Projection projection ( const ElementInfo &elementInfo, const int face ) const
534 {
535 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
536 if( index < std::numeric_limits< unsigned int >::max() )
537 {
538 const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
539 if( projection )
540 return Projection( projection );
541 }
542
543 assert( gridFactory().globalProjection_ );
544 return Projection( gridFactory().globalProjection_ );
545 };
546
547 Projection projection ( const ElementInfo &elementInfo ) const
548 {
549 assert( gridFactory().globalProjection_ );
550 return Projection( gridFactory().globalProjection_ );
551 };
552
553 const GridFactory &gridFactory () const
554 {
555 return gridFactory_;
556 }
557
558 private:
559 const GridFactory &gridFactory_;
560 };
561
562}
563
564#endif // #if HAVE_ALBERTA
565
566#endif // #ifndef DUNE_ALBERTA_GRIDFACTORY_HH
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:62
A dense n x m matrix.
Definition: fmatrix.hh:67
vector space out of a tensor product of fields.
Definition: fvector.hh:94
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:25
unsigned int dim() const
Return dimension of the type.
Definition: type.hh:321
bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:306
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:350
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:229
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:172
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:278
bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition: gridfactory.hh:326
virtual bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition: gridfactory.hh:344
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:154
GridFactory()
Definition: gridfactory.hh:105
static void destroyGrid(Grid *grid)
destroy a grid previously obtain from this factory
Definition: gridfactory.hh:312
virtual unsigned int insertionIndex(const typename Codim< dimension >::Entity &entity) const
obtain a vertex' insertion index
Definition: gridfactory.hh:356
virtual void insertBoundaryProjection(const DuneProjection *projection)
insert a global (boundary) projection into the macro grid
Definition: gridfactory.hh:204
Grid * createGrid()
finalize grid creation and hand over the grid
Definition: gridfactory.hh:295
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment into the macro grid
Definition: gridfactory.hh:217
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:388
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:486
@ dimensionworld
The dimension of the world the grid lives in.
Definition: grid.hh:408
@ dimension
The dimension of the grid.
Definition: grid.hh:402
This class provides access to geometric and topological properties of a reference element.
Definition: referenceelements.hh:55
Provide a generic factory class for unstructured grids.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
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:10
@ xdr
Definition: grapedataioformattypes.hh:16
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:30
Static tag representing a codimension.
Definition: dimension.hh:22
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:490
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 22, 23:33, 2024)