Dune Core Modules (2.8.0)

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
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
83 = Alberta::NumSubEntities< dimension, dimension >::value;
84
85 typedef Alberta::MacroElement< dimension > MacroElement;
86 typedef Alberta::ElementInfo< dimension > ElementInfo;
87 typedef Alberta::MacroData< dimension > MacroData;
88 typedef Alberta::NumberingMap< dimension, Alberta::Dune2AlbertaNumbering > NumberingMap;
89 typedef Alberta::DuneBoundaryProjection< dimension > Projection;
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
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 insertBoundaryProjection( GeometryTypes::simplex( dimension-1 ), vertices, 0 );
218 }
219
225 virtual void
226 insertBoundarySegment ( const std::vector< unsigned int > &vertices,
227 const std::shared_ptr< BoundarySegment > &boundarySegment )
228 {
230
231 if( !boundarySegment )
232 DUNE_THROW( GridError, "Trying to insert null as a boundary segment." );
233 if( (int)vertices.size() != refSimplex.size( dimension-1 ) )
234 DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
235
236 std::vector< WorldVector > coords( refSimplex.size( dimension-1 ) );
237 for( int i = 0; i < dimension; ++i )
238 {
239 Alberta::GlobalVector &x = macroData_.vertex( vertices[ i ] );
240 for( int j = 0; j < dimensionworld; ++j )
241 coords[ i ][ j ] = x[ j ];
242 if( ((*boundarySegment)( refSimplex.position( i, dimension-1 ) ) - coords[ i ]).two_norm() > 1e-6 )
243 DUNE_THROW( GridError, "Boundary segment does not interpolate the corners." );
244 }
245
246 const GeometryType gt = refSimplex.type( 0, 0 );
247 const DuneProjection *prj = new BoundarySegmentWrapper( gt, coords, boundarySegment );
248 insertBoundaryProjection( gt, vertices, prj );
249 }
250
264 void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
265
275 {
276 macroData_.markLongestEdge();
277 }
278
291 std::unique_ptr<Grid> createGrid ()
292 {
293 macroData_.finalize();
294 if( macroData_.elementCount() == 0 )
295 DUNE_THROW( GridError, "Cannot create empty AlbertaGrid." );
296 if( dimension < 3 )
297 macroData_.setOrientation( Alberta::Real( 1 ) );
298 assert( macroData_.checkNeighbors() );
299 macroData_.checkCycles();
300 ProjectionFactory projectionFactory( *this );
301 return std::make_unique<Grid>( macroData_, projectionFactory );
302 }
303
308 static void destroyGrid ( Grid *grid )
309 {
310 delete grid;
311 }
312
319 bool write ( const std::string &filename )
320 {
321 macroData_.finalize();
322 if( dimension < 3 )
323 macroData_.setOrientation( Alberta::Real( 1 ) );
324 assert( macroData_.checkNeighbors() );
325 return macroData_.write( filename, false );
326 }
327
328 virtual unsigned int
329 insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
330 {
331 return insertionIndex( entity.impl().elementInfo() );
332 }
333
334 virtual unsigned int
335 insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
336 {
337 const int elIndex = insertionIndex( entity.impl().elementInfo() );
338 const typename MacroData::ElementId &elementId = macroData_.element( elIndex );
339 return elementId[ entity.impl().subEntity() ];
340 }
341
342 virtual unsigned int
343 insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
344 {
345 const Grid &grid = intersection.impl().grid();
346 const ElementInfo &elementInfo = intersection.impl().elementInfo();
347 const int face = grid.generic2alberta( 1, intersection.indexInInside() );
348 return insertionIndex( elementInfo, face );
349 }
350
351 virtual bool
352 wasInserted ( const typename Grid::LeafIntersection &intersection ) const
353 {
354 return (insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
355 }
356
357 private:
358 unsigned int insertionIndex ( const ElementInfo &elementInfo ) const;
359 unsigned int insertionIndex ( const ElementInfo &elementInfo, const int face ) const;
360
361 FaceId faceId ( const ElementInfo &elementInfo, const int face ) const;
362
363 MacroData macroData_;
364 NumberingMap numberingMap_;
365 DuneProjectionPtr globalProjection_;
366 BoundaryMap boundaryMap_;
367 std::vector< DuneProjectionPtr > boundaryProjections_;
368 };
369
370
371 template< int dim, int dimworld >
372 GridFactory< AlbertaGrid< dim, dimworld > >::~GridFactory ()
373 {
374 macroData_.release();
375 }
376
377
378 template< int dim, int dimworld >
379 inline void
380 GridFactory< AlbertaGrid< dim, dimworld > >
381 ::insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift )
382 {
383 // make sure the matrix is orthogonal
384 for( int i = 0; i < dimworld; ++i )
385 for( int j = 0; j < dimworld; ++j )
386 {
387 const ctype delta = (i == j ? ctype( 1 ) : ctype( 0 ));
388 const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
389
390 if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
391 {
392 DUNE_THROW( AlbertaError,
393 "Matrix of face transformation is not orthogonal." );
394 }
395 }
396
397 // copy matrix
398 Alberta::GlobalMatrix M;
399 for( int i = 0; i < dimworld; ++i )
400 for( int j = 0; j < dimworld; ++j )
401 M[ i ][ j ] = matrix[ i ][ j ];
402
403 // copy shift
404 Alberta::GlobalVector t;
405 for( int i = 0; i < dimworld; ++i )
406 t[ i ] = shift[ i ];
407
408 // insert into ALBERTA macro data
409 macroData_.insertWallTrafo( M, t );
410 }
411
412
413 template< int dim, int dimworld >
414 inline unsigned int
416 ::insertionIndex ( const ElementInfo &elementInfo ) const
417 {
418 const MacroElement &macroElement = elementInfo.macroElement();
419 const unsigned int index = macroElement.index;
420
421#ifndef NDEBUG
422 const typename MacroData::ElementId &elementId = macroData_.element( index );
423 for( int i = 0; i <= dimension; ++i )
424 {
425 const Alberta::GlobalVector &x = macroData_.vertex( elementId[ i ] );
426 const Alberta::GlobalVector &y = macroElement.coordinate( i );
427 for( int j = 0; j < dimensionworld; ++j )
428 {
429 if( x[ j ] != y[ j ] )
430 DUNE_THROW( GridError, "Vertex in macro element does not coincide with same vertex in macro data structure." );
431 }
432 }
433#endif // #ifndef NDEBUG
434
435 return index;
436 }
437
438
439 template< int dim, int dimworld >
440 inline unsigned int
441 GridFactory< AlbertaGrid< dim, dimworld > >
442 ::insertionIndex ( const ElementInfo &elementInfo, const int face ) const
443 {
444 typedef typename BoundaryMap::const_iterator Iterator;
445 const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
446 if( it != boundaryMap_.end() )
447 return it->second;
448 else
450 }
451
452
453 template< int dim, int dimworld >
454 inline typename GridFactory< AlbertaGrid< dim, dimworld > >::FaceId
455 GridFactory< AlbertaGrid< dim, dimworld > >
456 ::faceId ( const ElementInfo &elementInfo, const int face ) const
457 {
458 const unsigned int index = insertionIndex( elementInfo );
459 const typename MacroData::ElementId &elementId = macroData_.element( index );
460
461 FaceId faceId;
462 for( size_t i = 0; i < faceId.size(); ++i )
463 {
464 const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
465 faceId[ i ] = elementId[ k ];
466 }
467 std::sort( faceId.begin(), faceId.end() );
468 return faceId;
469 }
470
471
472
473 // GridFactory::ProjectionFactory
474 // ------------------------------
475
476 template< int dim, int dimworld >
477 class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
478 : public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
479 {
480 typedef ProjectionFactory This;
481 typedef Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory > Base;
482
483 typedef typename Dune::GridFactory< AlbertaGrid< dim, dimworld > > Factory;
484
485 public:
486 typedef typename Base::Projection Projection;
487 typedef typename Base::ElementInfo ElementInfo;
488
489 typedef typename Projection::Projection DuneProjection;
490
491 ProjectionFactory( const GridFactory &gridFactory )
492 : gridFactory_( gridFactory )
493 {}
494
495 bool hasProjection ( const ElementInfo &elementInfo, const int face ) const
496 {
497 if( gridFactory().globalProjection_ )
498 return true;
499
500 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
502 return bool( gridFactory().boundaryProjections_[ index ] );
503 else
504 return false;
505 }
506
507 bool hasProjection ( const ElementInfo & /* elementInfo */ ) const
508 {
509 return bool( gridFactory().globalProjection_ );
510 }
511
512 Projection projection ( const ElementInfo &elementInfo, const int face ) const
513 {
514 const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
516 {
517 const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
518 if( projection )
519 return Projection( projection );
520 }
521
522 assert( gridFactory().globalProjection_ );
523 return Projection( gridFactory().globalProjection_ );
524 };
525
526 Projection projection ( const ElementInfo & /* elementInfo */ ) const
527 {
528 assert( gridFactory().globalProjection_ );
529 return Projection( gridFactory().globalProjection_ );
530 };
531
532 const GridFactory &gridFactory () const
533 {
534 return gridFactory_;
535 }
536
537 private:
538 const GridFactory &gridFactory_;
539 };
540
541}
542
543#endif // #if HAVE_ALBERTA
544
545#endif // #ifndef DUNE_ALBERTA_GRIDFACTORY_HH
provides the AlbertaGrid class
[ provides Dune::Grid ]
Definition: agrid.hh:107
Definition: boundaryprojection.hh:131
Wrapper class for entities.
Definition: entity.hh:64
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:369
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:328
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:18
Provide a generic factory class for unstructured grids.
Definition: gridfactory.hh:71
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition: gridfactory.hh:218
static const int dimension
dimension of the grid
Definition: gridfactory.hh:75
virtual bool wasInserted(const typename GridType::LeafIntersection &intersection) const
determine whether an intersection was inserted
Definition: gridfactory.hh:284
specialization of the generic GridFactory for AlbertaGrid
Definition: gridfactory.hh:49
FieldVector< ctype, dimensionworld > WorldVector
type of vector for world coordinates
Definition: gridfactory.hh:65
std::unique_ptr< Grid > createGrid()
finalize grid creation and hand over the grid
Definition: gridfactory.hh:291
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition: gridfactory.hh:329
AlbertaGrid< dim, dimworld > Grid
type of grid this factory is for
Definition: 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: gridfactory.hh:170
FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix
type of matrix from world coordinates to world coordinates
Definition: gridfactory.hh:67
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
insert an element into the macro grid
Definition: gridfactory.hh:125
void markLongestEdge()
mark the longest edge as refinemet edge
Definition: gridfactory.hh:274
Grid::ctype ctype
type of (scalar) coordinates
Definition: gridfactory.hh:57
virtual void insertVertex(const WorldVector &pos)
insert a vertex into the macro grid
Definition: gridfactory.hh:115
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:103
static void destroyGrid(Grid *grid)
destroy a grid previously obtain from this factory
Definition: gridfactory.hh:308
virtual unsigned int insertionIndex(const typename Codim< dimension >::Entity &entity) const
obtain a vertex' insertion index
Definition: gridfactory.hh:335
bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition: gridfactory.hh:319
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: gridfactory.hh:226
virtual void insertBoundaryProjection(const DuneProjection *projection)
insert a global (boundary) projection into the macro grid
Definition: gridfactory.hh:202
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment into the macro grid
Definition: gridfactory.hh:215
Provide a generic factory class for unstructured grids.
Definition: gridfactory.hh:312
GridFactory()
Default constructor.
Definition: gridfactory.hh:328
Grid abstract base class.
Definition: grid.hh:372
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:459
@ dimensionworld
The dimension of the world the grid lives in.
Definition: grid.hh:392
@ dimension
The dimension of the grid.
Definition: grid.hh:386
Provide a generic factory class for unstructured grids.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:156
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:461
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:11
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:92
Static tag representing a codimension.
Definition: dimension.hh:22
Interface class for vertex projection at the boundary.
Definition: boundaryprojection.hh:31
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:202
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)