dune-alugrid 3.0.0
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
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
20namespace 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#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
155 const typename Intersection::Entity & entity = intersection.inside();
156#else
157 typename Intersection::EntityPointer inside = intersection.inside();
158 const typename Intersection::Entity & entity = *inside;
159#endif
160 const int face = intersection.indexInInside();
161
162 const ReferenceElement< double, dimension > & refElem =
163 ReferenceElements< double, dimension >::general( entity.type() );
164 int corners = refElem.size( face, 1, dimension );
165 std :: vector< unsigned int > bound( corners );
166 for( int i=0; i < corners; ++i )
167 {
168 const int k = refElem.subEntity( face, 1, i, dimension );
169#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
170 bound[ i ] = factory_.insertionIndex( entity.template subEntity< dimension >( k ) );
171#else
172 bound[ i ] = factory_.insertionIndex( *entity.template subEntity< dimension >( k ) );
173#endif
174 }
175
176 DuneGridFormatParser::facemap_t::key_type key( bound, false );
177 const DuneGridFormatParser::facemap_t::const_iterator pos = dgf_.facemap.find( key );
178 if( pos != dgf_.facemap.end() )
179 return dgf_.facemap.find( key )->second.first;
180 else
181 return (intersection.boundary() ? 1 : 0);
182 }
183
184 template< class GG, class II >
185 const typename DGFBoundaryParameter::type &
186 boundaryParameter ( const Intersection< GG, II > & intersection ) const
187 {
188 typedef Dune::Intersection< GG, II > Intersection;
189#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
190 const typename Intersection::Entity & entity = intersection.inside();
191#else
192 typename Intersection::EntityPointer inside = intersection.inside();
193 const typename Intersection::Entity & entity = *inside;
194#endif
195 const int face = intersection.indexInInside();
196
197 const ReferenceElement< double, dimension > & refElem =
198 ReferenceElements< double, dimension >::general( entity.type() );
199 int corners = refElem.size( face, 1, dimension );
200 std :: vector< unsigned int > bound( corners );
201 for( int i=0; i < corners; ++i )
202 {
203 const int k = refElem.subEntity( face, 1, i, dimension );
204#if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
205 bound[ i ] = factory_.insertionIndex( entity.template subEntity< dimension >( k ) );
206#else
207 bound[ i ] = factory_.insertionIndex( *entity.template subEntity< dimension >( k ) );
208#endif
209 }
210
211 DuneGridFormatParser::facemap_t::key_type key( bound, false );
212 const DuneGridFormatParser::facemap_t::const_iterator pos = dgf_.facemap.find( key );
213 if( pos != dgf_.facemap.end() )
214 return dgf_.facemap.find( key )->second.second;
215 else
216 return DGFBoundaryParameter::defaultValue();
217 }
218
219 template< int codim >
220 int numParameters () const
221 {
222 if( codim == 0 )
223 return dgf_.nofelparams;
224 else if( codim == dimension )
225 return dgf_.nofvtxparams;
226 else
227 return 0;
228 }
229
230 // return true if boundary parameters found
232 {
233 return dgf_.haveBndParameters;
234 }
235
236 std::vector< double > &parameter ( const Element &element )
237 {
238 if( numParameters< 0 >() <= 0 )
239 {
240 DUNE_THROW( InvalidStateException,
241 "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
242 }
243 return dgf_.elParams[ factory_.insertionIndex( element ) ];
244 }
245
246 std::vector< double > &parameter ( const Vertex &vertex )
247 {
248 if( numParameters< dimension >() <= 0 )
249 {
250 DUNE_THROW( InvalidStateException,
251 "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
252 }
253 return dgf_.vtxParams[ factory_.insertionIndex( vertex ) ];
254 }
255
256 protected:
258 std::istream &file,
259 MPICommunicatorType communicator,
260 const std::string &filename );
261
262
263 static Grid* callDirectly( const std::string& gridname,
264 const int rank,
265 const char *filename,
266 MPICommunicatorType communicator )
267 {
268 if( !std::is_same< MPICommunicatorType, No_Comm >::value )
269 {
270 // in parallel runs add rank to filename
271 std :: stringstream tmps;
272 tmps << filename << "." << rank;
273 const std :: string &tmp = tmps.str();
274
275 // if file exits then use it
276 if( fileExists( tmp.c_str() ) )
277 return new Grid( tmp.c_str(), communicator );
278 }
279
280 // for rank 0 we also check the normal file name
281 if( rank == 0 )
282 {
283 if( fileExists( filename ) )
284 return new Grid( filename , communicator );
285
286 // only throw this exception on rank 0 because
287 // for the other ranks we can still create empty grids
288 DUNE_THROW( GridError, "Unable to create " << gridname << " from '"
289 << filename << "'." );
290 }
291 // don't create messages in every proc, this does not work for many cores.
292 //else
293 //{
294 // dwarn << "WARNING: P[" << rank << "]: Creating empty grid!" << std::endl;
295 //}
296
297 // return empty grid on all other processes
298 return new Grid( communicator );
299 }
300 static bool fileExists ( const char *fileName )
301 {
302 std :: ifstream testfile( fileName );
303 if( !testfile )
304 return false;
305 testfile.close();
306 return true;
307 }
308 static int rank( MPICommunicatorType MPICOMM )
309 {
310 int rank = 0;
311#if HAVE_MPI
312 MPI_Comm_rank( MPICOMM, &rank );
313#endif
314 return rank;
315 }
316 static int size( MPICommunicatorType MPICOMM )
317 {
318 int size = 1;
319#if HAVE_MPI
320 MPI_Comm_size( MPICOMM, &size );
321#endif
322 return size;
323 }
326 DuneGridFormatParser dgf_;
327 };
328
329 template <int dim, int dimw, ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm >
330 struct DGFGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > > :
331 public DGFBaseFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
332 {
335 typedef typename BaseType :: MPICommunicatorType MPICommunicatorType;
336 protected:
337 using BaseType :: grid_;
338 using BaseType :: callDirectly;
339 public:
340 explicit DGFGridFactory ( std::istream &input,
341 MPICommunicatorType comm = MPIHelper::getCommunicator() )
342 : BaseType( comm )
343 {
344 input.clear();
345 input.seekg( 0 );
346 if( !input )
347 DUNE_THROW( DGFException, "Error resetting input stream." );
348 generate( input, comm );
349 }
350
351 explicit DGFGridFactory ( const std::string &filename,
352 MPICommunicatorType comm = MPIHelper::getCommunicator())
353 : BaseType( comm )
354 {
355 std::ifstream input( filename.c_str() );
356 bool fileFound = input.is_open() ;
357 if( fileFound )
358 fileFound = generate( input, comm, filename );
359
360 if( ! fileFound )
361 {
362 std::stringstream gridname;
363 gridname << "ALUGrid< " << dim << ", " << dimw << ", eltype, ref, comm >";
364 grid_ = callDirectly( gridname.str(), this->rank( comm ), filename.c_str(), comm );
365 }
366 }
367
368 protected:
369 bool generate( std::istream &file, MPICommunicatorType comm, const std::string &filename = "" );
370 };
371
372
373 namespace dgf
374 {
375
377 : public GridParameterBlock
378 {
379 ALU2dGridParameterBlock( std::istream &in, const bool verbose )
380 : GridParameterBlock( in ),
381 tolerance_( 1e-8 )
382 {
383 if( findtoken( "tolerance" ) )
384 {
385 double x;
386 if( getnextentry(x) )
387 tolerance_ = x;
388 else
389 {
390 if( verbose )
391 {
392 dwarn << "GridParameterBlock: found keyword `tolerance' but no value, "
393 << "defaulting to `" << tolerance_ <<"'!" << std::endl;
394 }
395 }
396 if( tolerance_ <= 0 )
397 DUNE_THROW( DGFException, "Nonpositive tolerance specified!" );
398 }
399 else
400 {
401 if( verbose )
402 {
403 dwarn << "GridParameterBlock: Parameter 'tolerance' not specified, "
404 << "defaulting to `" << tolerance_ <<"'!" << std::endl;
405 }
406 }
407 }
408
409 double tolerance () const { return tolerance_; }
410
411 protected:
413 };
414
415 } //end namespace dgf
416
417 template < class G >
420 std::istream &file, MPICommunicatorType communicator,
421 const std::string &filename )
422 {
423 typedef G DGFGridType ;
424
425 const int dimworld = DGFGridType :: dimensionworld ;
426 const int dimgrid = DGFGridType :: dimension;
427 dgf_.element = ( eltype == simplex) ?
428 DuneGridFormatParser::Simplex :
429 DuneGridFormatParser::Cube ;
430 dgf_.dimgrid = dimgrid;
431 dgf_.dimw = dimworld;
432
433 const bool isDGF = dgf_.isDuneGridFormat( file );
434 file.seekg( 0 );
435 if( !isDGF )
436 return false;
437
438 int rank = 0;
439#if ALU3DGRID_PARALLEL
440 MPI_Comm_rank( communicator, &rank );
441#endif
442
443 dgf::GridParameterBlock parameter( file );
444
445 typedef FieldVector< typename DGFGridType :: ctype, dimworld > CoordinateType ;
446
447 ALUParallelBlock aluParallelBlock( file );
448 const bool readFromParallelDGF = aluParallelBlock.isactive();
449 bool parallelFileExists = false;
450
451 std::string newfilename;
452 if (readFromParallelDGF)
453 {
454 bool ok = true;
455 for (int p=0;p<=rank && ok;++p)
456 ok = aluParallelBlock.next(newfilename);
457 if (ok)
458 {
459 parallelFileExists = true;
460 std::ifstream newfile(newfilename.c_str());
461 if ( !newfile )
462 {
463 std::cout << "prozess " << rank << " failed to open file " << newfilename << " ... abort" << std::endl;
464 DUNE_THROW( InvalidStateException, "parallel DGF file could not opend" );
465 }
466 assert( newfile );
467 return generateALUGrid(eltype,newfile,communicator,filename);
468 }
469 }
470
471 GlobalVertexIndexBlock vertexIndex( file );
472 const bool globalVertexIndexFound = vertexIndex.isactive();
473 if( rank == 0 || globalVertexIndexFound )
474 {
475 if( !dgf_.readDuneGrid( file, dimgrid, dimworld ) )
476 DUNE_THROW( InvalidStateException, "DGF file not recognized on second call." );
477
478 if( eltype == simplex )
479 {
480 if(dimgrid == 3)
481 dgf_.setOrientation( 2, 3 );
482 else
483 dgf_.setRefinement( 0, 1);
484 }
485
486 if( parallelFileExists && !globalVertexIndexFound )
487 DUNE_THROW( DGFException, "Parallel DGF file requires GLOBALVERTEXINDEX block." );
488
489 for( int n = 0; n < dgf_.nofvtx; ++n )
490 {
491 CoordinateType pos;
492 for( int i = 0; i < dimworld; ++i )
493 pos[ i ] = dgf_.vtx[ n ][ i ];
494 if ( !globalVertexIndexFound )
495 factory_.insertVertex( pos );
496 else
497 {
498 int globalIndex;
499 bool ok = vertexIndex.next(globalIndex);
500 if (!ok)
501 DUNE_THROW( DGFException, "Not enough values in GlobalVertexIndex block" );
502 factory_.insertVertex( pos, globalIndex );
503 }
504 }
505
506 GeometryType elementType( (eltype == simplex) ?
507 GeometryType::simplex :
508 GeometryType::cube, dimgrid );
509
510 const int nFaces = (eltype == simplex) ? dimgrid+1 : 2*dimgrid;
511 for( int n = 0; n < dgf_.nofelements; ++n )
512 {
513 factory_.insertElement( elementType, dgf_.elements[ n ] );
514 for( int face = 0; face <nFaces; ++face )
515 {
516 typedef DuneGridFormatParser::facemap_t::key_type Key;
517 typedef DuneGridFormatParser::facemap_t::iterator Iterator;
518
519 const Key key = ElementFaceUtil::generateFace( dimgrid, dgf_.elements[ n ], face );
520 const Iterator it = dgf_.facemap.find( key );
521 if( it != dgf_.facemap.end() )
522 factory_.insertBoundary( n, face, it->second.first );
523 }
524 }
525
526 dgf::ProjectionBlock projectionBlock( file, dimworld );
527 const DuneBoundaryProjection< dimworld > *projection
528 = projectionBlock.defaultProjection< dimworld >();
529
530 if( projection != 0 )
531 factory_.insertBoundaryProjection( *projection );
532
533 const size_t numBoundaryProjections = projectionBlock.numBoundaryProjections();
534 for( size_t i = 0; i < numBoundaryProjections; ++i )
535 {
536 GeometryType type( (eltype == simplex) ?
537 GeometryType::simplex :
538 GeometryType::cube,
539 dimgrid-1);
540
541 const std::vector< unsigned int > &vertices = projectionBlock.boundaryFace( i );
542 const DuneBoundaryProjection< dimworld > *projection
543 = projectionBlock.boundaryProjection< dimworld >( i );
544 factory_.insertBoundaryProjection( type, vertices, projection );
545 }
546 }
547
548 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
549 dgf::PeriodicFaceTransformationBlock trafoBlock( file, dimworld );
550 const int size = trafoBlock.numTransformations();
551 for( int k = 0; k < size; ++k )
552 {
553 const Transformation &trafo = trafoBlock.transformation( k );
554
555 typename GridFactory::WorldMatrix matrix;
556 for( int i = 0; i < dimworld; ++i )
557 for( int j = 0; j < dimworld; ++j )
558 matrix[ i ][ j ] = trafo.matrix( i, j );
559
560 typename GridFactory::WorldVector shift;
561 for( int i = 0; i < dimworld; ++i )
562 shift[ i ] = trafo.shift[ i ];
563
564 factory_.insertFaceTransformation( matrix, shift );
565 }
566
567 int addMissingBoundariesLocal = (dgf_.nofelements > 0) && dgf_.facemap.empty();
568 int addMissingBoundariesGlobal = addMissingBoundariesLocal;
569#if ALU3DGRID_PARALLEL
570 MPI_Allreduce( &addMissingBoundariesLocal, &addMissingBoundariesGlobal, 1, MPI_INT, MPI_MAX, communicator );
571#endif
572
573 if( !parameter.dumpFileName().empty() )
574 grid_ = factory_.createGrid( addMissingBoundariesGlobal, false, parameter.dumpFileName() );
575 else
576 grid_ = factory_.createGrid( addMissingBoundariesGlobal, true, filename );
577 return true;
578 }
579
580 template <int dim, int dimw, ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm>
581 inline bool DGFGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
582 ::generate( std::istream &file, MPICommunicatorType communicator, const std::string &filename )
583 {
584 return BaseType :: generateALUGrid( eltype, file, communicator, filename );
585 }
586
587
588
589} // namespace Dune
590
591#endif // else if HAVE_ALUGRID
592
593#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:326
bool wasInserted(const Intersection &intersection) const
Definition dgf.hh:145
static int size(MPICommunicatorType MPICOMM)
Definition dgf.hh:316
bool generateALUGrid(const ALUGridElementType eltype, std::istream &file, MPICommunicatorType communicator, const std::string &filename)
Definition dgf.hh:419
DGFBaseFactory()
Definition dgf.hh:129
Grid * grid_
Definition dgf.hh:324
Grid::template Codim< 0 >::Entity Element
Definition dgf.hh:125
std::vector< double > & parameter(const Vertex &vertex)
Definition dgf.hh:246
static bool fileExists(const char *fileName)
Definition dgf.hh:300
bool haveBoundaryParameters() const
Definition dgf.hh:231
static const int dimension
Definition dgf.hh:123
Grid * grid() const
Definition dgf.hh:139
Grid::template Codim< dimension >::Entity Vertex
Definition dgf.hh:126
G Grid
Definition dgf.hh:122
DGFBaseFactory(MPICommunicatorType comm)
Definition dgf.hh:134
std::vector< double > & parameter(const Element &element)
Definition dgf.hh:236
int boundaryId(const Intersection< GG, II > &intersection) const
Definition dgf.hh:151
GridFactory factory_
Definition dgf.hh:325
int numParameters() const
Definition dgf.hh:220
static int rank(MPICommunicatorType MPICOMM)
Definition dgf.hh:308
const DGFBoundaryParameter::type & boundaryParameter(const Intersection< GG, II > &intersection) const
Definition dgf.hh:186
static Grid * callDirectly(const std::string &gridname, const int rank, const char *filename, MPICommunicatorType communicator)
Definition dgf.hh:263
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:340
DGFBaseFactory< DGFGridType > BaseType
Definition dgf.hh:334
DGFGridFactory(const std::string &filename, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition dgf.hh:351
ALUGrid< dim, dimw, eltype, refinementtype, Comm > DGFGridType
Definition dgf.hh:333
BaseType::MPICommunicatorType MPICommunicatorType
Definition dgf.hh:335
double tolerance() const
Definition dgf.hh:409
double tolerance_
Definition dgf.hh:412
ALU2dGridParameterBlock(std::istream &in, const bool verbose)
Definition dgf.hh:379