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