dune-alugrid  2.6-git
dgf.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALUGRID_DGF_HH
2 #define DUNE_ALUGRID_DGF_HH
3 
4 #include <type_traits>
5 
6 #if HAVE_ALUGRID
7 #include <dune/alugrid/grid.hh>
8 #include <dune/grid/io/file/dgfparser/dgfalu.hh>
9 #else
10 
11 // include grid first to avoid includes from dune-grid/grid/alugrid
12 #include <dune/alugrid/grid.hh>
13 
14 #include <dune/grid/common/intersection.hh>
15 #include <dune/grid/io/file/dgfparser/dgfparser.hh>
16 #include <dune/grid/io/file/dgfparser/parser.hh>
17 #include <dune/grid/io/file/dgfparser/blocks/projection.hh>
18 
19 
20 namespace Dune
21 {
22 
23  namespace
24  {
25 
26  // GlobalVertexIndexBlock
27  // ----------------------
28 
29  class GlobalVertexIndexBlock
30  : public dgf::BasicBlock
31  {
32  bool goodline;
33 
34  public:
35  GlobalVertexIndexBlock ( std :: istream &in )
36  : dgf::BasicBlock( in, "GlobalVertexIndex" ),
37  goodline( true )
38  {}
39 
40  bool next ( int &index )
41  {
42  assert( ok() );
43  if( !getnextline() )
44  return (goodline = false);
45 
46  if( !getnextentry( index ) )
47  {
48  DUNE_THROW ( DGFException, "Error in " << *this << ": "
49  << "Wrong global vertex indices " );
50  }
51  return (goodline = true);
52  }
53 
54  // some information
55  bool ok ()
56  {
57  return goodline;
58  }
59  };
60 
61 
62 
63  // ALUParallelBlock
64  // ----------------
65 
66  class ALUParallelBlock
67  : public dgf::BasicBlock
68  {
69  bool goodline;
70 
71  public:
72  ALUParallelBlock ( std :: istream &in )
73  : dgf::BasicBlock( in, "ALUParallel" ),
74  goodline( true )
75  {}
76 
77  bool next ( std::string &name )
78  {
79  assert( ok() );
80  if( !getnextline() )
81  return (goodline = false);
82 
83  if( !getnextentry( name ) )
84  {
85  DUNE_THROW ( DGFException, "Error in " << *this << ": "
86  << "Wrong global vertex indices " );
87  }
88  return (goodline = true);
89  }
90 
91  // some information
92  bool ok ()
93  {
94  return goodline;
95  }
96  };
97 
98  } // end empty namespace
99 
100 
101 
102  // DGFGridInfo (specialization for ALUGrid)
103  // ----------------------------------------
104 
105  template<int dimg, int dimw, ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm >
106  struct DGFGridInfo< Dune::ALUGrid< dimg, dimw, eltype, refinementtype, Comm > >
107  {
108  static int refineStepsForHalf () { return ( refinementtype == conforming ) ? dimg : 1; }
109  static double refineWeight () { return ( refinementtype == conforming ) ? 0.5 : 1.0/(std::pow( 2.0, double(dimg))); }
110  };
115  // DGFGridFactory for AluGrid
116  // --------------------------
117 
118  // template< int dim, int dimworld > // for a first version
119  template< class G >
121  {
122  typedef G Grid;
123  const static int dimension = Grid::dimension;
124  typedef MPIHelper::MPICommunicator MPICommunicatorType;
125  typedef typename Grid::template Codim<0>::Entity Element;
126  typedef typename Grid::template Codim<dimension>::Entity Vertex;
127  typedef Dune::GridFactory<Grid> GridFactory;
128 
130  : factory_( ),
131  dgf_( 0, 1 )
132  {}
133 
135  : factory_(comm),
136  dgf_( rank(comm), size(comm) )
137  {}
138 
139  Grid *grid () const
140  {
141  return grid_;
142  }
143 
144  template< class Intersection >
145  bool wasInserted ( const Intersection &intersection ) const
146  {
147  return factory_.wasInserted( intersection );
148  }
149 
150  template< class GG, class II >
151  int boundaryId ( const Intersection< GG, II > & intersection ) const
152  {
153  typedef Dune::Intersection< GG, II > Intersection;
154  const typename Intersection::Entity & entity = intersection.inside();
155 
156  const int face = intersection.indexInInside();
157 
158  const auto& refElem =
159  ReferenceElements< double, dimension >::general( entity.type() );
160  int corners = refElem.size( face, 1, dimension );
161  std :: vector< unsigned int > bound( corners );
162  for( int i=0; i < corners; ++i )
163  {
164  const int k = refElem.subEntity( face, 1, i, dimension );
165  bound[ i ] = factory_.insertionIndex( entity.template subEntity< dimension >( k ) );
166  }
167 
168  DuneGridFormatParser::facemap_t::key_type key( bound, false );
169  const DuneGridFormatParser::facemap_t::const_iterator pos = dgf_.facemap.find( key );
170  if( pos != dgf_.facemap.end() )
171  return dgf_.facemap.find( key )->second.first;
172  else
173  return (intersection.boundary() ? 1 : 0);
174  }
175 
176  template< class GG, class II >
177  const typename DGFBoundaryParameter::type &
178  boundaryParameter ( const Intersection< GG, II > & intersection ) const
179  {
180  typedef Dune::Intersection< GG, II > Intersection;
181  const typename Intersection::Entity & entity = intersection.inside();
182 
183  const int face = intersection.indexInInside();
184 
185  const auto& refElem =
186  ReferenceElements< double, dimension >::general( entity.type() );
187  int corners = refElem.size( face, 1, dimension );
188  std :: vector< unsigned int > bound( corners );
189  for( int i=0; i < corners; ++i )
190  {
191  const int k = refElem.subEntity( face, 1, i, dimension );
192  bound[ i ] = factory_.insertionIndex( entity.template subEntity< dimension >( k ) );
193  }
194 
195  DuneGridFormatParser::facemap_t::key_type key( bound, false );
196  const DuneGridFormatParser::facemap_t::const_iterator pos = dgf_.facemap.find( key );
197  if( pos != dgf_.facemap.end() )
198  return dgf_.facemap.find( key )->second.second;
199  else
200  return DGFBoundaryParameter::defaultValue();
201  }
202 
203  template< int codim >
204  int numParameters () const
205  {
206  if( codim == 0 )
207  return dgf_.nofelparams;
208  else if( codim == dimension )
209  return dgf_.nofvtxparams;
210  else
211  return 0;
212  }
213 
214  // return true if boundary parameters found
216  {
217  return dgf_.haveBndParameters;
218  }
219 
220  std::vector< double > &parameter ( const Element &element )
221  {
222  if( numParameters< 0 >() <= 0 )
223  {
224  DUNE_THROW( InvalidStateException,
225  "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
226  }
227  return dgf_.elParams[ factory_.insertionIndex( element ) ];
228  }
229 
230  std::vector< double > &parameter ( const Vertex &vertex )
231  {
232  if( numParameters< dimension >() <= 0 )
233  {
234  DUNE_THROW( InvalidStateException,
235  "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
236  }
237  return dgf_.vtxParams[ factory_.insertionIndex( vertex ) ];
238  }
239 
240  protected:
242  std::istream &file,
243  MPICommunicatorType communicator,
244  const std::string &filename );
245 
246 
247  static Grid* callDirectly( const std::string& gridname,
248  const int rank,
249  const char *filename,
250  MPICommunicatorType communicator )
251  {
252  if( !std::is_same< MPICommunicatorType, No_Comm >::value )
253  {
254  // in parallel runs add rank to filename
255  std :: stringstream tmps;
256  tmps << filename << "." << rank;
257  const std :: string &tmp = tmps.str();
258 
259  // if file exits then use it
260  if( fileExists( tmp.c_str() ) )
261  return new Grid( tmp.c_str(), communicator );
262  }
263 
264  // for rank 0 we also check the normal file name
265  if( rank == 0 )
266  {
267  if( fileExists( filename ) )
268  return new Grid( filename , communicator );
269 
270  // only throw this exception on rank 0 because
271  // for the other ranks we can still create empty grids
272  DUNE_THROW( GridError, "Unable to create " << gridname << " from '"
273  << filename << "'." );
274  }
275  // don't create messages in every proc, this does not work for many cores.
276  //else
277  //{
278  // dwarn << "WARNING: P[" << rank << "]: Creating empty grid!" << std::endl;
279  //}
280 
281  // return empty grid on all other processes
282  return new Grid( communicator );
283  }
284  static bool fileExists ( const char *fileName )
285  {
286  std :: ifstream testfile( fileName );
287  if( !testfile )
288  return false;
289  testfile.close();
290  return true;
291  }
292  static int rank( MPICommunicatorType MPICOMM )
293  {
294  int rank = 0;
295 #if HAVE_MPI
296  MPI_Comm_rank( MPICOMM, &rank );
297 #endif
298  return rank;
299  }
300  static int size( MPICommunicatorType MPICOMM )
301  {
302  int size = 1;
303 #if HAVE_MPI
304  MPI_Comm_size( MPICOMM, &size );
305 #endif
306  return size;
307  }
310  DuneGridFormatParser dgf_;
311  };
312 
313  template <int dim, int dimw, ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm >
314  struct DGFGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > > :
315  public DGFBaseFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
316  {
319  typedef typename BaseType :: MPICommunicatorType MPICommunicatorType;
320  protected:
321  using BaseType :: grid_;
322  using BaseType :: callDirectly;
323  public:
324  explicit DGFGridFactory ( std::istream &input,
325  MPICommunicatorType comm = MPIHelper::getCommunicator() )
326  : BaseType( comm )
327  {
328  input.clear();
329  input.seekg( 0 );
330  if( !input )
331  DUNE_THROW( DGFException, "Error resetting input stream." );
332  generate( input, comm );
333  }
334 
335  explicit DGFGridFactory ( const std::string &filename,
336  MPICommunicatorType comm = MPIHelper::getCommunicator())
337  : BaseType( comm )
338  {
339  std::ifstream input( filename.c_str() );
340  bool fileFound = input.is_open() ;
341  if( fileFound )
342  fileFound = generate( input, comm, filename );
343 
344  if( ! fileFound )
345  {
346  std::stringstream gridname;
347  gridname << "ALUGrid< " << dim << ", " << dimw << ", eltype, ref, comm >";
348  grid_ = callDirectly( gridname.str(), this->rank( comm ), filename.c_str(), comm );
349  }
350  }
351 
352  protected:
353  bool generate( std::istream &file, MPICommunicatorType comm, const std::string &filename = "" );
354  };
355 
356 
357  namespace dgf
358  {
359 
361  : public GridParameterBlock
362  {
363  ALU2dGridParameterBlock( std::istream &in, const bool verbose )
364  : GridParameterBlock( in ),
365  tolerance_( 1e-8 )
366  {
367  if( findtoken( "tolerance" ) )
368  {
369  double x;
370  if( getnextentry(x) )
371  tolerance_ = x;
372  else
373  {
374  if( verbose )
375  {
376  dwarn << "GridParameterBlock: found keyword `tolerance' but no value, "
377  << "defaulting to `" << tolerance_ <<"'!" << std::endl;
378  }
379  }
380  if( tolerance_ <= 0 )
381  DUNE_THROW( DGFException, "Nonpositive tolerance specified!" );
382  }
383  else
384  {
385  if( verbose )
386  {
387  dwarn << "GridParameterBlock: Parameter 'tolerance' not specified, "
388  << "defaulting to `" << tolerance_ <<"'!" << std::endl;
389  }
390  }
391  }
392 
393  double tolerance () const { return tolerance_; }
394 
395  protected:
396  double tolerance_;
397  };
398 
399  } //end namespace dgf
400 
401  template < class G >
404  std::istream &file, MPICommunicatorType communicator,
405  const std::string &filename )
406  {
407  typedef G DGFGridType ;
408 
409  const int dimworld = DGFGridType :: dimensionworld ;
410  const int dimgrid = DGFGridType :: dimension;
411  dgf_.element = ( eltype == simplex) ?
412  DuneGridFormatParser::Simplex :
413  DuneGridFormatParser::Cube ;
414  dgf_.dimgrid = dimgrid;
415  dgf_.dimw = dimworld;
416 
417  const bool isDGF = dgf_.isDuneGridFormat( file );
418  file.seekg( 0 );
419  if( !isDGF )
420  return false;
421 
422  int rank = 0;
423 #if ALU3DGRID_PARALLEL
424  MPI_Comm_rank( communicator, &rank );
425 #endif
426 
427  dgf::GridParameterBlock parameter( file );
428 
429  typedef FieldVector< typename DGFGridType :: ctype, dimworld > CoordinateType ;
430 
431  ALUParallelBlock aluParallelBlock( file );
432  const bool readFromParallelDGF = aluParallelBlock.isactive();
433  bool parallelFileExists = false;
434 
435  std::string newfilename;
436  if (readFromParallelDGF)
437  {
438  bool ok = true;
439  for (int p=0;p<=rank && ok;++p)
440  ok = aluParallelBlock.next(newfilename);
441  if (ok)
442  {
443  parallelFileExists = true;
444  std::ifstream newfile(newfilename.c_str());
445  if ( !newfile )
446  {
447  std::cout << "prozess " << rank << " failed to open file " << newfilename << " ... abort" << std::endl;
448  DUNE_THROW( InvalidStateException, "parallel DGF file could not opend" );
449  }
450  assert( newfile );
451  return generateALUGrid(eltype,newfile,communicator,filename);
452  }
453  }
454 
455  GlobalVertexIndexBlock vertexIndex( file );
456  const bool globalVertexIndexFound = vertexIndex.isactive();
457  if( rank == 0 || globalVertexIndexFound )
458  {
459  if( !dgf_.readDuneGrid( file, dimgrid, dimworld ) )
460  DUNE_THROW( InvalidStateException, "DGF file not recognized on second call." );
461 
462  if( eltype == simplex )
463  {
464  if(dimgrid == 3)
465  dgf_.setOrientation( 2, 3 );
466  else
467  dgf_.setRefinement( 0, 1);
468  }
469 
470  if( parallelFileExists && !globalVertexIndexFound )
471  DUNE_THROW( DGFException, "Parallel DGF file requires GLOBALVERTEXINDEX block." );
472 
473  for( int n = 0; n < dgf_.nofvtx; ++n )
474  {
475  CoordinateType pos;
476  for( int i = 0; i < dimworld; ++i )
477  pos[ i ] = dgf_.vtx[ n ][ i ];
478  if ( !globalVertexIndexFound )
479  factory_.insertVertex( pos );
480  else
481  {
482  int globalIndex;
483  bool ok = vertexIndex.next(globalIndex);
484  if (!ok)
485  DUNE_THROW( DGFException, "Not enough values in GlobalVertexIndex block" );
486  factory_.insertVertex( pos, globalIndex );
487  }
488  }
489 
490  const unsigned int elemTopoId = (eltype == simplex) ?
491  Dune::Impl::SimplexTopology< dimgrid >::type::id : Dune::Impl::CubeTopology< dimgrid >::type::id ;
492  const unsigned int faceTopoId = (eltype == simplex) ?
493  Dune::Impl::SimplexTopology< dimgrid-1 >::type::id : Dune::Impl::CubeTopology< dimgrid-1 >::type::id ;
494  GeometryType elementType( elemTopoId, dimgrid );
495 
496  const int nFaces = (eltype == simplex) ? dimgrid+1 : 2*dimgrid;
497  for( int n = 0; n < dgf_.nofelements; ++n )
498  {
499  factory_.insertElement( elementType, dgf_.elements[ n ] );
500  for( int face = 0; face <nFaces; ++face )
501  {
502  typedef DuneGridFormatParser::facemap_t::key_type Key;
503  typedef DuneGridFormatParser::facemap_t::iterator Iterator;
504 
505  const Key key = ElementFaceUtil::generateFace( dimgrid, dgf_.elements[ n ], face );
506  const Iterator it = dgf_.facemap.find( key );
507  if( it != dgf_.facemap.end() )
508  factory_.insertBoundary( n, face, it->second.first );
509  }
510  }
511 
512  dgf::ProjectionBlock projectionBlock( file, dimworld );
513  const DuneBoundaryProjection< dimworld > *projection
514  = projectionBlock.defaultProjection< dimworld >();
515 
516  if( projection != 0 )
517  factory_.insertBoundaryProjection( *projection );
518 
519  const size_t numBoundaryProjections = projectionBlock.numBoundaryProjections();
520  GeometryType type( faceTopoId, dimgrid-1 );
521  for( size_t i = 0; i < numBoundaryProjections; ++i )
522  {
523 
524  const std::vector< unsigned int > &vertices = projectionBlock.boundaryFace( i );
525  const DuneBoundaryProjection< dimworld > *projection
526  = projectionBlock.boundaryProjection< dimworld >( i );
527  factory_.insertBoundaryProjection( type, vertices, projection );
528  }
529  }
530 
531  typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
532  dgf::PeriodicFaceTransformationBlock trafoBlock( file, dimworld );
533  const int size = trafoBlock.numTransformations();
534  for( int k = 0; k < size; ++k )
535  {
536  const Transformation &trafo = trafoBlock.transformation( k );
537 
538  typename GridFactory::WorldMatrix matrix;
539  for( int i = 0; i < dimworld; ++i )
540  for( int j = 0; j < dimworld; ++j )
541  matrix[ i ][ j ] = trafo.matrix( i, j );
542 
543  typename GridFactory::WorldVector shift;
544  for( int i = 0; i < dimworld; ++i )
545  shift[ i ] = trafo.shift[ i ];
546 
547  factory_.insertFaceTransformation( matrix, shift );
548  }
549 
550  int addMissingBoundariesLocal = (dgf_.nofelements > 0) && dgf_.facemap.empty();
551  int addMissingBoundariesGlobal = addMissingBoundariesLocal;
552 #if ALU3DGRID_PARALLEL
553  MPI_Allreduce( &addMissingBoundariesLocal, &addMissingBoundariesGlobal, 1, MPI_INT, MPI_MAX, communicator );
554 #endif
555 
556  if( !parameter.dumpFileName().empty() )
557  grid_ = factory_.createGrid( addMissingBoundariesGlobal, false, parameter.dumpFileName() );
558  else
559  grid_ = factory_.createGrid( addMissingBoundariesGlobal, true, filename );
560  return true;
561  }
562 
563  template <int dim, int dimw, ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm>
564  inline bool DGFGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
565  ::generate( std::istream &file, MPICommunicatorType communicator, const std::string &filename )
566  {
567  return BaseType :: generateALUGrid( eltype, file, communicator, filename );
568  }
569 
570 
571 
572 } // namespace Dune
573 
574 #endif // else if HAVE_ALUGRID
575 
576 #endif // #ifndef DUNE_ALUGRID_DGF_HH
Definition: alu3dinclude.hh:50
Definition: alu3dinclude.hh:80
ALUGridElementType
basic element types for ALUGrid
Definition: declaration.hh:17
@ simplex
use only simplex elements (i.e., triangles or tetrahedra)
Definition: declaration.hh:18
@ conforming
use conforming bisection refinement
Definition: declaration.hh:25
unstructured parallel implementation of the DUNE grid interface
Definition: alugrid.hh:31
Definition: dgf.hh:121
DuneGridFormatParser dgf_
Definition: dgf.hh:310
Grid * grid() const
Definition: dgf.hh:139
bool wasInserted(const Intersection &intersection) const
Definition: dgf.hh:145
static int size(MPICommunicatorType MPICOMM)
Definition: dgf.hh:300
static Grid * callDirectly(const std::string &gridname, const int rank, const char *filename, MPICommunicatorType communicator)
Definition: dgf.hh:247
bool generateALUGrid(const ALUGridElementType eltype, std::istream &file, MPICommunicatorType communicator, const std::string &filename)
Definition: dgf.hh:403
DGFBaseFactory()
Definition: dgf.hh:129
const DGFBoundaryParameter::type & boundaryParameter(const Intersection< GG, II > &intersection) const
Definition: dgf.hh:178
Grid * grid_
Definition: dgf.hh:308
Grid::template Codim< 0 >::Entity Element
Definition: dgf.hh:125
static bool fileExists(const char *fileName)
Definition: dgf.hh:284
bool haveBoundaryParameters() const
Definition: dgf.hh:215
static const int dimension
Definition: dgf.hh:123
std::vector< double > & parameter(const Element &element)
Definition: dgf.hh:220
Grid::template Codim< dimension >::Entity Vertex
Definition: dgf.hh:126
G Grid
Definition: dgf.hh:122
DGFBaseFactory(MPICommunicatorType comm)
Definition: dgf.hh:134
int boundaryId(const Intersection< GG, II > &intersection) const
Definition: dgf.hh:151
GridFactory factory_
Definition: dgf.hh:309
int numParameters() const
Definition: dgf.hh:204
static int rank(MPICommunicatorType MPICOMM)
Definition: dgf.hh:292
std::vector< double > & parameter(const Vertex &vertex)
Definition: dgf.hh:230
MPIHelper::MPICommunicator MPICommunicatorType
Definition: dgf.hh:124
Dune::GridFactory< Grid > GridFactory
Definition: dgf.hh:127
DGFGridFactory(std::istream &input, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition: dgf.hh:324
DGFBaseFactory< DGFGridType > BaseType
Definition: dgf.hh:318
DGFGridFactory(const std::string &filename, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition: dgf.hh:335
ALUGrid< dim, dimw, eltype, refinementtype, Comm > DGFGridType
Definition: dgf.hh:317
BaseType ::MPICommunicatorType MPICommunicatorType
Definition: dgf.hh:319
double tolerance() const
Definition: dgf.hh:393
double tolerance_
Definition: dgf.hh:396
ALU2dGridParameterBlock(std::istream &in, const bool verbose)
Definition: dgf.hh:363