dune-grid  2.3.1
alu3dgridfactory.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_ALU3DGRID_FACTORY_HH
5 #define DUNE_ALU3DGRID_FACTORY_HH
6 
7 #include <map>
8 #include <vector>
9 
10 #if HAVE_ALUGRID
11 
12 #include <dune/common/array.hh>
13 #include <dune/common/parallel/mpihelper.hh>
14 
15 #include <dune/geometry/referenceelements.hh>
16 
19 
22 
23 namespace Dune
24 {
25 
27  template< class ALUGrid >
28  class ALU3dGridFactory
29  : public GridFactoryInterface< ALUGrid >
30  {
31  typedef ALU3dGridFactory< ALUGrid > ThisType;
32  typedef GridFactoryInterface< ALUGrid > BaseType;
33 
34  public:
35  typedef ALUGrid Grid;
36 
37  typedef typename Grid::ctype ctype;
38 
39  static const ALU3dGridElementType elementType = Grid::elementType;
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 
47  typedef DuneBoundaryProjection< 3 > DuneBoundaryProjectionType;
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 
57  typedef ALUGridTransformation< ctype, dimensionworld > Transformation;
58 
60  typedef typename Transformation::WorldVector WorldVector;
62  typedef typename Transformation::WorldMatrix WorldMatrix;
63 
64  private:
65  dune_static_assert( (elementType == tetra || elementType == hexa),
66  "ALU3dGridFactory supports only grids containing "
67  "tetrahedrons or hexahedrons exclusively." );
68 
69  typedef Dune::BoundarySegmentWrapper< dimension, dimensionworld > BoundarySegmentWrapperType;
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 
75  typedef ElementTopologyMapping< elementType > ElementTopologyMappingType;
76  typedef FaceTopologyMapping< elementType > FaceTopologyMappingType;
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  {
178  insertBoundary( element, face, ALU3DSPACE ProcessorBoundary_t );
179  }
180 
189  virtual void
190  insertBoundaryProjection ( const GeometryType &type,
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 
339  template< int dimw >
340  class GridFactory< ALUSimplexGrid< 3, dimw > >
341  : public ALU3dGridFactory< ALUSimplexGrid< 3, dimw > >
342  {
343  typedef GridFactory< ALUSimplexGrid< 3, dimw > > ThisType;
344  typedef ALU3dGridFactory< ALUSimplexGrid< 3, dimw > > BaseType;
345 
346  public:
347  typedef typename BaseType::Grid Grid;
348 
349  typedef typename BaseType::MPICommunicatorType MPICommunicatorType;
350 
352  explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
353  : BaseType( communicator )
354  {}
355 
357  GridFactory ( const std::string &filename,
358  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
359  : BaseType( filename, communicator )
360  {}
361 
362  protected:
363  template< class, class, int > friend class ALULocalGeometryStorage;
365  GridFactory ( const bool realGrid,
366  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
367  : BaseType( realGrid, communicator )
368  {}
369  };
370 
371 
372 
378  template< int dimw >
379  class GridFactory< ALUCubeGrid< 3, dimw > >
380  : public ALU3dGridFactory< ALUCubeGrid< 3, dimw > >
381  {
382  typedef GridFactory< ALUCubeGrid< 3, dimw > > ThisType;
383  typedef ALU3dGridFactory< ALUCubeGrid< 3, dimw > > BaseType;
384 
385  public:
386  typedef typename BaseType::Grid Grid;
387 
388  typedef typename BaseType::MPICommunicatorType MPICommunicatorType;
389 
391  explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
392  : BaseType( communicator )
393  {}
394 
396  GridFactory ( const std::string &filename,
397  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
398  : BaseType( filename, communicator )
399  {}
400 
401  protected:
402  template< class, class, int > friend class ALULocalGeometryStorage;
404  GridFactory ( const bool realGrid,
405  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
406  : BaseType( realGrid, communicator )
407  {}
408  };
409 
410 
414  template<ALUGridElementType eltype, ALUGridRefinementType refinementtype , class Comm >
415  class GridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > >
416  : public ALU3dGridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > >
417  {
418  typedef GridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > ThisType;
419  typedef ALU3dGridFactory< ALUGrid< 3, 3, eltype, refinementtype, Comm > > BaseType;
420 
421  public:
422  typedef typename BaseType::Grid Grid;
423 
424  typedef typename BaseType::MPICommunicatorType MPICommunicatorType;
425 
427  explicit GridFactory ( const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
428  : BaseType( communicator )
429  {}
430 
432  GridFactory ( const std::string &filename,
433  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
434  : BaseType( filename, communicator )
435  {}
436 
437  protected:
438  template< class, class, int > friend class ALULocalGeometryStorage;
440  GridFactory ( const bool realGrid,
441  const MPICommunicatorType &communicator = Grid::defaultCommunicator() )
442  : BaseType( realGrid, communicator )
443  {}
444  };
445 
446 
447 
448  // Implementation of ALU3dGridFactory
449  // ----------------------------------
450 
451  template< class ALUGrid >
452  inline
453  ALU3dGridFactory< ALUGrid >
454  :: ALU3dGridFactory ( const MPICommunicatorType &communicator,
455  bool removeGeneratedFile )
456  : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
457  globalProjection_ ( 0 ),
458  numFacesInserted_ ( 0 ),
459  realGrid_( true ),
460  allowGridGeneration_( rank_ == 0 ),
461  communicator_( communicator )
462  {}
463 
464  template< class ALUGrid >
465  inline
466  ALU3dGridFactory< ALUGrid >
467  :: ALU3dGridFactory ( const std::string &filename,
468  const MPICommunicatorType &communicator )
469  : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
470  globalProjection_ ( 0 ),
471  numFacesInserted_ ( 0 ),
472  realGrid_( true ),
473  allowGridGeneration_( rank_ == 0 ),
474  communicator_( communicator )
475  {}
476 
477  template< class ALUGrid >
478  inline
479  ALU3dGridFactory< ALUGrid >
480  :: ALU3dGridFactory ( const bool realGrid,
481  const MPICommunicatorType &communicator )
482  : rank_( ALU3dGridCommunications< elementType, MPICommunicatorType >::getRank( communicator ) ),
483  globalProjection_ ( 0 ),
484  numFacesInserted_ ( 0 ),
485  realGrid_( realGrid ),
486  allowGridGeneration_( true ),
487  communicator_( communicator )
488  {}
489 
490  template< class ALUGrid >
491  inline void ALU3dGridFactory< ALUGrid > ::
492  insertBoundarySegment ( const std::vector< unsigned int >& vertices )
493  {
494  FaceType faceId;
495  copyAndSort( vertices, faceId );
496 
497  if( vertices.size() != numFaceCorners )
498  DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
499 
500  if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
501  DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
502  // DUNE_THROW( NotImplemented, "insertBoundarySegment with a single argument" );
503 
504  boundaryProjections_[ faceId ] = 0;
505 
506  BndPair boundaryId;
507  for( unsigned int i = 0; i < numFaceCorners; ++i )
508  {
509  const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i );
510  boundaryId.first[ j ] = vertices[ i ];
511  }
512  boundaryId.second = 1;
513  boundaryIds_.insert( boundaryId );
514  }
515 
516  template< class ALUGrid >
517  inline void ALU3dGridFactory< ALUGrid > ::
518  insertBoundarySegment ( const std::vector< unsigned int >& vertices,
519  const shared_ptr<BoundarySegment<3,3> >& boundarySegment )
520  {
521  FaceType faceId;
522  copyAndSort( vertices, faceId );
523 
524  if( vertices.size() != numFaceCorners )
525  DUNE_THROW( GridError, "Wrong number of face vertices passed: " << vertices.size() << "." );
526 
527  if( boundaryProjections_.find( faceId ) != boundaryProjections_.end() )
528  DUNE_THROW( GridError, "Only one boundary projection can be attached to a face." );
529 
530  const size_t numVx = vertices.size();
531  GeometryType type;
532  if( numVx == 3 )
533  type.makeSimplex( dimension-1 );
534  else
535  type.makeCube( dimension-1 );
536 
537  // we need double here because of the structure of BoundarySegment
538  // and BoundarySegmentWrapper which have double as coordinate type
539  typedef FieldVector< double, dimensionworld > CoordType;
540  std::vector< CoordType > coords( numVx );
541  for( size_t i = 0; i < numVx; ++i )
542  {
543  // if this assertions is thrown vertices were not inserted at first
544  assert( vertices_.size() > vertices[ i ] );
545 
546  // get global coordinate and copy it
547  const VertexType &x = position( vertices[ i ] );
548  for( unsigned int j = 0; j < dimensionworld; ++j )
549  coords[ i ][ j ] = x[ j ];
550  }
551 
552  BoundarySegmentWrapperType* prj
553  = new BoundarySegmentWrapperType( type, coords, boundarySegment );
554  boundaryProjections_[ faceId ] = prj;
555 #ifndef NDEBUG
556  // consistency check
557  for( size_t i = 0; i < numVx; ++i )
558  {
559  CoordType global = (*prj)( coords [ i ] );
560  if( (global - coords[ i ]).two_norm() > 1e-6 )
561  DUNE_THROW(GridError,"BoundarySegment does not map face vertices to face vertices.");
562  }
563 #endif
564 
565  BndPair boundaryId;
566  for( unsigned int i = 0; i < numFaceCorners; ++i )
567  {
568  const unsigned int j = FaceTopologyMappingType::dune2aluVertex( i );
569  boundaryId.first[ j ] = vertices[ i ];
570  }
571  boundaryId.second = 1;
572  boundaryIds_.insert( boundaryId );
573  }
574 
575 
576  template< class ALUGrid >
577  inline void ALU3dGridFactory< ALUGrid >
578  ::generateFace ( const SubEntity &subEntity, FaceType &face ) const
579  {
580  generateFace( elements_[ subEntity.first ], subEntity.second, face );
581  }
582 
583 } // end namespace Dune
584 
585 #endif // #if HAVE_ALUGRID
586 
587 #if COMPILE_ALUGRID_INLINE
588  #include "alu3dgridfactory.cc"
589 #endif
590 #endif