dune-grid  2.2.0
alu3dgridfactory.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef DUNE_ALU3DGRID_FACTORY_HH
5 #define DUNE_ALU3DGRID_FACTORY_HH
6 
7 #include <map>
8 #include <vector>
9 
10 #ifdef ENABLE_ALUGRID
11 
12 #include <dune/common/array.hh>
13 #include <dune/common/mpihelper.hh>
14 
15 #include <dune/geometry/referenceelements.hh>
16 
19 
22 
23 namespace Dune
24 {
25 
27  template< class ALUGrid >
29  : public GridFactoryInterface< ALUGrid >
30  {
33 
34  public:
35  typedef ALUGrid Grid;
36 
37  typedef typename Grid::ctype ctype;
38 
40 
41  static const unsigned int dimension = Grid::dimension;
42  static const unsigned int dimensionworld = Grid::dimensionworld;
43 
44  typedef typename Grid::MPICommunicatorType MPICommunicatorType;
45 
48 
49  template< int codim >
50  struct Codim
51  {
52  typedef typename Grid::template Codim< codim >::Entity Entity;
53  };
54 
55  typedef unsigned int VertexId;
56 
58 
63 
64  private:
65  dune_static_assert( (elementType == tetra || elementType == hexa),
66  "ALU3dGridFactory supports only grids containing "
67  "tetrahedrons or hexahedrons exclusively." );
68 
70 
71  static const unsigned int numCorners = EntityCount< elementType >::numVertices;
72  static const unsigned int numFaces = EntityCount< elementType >::numFaces;
73  static const unsigned int numFaceCorners = EntityCount< elementType >::numVerticesPerFace;
74 
77 
78  typedef FieldVector< ctype, dimensionworld > VertexType;
79  typedef std::vector< unsigned int > ElementType;
80  typedef array< unsigned int, numFaceCorners > FaceType;
81 
82  struct FaceLess;
83 
84  typedef std::vector< std::pair< VertexType, size_t > > VertexVector;
85  typedef std::vector< ElementType > ElementVector;
86  typedef std::pair< FaceType, int > BndPair ;
87  typedef std::map< FaceType, int > BoundaryIdMap;
88  typedef std::vector< std::pair< BndPair, BndPair > > PeriodicBoundaryVector;
89  typedef std::pair< unsigned int, int > SubEntity;
90  typedef std::map< FaceType, SubEntity, FaceLess > FaceMap;
91 
92  typedef std::map< FaceType, const DuneBoundaryProjectionType* > BoundaryProjectionMap;
93  typedef std::vector< const DuneBoundaryProjectionType* > BoundaryProjectionVector;
94 
95  typedef std::vector< Transformation > FaceTransformationVector;
96 
97  // copy vertex numbers and store smalled #dimension ones
98  void copyAndSort ( const std::vector< unsigned int > &vertices, FaceType &faceId ) const
99  {
100  std::vector<unsigned int> tmp( vertices );
101  std::sort( tmp.begin(), tmp.end() );
102 
103  // copy only the first dimension vertices (enough for key)
104  for( size_t i = 0; i < faceId.size(); ++i ) faceId[ i ] = tmp[ i ];
105  }
106 
107  private:
108  // return grid object
109  virtual Grid* createGridObj( BoundaryProjectionVector* bndProjections, const std::string& name ) const
110  {
111  return ( allowGridGeneration_ ) ?
112  new Grid( communicator_, globalProjection_, bndProjections , name, realGrid_ ) :
113  new Grid( communicator_ );
114  }
115 
116  public:
118  explicit ALU3dGridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator(),
119  bool removeGeneratedFile = true );
120 
122  explicit ALU3dGridFactory ( const std::string &filename,
123  const MPICommunicatorType &communicator = Grid::defaultCommunicator() );
124 
126  explicit ALU3dGridFactory ( const bool verbose, const MPICommunicatorType &communicator );
127 
129  virtual ~ALU3dGridFactory ();
130 
135  virtual void insertVertex ( const VertexType &pos );
136 
137  // for testing parallel GridFactory
138  VertexId insertVertex ( const VertexType &pos, const size_t globalId );
139 
148  virtual void
149  insertElement ( const GeometryType &geometry,
150  const std::vector< VertexId > &vertices );
151 
162  virtual void
163  insertBoundary ( const GeometryType &geometry,
164  const std::vector< VertexId > &faceVertices,
165  const int id );
166 
173  virtual void insertBoundary ( const int element, const int face, const int id );
174 
175  // for testing parallel GridFactory
176  void insertProcessBorder ( const int element, const int face )
177  {
179  }
180 
189  virtual void
191  const std::vector< VertexId > &vertices,
192  const DuneBoundaryProjectionType *projection );
193 
198  virtual void
199  insertBoundarySegment ( const std::vector< VertexId >& vertices ) ;
200 
206  virtual void
207  insertBoundarySegment ( const std::vector< VertexId >& vertices,
208  const shared_ptr<BoundarySegment<3,3> >& boundarySegment ) ;
209 
214  virtual void insertBoundaryProjection ( const DuneBoundaryProjectionType& bndProjection );
215 
225  void insertFaceTransformation ( const WorldMatrix &matrix, const WorldVector &shift );
226 
231  Grid *createGrid ();
232 
233  Grid *createGrid ( const bool addMissingBoundaries, const std::string dgfName = "" );
234 
235  Grid *createGrid ( const bool addMissingBoundaries, bool temporary, const std::string dgfName = "" );
236 
237  virtual unsigned int
238  insertionIndex ( const typename Codim< 0 >::Entity &entity ) const
239  {
240  return Grid::getRealImplementation( entity ).getIndex();
241  }
242  virtual unsigned int
243  insertionIndex ( const typename Codim< dimension >::Entity &entity ) const
244  {
245  return Grid::getRealImplementation( entity ).getIndex();
246  }
247  virtual unsigned int
248  insertionIndex ( const typename Grid::LeafIntersection &intersection ) const
249  {
250  return intersection.boundarySegmentIndex();
251  }
252  virtual bool
253  wasInserted ( const typename Grid::LeafIntersection &intersection ) const
254  {
255  return intersection.boundary() &&
256  ( insertionIndex(intersection) < numFacesInserted_ );
257  }
258 
259  private:
260  size_t globalId ( const VertexId &id ) const
261  {
262  assert( id < vertices_.size() );
263  return vertices_[ id ].second;
264  }
265 
266  const VertexType &position ( const VertexId &id ) const
267  {
268  assert( id < vertices_.size() );
269  return vertices_[ id ].first;
270  }
271 
272  void assertGeometryType( const GeometryType &geometry );
273  static void generateFace ( const ElementType &element, const int f, FaceType &face );
274  void generateFace ( const SubEntity &subEntity, FaceType &face ) const;
275  void correctElementOrientation ();
276  bool identifyFaces ( const Transformation &transformation, const FaceType &key1, const FaceType &key2, const int defaultId );
277  void searchPeriodicNeighbor ( FaceMap &faceMap, const typename FaceMap::iterator &pos, const int defaultId );
278  void reinsertBoundary ( const FaceMap &faceMap, const typename FaceMap::const_iterator &pos, const int id );
279  void recreateBoundaryIds ( const int defaultId = 1 );
280 
281  int rank_;
282 
283  VertexVector vertices_;
284  ElementVector elements_;
285  BoundaryIdMap boundaryIds_;
286  PeriodicBoundaryVector periodicBoundaries_;
287  const DuneBoundaryProjectionType* globalProjection_ ;
288  BoundaryProjectionMap boundaryProjections_;
289  FaceTransformationVector faceTransformations_;
290  unsigned int numFacesInserted_;
291  bool realGrid_;
292  const bool allowGridGeneration_;
293 
294  MPICommunicatorType communicator_;
295  };
296 
297 
298 
299  template< class ALUGrid >
300  struct ALU3dGridFactory< ALUGrid >::FaceLess
301  : public std::binary_function< FaceType, FaceType, bool >
302  {
303  bool operator() ( const FaceType &a, const FaceType &b ) const
304  {
305  for( unsigned int i = 0; i < numFaceCorners; ++i )
306  {
307  if( a[ i ] != b[ i ] )
308  return (a[ i ] < b[ i ]);
309  }
310  return false;
311  }
312  };
313 
314 
315  template< class ALUGrid >
316  inline void ALU3dGridFactory< ALUGrid >
317  ::assertGeometryType( const GeometryType &geometry )
318  {
319  if( elementType == tetra )
320  {
321  if( !geometry.isSimplex() )
322  DUNE_THROW( GridError, "Only simplex geometries can be inserted into "
323  "ALUSimplexGrid< 3, 3 >." );
324  }
325  else
326  {
327  if( !geometry.isCube() )
328  DUNE_THROW( GridError, "Only cube geometries can be inserted into "
329  "ALUCubeGrid< 3, 3 >." );
330  }
331  }
332 
333 
337  template<>
338  class GridFactory< ALUSimplexGrid< 3, 3 > >
339  : public ALU3dGridFactory< ALUSimplexGrid< 3, 3 > >
340  {
343 
344  public:
346 
348 
350  explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
351  : BaseType( communicator )
352  {}
353 
355  GridFactory ( const std::string &filename,
356  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
357  : BaseType( filename, communicator )
358  {}
359 
360  protected:
361  template< class, class, int > friend class ALULocalGeometryStorage;
363  GridFactory ( const bool realGrid,
364  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
365  : BaseType( realGrid, communicator )
366  {}
367  };
368 
369 
370 
374  template<>
375  class GridFactory< ALUCubeGrid< 3, 3 > >
376  : public ALU3dGridFactory< ALUCubeGrid< 3, 3 > >
377  {
380 
381  public:
383 
385 
387  explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
388  : BaseType( communicator )
389  {}
390 
392  GridFactory ( const std::string &filename,
393  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
394  : BaseType( filename, communicator )
395  {}
396 
397  protected:
398  template< class, class, int > friend class ALULocalGeometryStorage;
400  GridFactory ( const bool realGrid,
401  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
402  : BaseType( realGrid, communicator )
403  {}
404  };
405 
406 
410  template<ALUGridElementType eltype, ALUGridRefinementType refinementtype , class Comm >
411  class GridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > >
412  : public ALU3dGridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > >
413  {
416 
417  public:
418  typedef typename BaseType::Grid Grid;
419 
421 
423  explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
424  : BaseType( communicator )
425  {}
426 
428  GridFactory ( const std::string &filename,
429  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
430  : BaseType( filename, communicator )
431  {}
432 
433  protected:
434  template< class, class, int > friend class ALULocalGeometryStorage;
436  GridFactory ( const bool realGrid,
437  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
438  : BaseType( realGrid, communicator )
439  {}
440  };
441 
442 
443 
444  // Implementation of ALU3dGridFactory
445  // ----------------------------------
446 
447  template< class ALUGrid >
448  inline
451  bool removeGeneratedFile )
452  : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
453  globalProjection_ ( 0 ),
454  numFacesInserted_ ( 0 ),
455  realGrid_( true ),
456  allowGridGeneration_( rank_ == 0 ),
457  communicator_( communicator )
458  {}
459 
460  template< class ALUGrid >
461  inline
463  :: ALU3dGridFactory ( const std::string &filename,
464  const MPICommunicatorType &communicator )
465  : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
466  globalProjection_ ( 0 ),
467  numFacesInserted_ ( 0 ),
468  realGrid_( true ),
469  allowGridGeneration_( rank_ == 0 ),
470  communicator_( communicator )
471  {}
472 
473  template< class ALUGrid >
474  inline
476  :: ALU3dGridFactory ( const bool realGrid,
477  const MPICommunicatorType &communicator )
478  : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
479  globalProjection_ ( 0 ),
480  numFacesInserted_ ( 0 ),
481  realGrid_( realGrid ),
482  allowGridGeneration_( true ),
483  communicator_( communicator )
484  {}
485 
486  template< class ALUGrid >
488  insertBoundarySegment ( const std::vector< unsigned int >& vertices )
489  {
490  FaceType faceId;
491  copyAndSort( vertices, faceId );
492 
493  if( vertices.size() != numFaceCorners )
494  DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
495 
496  if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
497  DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
498  // DUNE_THROW( NotImplemented, "insertBoundarySegment with a single argument" );
499 
500  boundaryProjections_[ faceId ] = 0;
501 
502  BndPair boundaryId;
503  for( unsigned int i = 0; i < numFaceCorners; ++i )
504  {
505  const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i );
506  boundaryId.first[ j ] = vertices[ i ];
507  }
508  boundaryId.second = 1;
509  boundaryIds_.insert( boundaryId );
510  }
511 
512  template< class ALUGrid >
514  insertBoundarySegment ( const std::vector< unsigned int >& vertices,
515  const shared_ptr<BoundarySegment<3,3> >& boundarySegment )
516  {
517  FaceType faceId;
518  copyAndSort( vertices, faceId );
519 
520  if( vertices.size() != numFaceCorners )
521  DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
522 
523  if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
524  DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
525 
526  const size_t numVx = vertices.size();
527  GeometryType type;
528  if( numVx == 3 )
529  type.makeSimplex( dimension-1 );
530  else
531  type.makeCube( dimension-1 );
532 
533  // we need double here because of the structure of BoundarySegment
534  // and BoundarySegmentWrapper which have double as coordinate type
535  typedef FieldVector< double, dimensionworld > CoordType;
536  std::vector< CoordType > coords( numVx );
537  for( size_t i = 0; i < numVx; ++i )
538  {
539  // if this assertions is thrown vertices were not inserted at first
540  assert( vertices_.size() > vertices[ i ] );
541 
542  // get global coordinate and copy it
543  const VertexType &x = position( vertices[ i ] );
544  for( unsigned int j = 0; j < dimensionworld; ++j )
545  coords[ i ][ j ] = x[ j ];
546  }
547 
549  = new BoundarySegmentWrapperType( type, coords, boundarySegment );
550  boundaryProjections_[ faceId ] = prj;
551 #ifndef NDEBUG
552  // consistency check
553  for( size_t i = 0; i < numVx; ++i )
554  {
555  CoordType global = (*prj)( coords [ i ] );
556  if( (global - coords[ i ]).two_norm() > 1e-6 )
557  DUNE_THROW(GridError,"BoundarySegment does not map face vertices to face vertices.");
558  }
559 #endif
560 
561  BndPair boundaryId;
562  for( unsigned int i = 0; i < numFaceCorners; ++i )
563  {
564  const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i );
565  boundaryId.first[ j ] = vertices[ i ];
566  }
567  boundaryId.second = 1;
568  boundaryIds_.insert( boundaryId );
569  }
570 
571 
572  template< class ALUGrid >
573  inline void ALU3dGridFactory< ALUGrid >
574  ::generateFace ( const SubEntity &subEntity, FaceType &face ) const
575  {
576  generateFace( elements_[ subEntity.first ], subEntity.second, face );
577  }
578 
579 } // end namespace Dune
580 
581 #endif // #ifdef ENABLE_ALUGRID
582 
583 #if COMPILE_ALUGRID_INLINE
584  #include "alu3dgridfactory.cc"
585 #endif
586 #endif