dune-grid  2.3.1
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 <limits>
14 #include <map>
15 
16 #include <dune/common/array.hh>
17 
18 #include <dune/geometry/referenceelements.hh>
19 
21 
23 
25 
26 #if HAVE_ALBERTA
27 
28 namespace 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 
67  typedef FieldVector< ctype, dimensionworld > WorldVector;
69  typedef FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix;
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
86 
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 
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 
167  virtual void
169  const std::vector< unsigned int > &vertices,
170  const DuneProjection *projection )
171  {
172  if( (int)type.dim() != dimension-1 )
173  DUNE_THROW( AlbertaError, "Inserting boundary face of wrong dimension: " << type.dim() );
174  if( !type.isSimplex() )
175  DUNE_THROW( AlbertaError, "Alberta supports only simplices." );
176 
177  FaceId faceId;
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() );
183 
184  typedef std::pair< typename BoundaryMap::iterator, bool > InsertResult;
185  const InsertResult result = boundaryMap_.insert( std::make_pair( faceId, boundaryProjections_.size() ) );
186  if( !result.second )
187  DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
188  boundaryProjections_.push_back( DuneProjectionPtr( projection ) );
189  }
190 
191 
198  virtual void insertBoundaryProjection ( const DuneProjection *projection )
199  {
200  if( globalProjection_ )
201  DUNE_THROW( GridError, "Only one global boundary projection can be attached to a grid." );
202  globalProjection_ = DuneProjectionPtr( projection );
203  }
204 
210  virtual void
211  insertBoundarySegment ( const std::vector< unsigned int >& vertices )
212  {
213  typedef typename GenericGeometry::SimplexTopology< dimension-1 >::type Topology;
214  insertBoundaryProjection( GeometryType( Topology() ), vertices, 0 );
215  }
216 
222  virtual void
223  insertBoundarySegment ( const std::vector< unsigned int > &vertices,
224  const shared_ptr< BoundarySegment > &boundarySegment )
225  {
226  const ReferenceElement< ctype, dimension-1 > &refSimplex
228 
229  if( !boundarySegment )
230  DUNE_THROW( GridError, "Trying to insert null as a boundary segment." );
231  if( (int)vertices.size() != refSimplex.size( dimension-1 ) )
232  DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
233 
234  std::vector< WorldVector > coords( refSimplex.size( dimension-1 ) );
235  for( int i = 0; i < dimension; ++i )
236  {
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 )
241  DUNE_THROW( GridError, "Boundary segment does not interpolate the corners." );
242  }
243 
244  const GeometryType gt = refSimplex.type( 0, 0 );
245  const DuneProjection *prj = new BoundarySegmentWrapper( gt, coords, boundarySegment );
246  insertBoundaryProjection( gt, vertices, prj );
247  }
248 
262  void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
263 
273  {
274  macroData_.markLongestEdge();
275  }
276 
290  {
291  macroData_.finalize();
292  if( macroData_.elementCount() == 0 )
293  DUNE_THROW( GridError, "Cannot create empty AlbertaGrid." );
294  if( dimension < 3 )
295  macroData_.setOrientation( Alberta::Real( 1 ) );
296  assert( macroData_.checkNeighbors() );
297  macroData_.checkCycles();
298  ProjectionFactory projectionFactory( *this );
299  return new Grid( macroData_, projectionFactory );
300  }
301 
306  static void destroyGrid ( Grid *grid )
307  {
308  delete grid;
309  }
310 
319  template< GrapeIOFileFormatType type >
320  bool write ( const std::string &filename )
321  {
322  dune_static_assert( type != pgm, "AlbertaGridFactory: writing pgm format is not supported." );
323  macroData_.finalize();
324  if( dimension < 3 )
325  macroData_.setOrientation( Alberta::Real( 1 ) );
326  assert( macroData_.checkNeighbors() );
327  return macroData_.write( filename, (type == xdr) );
328  }
329 
338  virtual bool write ( const std::string &filename )
339  {
340  return write< ascii >( filename );
341  }
342 
343  virtual unsigned int
344  insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
345  {
346  return insertionIndex( Grid::getRealImplementation( entity ).elementInfo() );
347  }
348 
349  virtual unsigned int
350  insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
351  {
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() ];
355  }
356 
357  virtual unsigned int
358  insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
359  {
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() );
363  return insertionIndex( elementInfo, face );
364  }
365 
366  virtual bool
367  wasInserted ( const typename Grid::LeafIntersection &intersection ) const
368  {
369  return (insertionIndex( intersection ) < std::numeric_limits< unsigned int >::max());
370  }
371 
372  private:
373  unsigned int insertionIndex ( const ElementInfo &elementInfo ) const;
374  unsigned int insertionIndex ( const ElementInfo &elementInfo, const int face ) const;
375 
376  FaceId faceId ( const ElementInfo &elementInfo, const int face ) const;
377 
378  MacroData macroData_;
379  NumberingMap numberingMap_;
380  DuneProjectionPtr globalProjection_;
381  BoundaryMap boundaryMap_;
382  std::vector< DuneProjectionPtr > boundaryProjections_;
383  };
384 
385 
386  template< int dim, int dimworld >
388  {
389  macroData_.release();
390  }
391 
392 
393  template< int dim, int dimworld >
394  inline void
396  ::insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift )
397  {
398  // make sure the matrix is orthogonal
399  for( int i = 0; i < dimworld; ++i )
400  for( int j = 0; j < dimworld; ++j )
401  {
402  const ctype delta = (i == j ? ctype( 1 ) : ctype( 0 ));
403  const ctype epsilon = (8*dimworld)*std::numeric_limits< ctype >::epsilon();
404 
405  if( std::abs( matrix[ i ] * matrix[ j ] - delta ) > epsilon )
406  {
407  DUNE_THROW( AlbertaError,
408  "Matrix of face transformation is not orthogonal." );
409  }
410  }
411 
412  // copy matrix
414  for( int i = 0; i < dimworld; ++i )
415  for( int j = 0; j < dimworld; ++j )
416  M[ i ][ j ] = matrix[ i ][ j ];
417 
418  // copy shift
420  for( int i = 0; i < dimworld; ++i )
421  t[ i ] = shift[ i ];
422 
423  // insert into ALBERTA macro data
424  macroData_.insertWallTrafo( M, t );
425  }
426 
427 
428  template< int dim, int dimworld >
429  inline unsigned int
431  ::insertionIndex ( const ElementInfo &elementInfo ) const
432  {
433  const MacroElement &macroElement = elementInfo.macroElement();
434  const unsigned int index = macroElement.index;
435 
436 #ifndef NDEBUG
437  const typename MacroData::ElementId &elementId = macroData_.element( index );
438  for( int i = 0; i <= dimension; ++i )
439  {
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 )
443  {
444  if( x[ j ] != y[ j ] )
445  DUNE_THROW( GridError, "Vertex in macro element does not coincide with same vertex in macro data structure." );
446  }
447  }
448 #endif // #ifndef NDEBUG
449 
450  return index;
451  }
452 
453 
454  template< int dim, int dimworld >
455  inline unsigned int
456  GridFactory< AlbertaGrid< dim, dimworld > >
457  ::insertionIndex ( const ElementInfo &elementInfo, const int face ) const
458  {
459  typedef typename BoundaryMap::const_iterator Iterator;
460  const Iterator it = boundaryMap_.find( faceId( elementInfo, face ) );
461  if( it != boundaryMap_.end() )
462  return it->second;
463  else
465  }
466 
467 
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
472  {
473  const unsigned int index = insertionIndex( elementInfo );
474  const typename MacroData::ElementId &elementId = macroData_.element( index );
475 
476  FaceId faceId;
477  for( size_t i = 0; i < faceId.size(); ++i )
478  {
479  const int k = Alberta::MapVertices< dimension, 1 >::apply( face, i );
480  faceId[ i ] = elementId[ k ];
481  }
482  std::sort( faceId.begin(), faceId.end() );
483  return faceId;
484  }
485 
486 
487 
488  // GridFactory::ProjectionFactory
489  // ------------------------------
490 
491  template< int dim, int dimworld >
492  class GridFactory< AlbertaGrid< dim, dimworld > >::ProjectionFactory
493  : public Alberta::ProjectionFactory< Alberta::DuneBoundaryProjection< dim >, ProjectionFactory >
494  {
495  typedef ProjectionFactory This;
497 
499 
500  public:
501  typedef typename Base::Projection Projection;
502  typedef typename Base::ElementInfo ElementInfo;
503 
505 
506  ProjectionFactory( const GridFactory &gridFactory )
507  : gridFactory_( gridFactory )
508  {}
509 
510  bool hasProjection ( const ElementInfo &elementInfo, const int face ) const
511  {
512  if( gridFactory().globalProjection_ )
513  return true;
514 
515  const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
517  return bool( gridFactory().boundaryProjections_[ index ] );
518  else
519  return false;
520  }
521 
522  bool hasProjection ( const ElementInfo &elementInfo ) const
523  {
524  return bool( gridFactory().globalProjection_ );
525  }
526 
527  Projection projection ( const ElementInfo &elementInfo, const int face ) const
528  {
529  const unsigned int index = gridFactory().insertionIndex( elementInfo, face );
531  {
532  const DuneProjectionPtr &projection = gridFactory().boundaryProjections_[ index ];
533  if( projection )
534  return Projection( projection );
535  }
536 
537  assert( gridFactory().globalProjection_ );
538  return Projection( gridFactory().globalProjection_ );
539  };
540 
541  Projection projection ( const ElementInfo &elementInfo ) const
542  {
543  assert( gridFactory().globalProjection_ );
544  return Projection( gridFactory().globalProjection_ );
545  };
546 
547  const GridFactory &gridFactory () const
548  {
549  return gridFactory_;
550  }
551 
552  private:
553  const GridFactory &gridFactory_;
554  };
555 
556 }
557 
558 #endif // #if HAVE_ALBERTA
559 
560 #endif // #ifndef DUNE_ALBERTA_GRIDFACTORY_HH
Dune::shared_ptr< const DuneProjection > DuneProjectionPtr
Definition: albertagrid/gridfactory.hh:72
Base::ElementInfo ElementInfo
Definition: albertagrid/gridfactory.hh:502
Definition: grapedataioformattypes.hh:16
int max(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:341
GridFamily::ctype ctype
Definition: agrid.hh:175
Definition: grapedataioformattypes.hh:18
ProjectionFactory(const GridFactory &gridFactory)
Definition: albertagrid/gridfactory.hh:506
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:178
Definition: albertagrid/projection.hh:25
Projection projection(const ElementInfo &elementInfo) const
Definition: albertagrid/gridfactory.hh:541
Base class for classes implementing geometries of boundary segments.
Definition: boundarysegment.hh:29
virtual void insertBoundaryProjection(const DuneProjection *projection)
insert a global (boundary) projection into the macro grid
Definition: albertagrid/gridfactory.hh:198
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition: common/gridfactory.hh:179
FieldVector< ctype, dimensionworld > WorldVector
type of vector for world coordinates
Definition: albertagrid/gridfactory.hh:67
Provide a generic factory class for unstructured grids.
Definition: common/gridfactory.hh:263
Definition: albertagrid/projection.hh:33
void markLongestEdge()
mark the longest edge as refinemet edge
Definition: albertagrid/gridfactory.hh:272
Provide a generic factory class for unstructured grids.
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
insert an element into the macro grid
Definition: albertagrid/gridfactory.hh:127
Base::Projection Projection
Definition: albertagrid/gridfactory.hh:501
virtual void insertBoundarySegment(const std::vector< unsigned int > &vertices)
insert a boundary segment into the macro grid
Definition: albertagrid/gridfactory.hh:211
virtual unsigned int insertionIndex(const typename Codim< 0 >::Entity &entity) const
obtain an element's insertion index
Definition: albertagrid/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: albertagrid/gridfactory.hh:223
Definition: macroelement.hh:20
static const int dimension
dimension of the grid
Definition: common/gridfactory.hh:78
const GridFactory & gridFactory() const
Definition: albertagrid/gridfactory.hh:547
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:484
FieldMatrix< ctype, dimensionworld, dimensionworld > WorldMatrix
type of matrix from world coordinates to world coordinates
Definition: albertagrid/gridfactory.hh:69
Projection projection(const ElementInfo &elementInfo, const int face) const
Definition: albertagrid/gridfactory.hh:527
ALBERTA REAL_D GlobalVector
Definition: misc.hh:47
int generic2alberta(int codim, int i) const
Definition: agrid.hh:542
The dimension of the world the grid lives in.
Definition: common/grid.hh:406
bool hasProjection(const ElementInfo &elementInfo) const
Definition: albertagrid/gridfactory.hh:522
bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition: albertagrid/gridfactory.hh:320
AlbertaGrid< dim, dimworld > Grid
type of grid this factory is for
Definition: albertagrid/gridfactory.hh:56
Grid * createGrid()
finalize grid creation and hand over the grid
Definition: albertagrid/gridfactory.hh:289
Provide a generic factory class for unstructured grids.
Definition: common/gridfactory.hh:73
The dimension of the grid.
Definition: common/grid.hh:400
virtual bool write(const std::string &filename)
write out the macro triangulation in native grid file format
Definition: albertagrid/gridfactory.hh:338
[ provides Dune::Grid ]
Definition: agrid.hh:137
virtual bool wasInserted(const typename Grid::LeafIntersection &intersection) const
Definition: albertagrid/gridfactory.hh:367
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:16
virtual void insertBoundary(int element, int face, int id)
mark a face as boundary (and assign a boundary id)
Definition: albertagrid/gridfactory.hh:152
Interface class for vertex projection at the boundary.
Definition: boundaryprojection.hh:23
Grid::template Codim< codim >::Entity Entity
Definition: albertagrid/gridfactory.hh:78
Definition: misc.hh:159
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:168
Definition: misc.hh:27
bool hasProjection(const ElementInfo &elementInfo, const int face) const
Definition: albertagrid/gridfactory.hh:510
provides the AlbertaGrid class
DuneBoundaryProjection< dimensionworld > DuneProjection
Definition: albertagrid/gridfactory.hh:71
Dune::BoundarySegment< dimension, dimensionworld > BoundarySegment
Definition: albertagrid/gridfactory.hh:73
Definition: alugrid/common/declaration.hh:18
Grid::ctype ctype
type of (scalar) coordinates
Definition: albertagrid/gridfactory.hh:59
ALBERTA REAL_DD GlobalMatrix
Definition: misc.hh:48
ALBERTA REAL Real
Definition: misc.hh:45
static void destroyGrid(Grid *grid)
destroy a grid previously obtain from this factory
Definition: albertagrid/gridfactory.hh:306
Projection::Projection DuneProjection
Definition: albertagrid/gridfactory.hh:504
GridFactory()
Definition: albertagrid/gridfactory.hh:105
virtual unsigned int insertionIndex(const typename Grid::LeafIntersection &intersection) const
Definition: albertagrid/gridfactory.hh:358
Definition: boundaryprojection.hh:65
virtual unsigned int insertionIndex(const typename Codim< dimension >::Entity &entity) const
obtain a vertex' insertion index
Definition: albertagrid/gridfactory.hh:350
void abs(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:332
specialization of the generic GridFactory for AlbertaGrid
Definition: albertagrid/gridfactory.hh:49
virtual void insertVertex(const WorldVector &pos)
insert a vertex into the macro grid
Definition: albertagrid/gridfactory.hh:117