1 #ifndef DUNE_ALUGRID_DGF_HH
2 #define DUNE_ALUGRID_DGF_HH
8 #include <dune/grid/io/file/dgfparser/dgfalu.hh>
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>
29 class GlobalVertexIndexBlock
30 :
public dgf::BasicBlock
35 GlobalVertexIndexBlock ( std :: istream &in )
36 : dgf::BasicBlock( in,
"GlobalVertexIndex" ),
40 bool next (
int &index )
44 return (goodline =
false);
46 if( !getnextentry( index ) )
48 DUNE_THROW ( DGFException,
"Error in " << *
this <<
": "
49 <<
"Wrong global vertex indices " );
51 return (goodline =
true);
66 class ALUParallelBlock
67 :
public dgf::BasicBlock
72 ALUParallelBlock ( std :: istream &in )
73 : dgf::BasicBlock( in,
"ALUParallel" ),
77 bool next ( std::string &name )
81 return (goodline =
false);
83 if( !getnextentry( name ) )
85 DUNE_THROW ( DGFException,
"Error in " << *
this <<
": "
86 <<
"Wrong global vertex indices " );
88 return (goodline =
true);
105 template<
int dimg,
int dimw, ALUGr
idElementType eltype, ALUGr
idRefinementType refinementtype,
class Comm >
106 struct DGFGridInfo<
Dune::
ALUGrid< dimg, dimw, eltype, refinementtype, Comm > >
125 typedef typename Grid::template Codim<0>::Entity
Element;
126 typedef typename Grid::template Codim<dimension>::Entity
Vertex;
144 template<
class Intersection >
147 return factory_.wasInserted( intersection );
150 template<
class GG,
class II >
151 int boundaryId (
const Intersection< GG, II > & intersection )
const
153 typedef Dune::Intersection< GG, II > Intersection;
154 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
155 const typename Intersection::Entity & entity = intersection.inside();
157 typename Intersection::EntityPointer inside = intersection.inside();
158 const typename Intersection::Entity & entity = *inside;
160 const int face = intersection.indexInInside();
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 )
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 ) );
172 bound[ i ] =
factory_.insertionIndex( *entity.template subEntity< dimension >( k ) );
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;
181 return (intersection.boundary() ? 1 : 0);
184 template<
class GG,
class II >
185 const typename DGFBoundaryParameter::type &
188 typedef Dune::Intersection< GG, II > Intersection;
189 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
190 const typename Intersection::Entity & entity = intersection.inside();
192 typename Intersection::EntityPointer inside = intersection.inside();
193 const typename Intersection::Entity & entity = *inside;
195 const int face = intersection.indexInInside();
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 )
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 ) );
207 bound[ i ] =
factory_.insertionIndex( *entity.template subEntity< dimension >( k ) );
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;
216 return DGFBoundaryParameter::defaultValue();
219 template<
int codim >
223 return dgf_.nofelparams;
225 return dgf_.nofvtxparams;
233 return dgf_.haveBndParameters;
238 if( numParameters< 0 >() <= 0 )
240 DUNE_THROW( InvalidStateException,
241 "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
243 return dgf_.elParams[
factory_.insertionIndex( element ) ];
248 if( numParameters< dimension >() <= 0 )
250 DUNE_THROW( InvalidStateException,
251 "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
253 return dgf_.vtxParams[
factory_.insertionIndex( vertex ) ];
260 const std::string &filename );
265 const char *filename,
268 if( !std::is_same< MPICommunicatorType, No_Comm >::value )
271 std :: stringstream tmps;
272 tmps << filename <<
"." <<
rank;
273 const std :: string &tmp = tmps.str();
277 return new Grid( tmp.c_str(), communicator );
284 return new Grid( filename , communicator );
288 DUNE_THROW( GridError,
"Unable to create " << gridname <<
" from '"
289 << filename <<
"'." );
298 return new Grid( communicator );
302 std :: ifstream testfile( fileName );
312 MPI_Comm_rank( MPICOMM, &
rank );
320 MPI_Comm_size( MPICOMM, &
size );
329 template <
int dim,
int dimw, ALUGr
idElementType eltype, ALUGr
idRefinementType refinementtype,
class Comm >
330 struct DGFGridFactory<
ALUGrid< dim, dimw, eltype, refinementtype, Comm > > :
331 public DGFBaseFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
337 using BaseType :: grid_;
338 using BaseType :: callDirectly;
347 DUNE_THROW( DGFException,
"Error resetting input stream." );
348 generate( input, comm );
355 std::ifstream input( filename.c_str() );
356 bool fileFound = input.is_open() ;
358 fileFound = generate( input, comm, filename );
362 std::stringstream gridname;
363 gridname <<
"ALUGrid< " << dim <<
", " << dimw <<
", eltype, ref, comm >";
364 grid_ = callDirectly( gridname.str(), this->rank( comm ), filename.c_str(), comm );
369 bool generate( std::istream &file, MPICommunicatorType comm,
const std::string &filename =
"" );
377 :
public GridParameterBlock
380 : GridParameterBlock( in ),
383 if( findtoken(
"tolerance" ) )
386 if( getnextentry(x) )
392 dwarn <<
"GridParameterBlock: found keyword `tolerance' but no value, "
393 <<
"defaulting to `" <<
tolerance_ <<
"'!" << std::endl;
397 DUNE_THROW( DGFException,
"Nonpositive tolerance specified!" );
403 dwarn <<
"GridParameterBlock: Parameter 'tolerance' not specified, "
404 <<
"defaulting to `" <<
tolerance_ <<
"'!" << std::endl;
421 const std::string &filename )
423 typedef G DGFGridType ;
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;
433 const bool isDGF = dgf_.isDuneGridFormat( file );
439 #if ALU3DGRID_PARALLEL
440 MPI_Comm_rank( communicator, &rank );
443 dgf::GridParameterBlock parameter( file );
445 typedef FieldVector< typename DGFGridType :: ctype, dimworld > CoordinateType ;
447 ALUParallelBlock aluParallelBlock( file );
448 const bool readFromParallelDGF = aluParallelBlock.isactive();
449 bool parallelFileExists =
false;
451 std::string newfilename;
452 if (readFromParallelDGF)
455 for (
int p=0;p<=rank && ok;++p)
456 ok = aluParallelBlock.next(newfilename);
459 parallelFileExists =
true;
460 std::ifstream newfile(newfilename.c_str());
463 std::cout <<
"prozess " << rank <<
" failed to open file " << newfilename <<
" ... abort" << std::endl;
464 DUNE_THROW( InvalidStateException,
"parallel DGF file could not opend" );
467 return generateALUGrid(eltype,newfile,communicator,filename);
471 GlobalVertexIndexBlock vertexIndex( file );
472 const bool globalVertexIndexFound = vertexIndex.isactive();
473 if( rank == 0 || globalVertexIndexFound )
475 if( !dgf_.readDuneGrid( file, dimgrid, dimworld ) )
476 DUNE_THROW( InvalidStateException,
"DGF file not recognized on second call." );
481 dgf_.setOrientation( 2, 3 );
483 dgf_.setRefinement( 0, 1);
486 if( parallelFileExists && !globalVertexIndexFound )
487 DUNE_THROW( DGFException,
"Parallel DGF file requires GLOBALVERTEXINDEX block." );
489 for(
int n = 0; n < dgf_.nofvtx; ++n )
492 for(
int i = 0; i < dimworld; ++i )
493 pos[ i ] = dgf_.vtx[ n ][ i ];
494 if ( !globalVertexIndexFound )
495 factory_.insertVertex( pos );
499 bool ok = vertexIndex.next(globalIndex);
501 DUNE_THROW( DGFException,
"Not enough values in GlobalVertexIndex block" );
502 factory_.insertVertex( pos, globalIndex );
506 GeometryType elementType( (eltype ==
simplex) ?
510 const int nFaces = (eltype ==
simplex) ? dimgrid+1 : 2*dimgrid;
511 for(
int n = 0; n < dgf_.nofelements; ++n )
513 factory_.insertElement( elementType, dgf_.elements[ n ] );
514 for(
int face = 0; face <nFaces; ++face )
516 typedef DuneGridFormatParser::facemap_t::key_type Key;
517 typedef DuneGridFormatParser::facemap_t::iterator Iterator;
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 );
526 dgf::ProjectionBlock projectionBlock( file, dimworld );
527 const DuneBoundaryProjection< dimworld > *projection
528 = projectionBlock.defaultProjection< dimworld >();
530 if( projection != 0 )
531 factory_.insertBoundaryProjection( *projection );
533 const size_t numBoundaryProjections = projectionBlock.numBoundaryProjections();
534 for(
size_t i = 0; i < numBoundaryProjections; ++i )
536 GeometryType type( (eltype ==
simplex) ?
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 );
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 )
553 const Transformation &trafo = trafoBlock.transformation( k );
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 );
560 typename GridFactory::WorldVector shift;
561 for(
int i = 0; i < dimworld; ++i )
562 shift[ i ] = trafo.shift[ i ];
564 factory_.insertFaceTransformation( matrix, shift );
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 );
573 if( !parameter.dumpFileName().empty() )
574 grid_ = factory_.createGrid( addMissingBoundariesGlobal,
false, parameter.dumpFileName() );
576 grid_ = factory_.createGrid( addMissingBoundariesGlobal,
true, filename );
580 template <
int dim,
int dimw, ALUGr
idElementType eltype, ALUGr
idRefinementType refinementtype,
class Comm>
581 inline bool DGFGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
584 return BaseType :: generateALUGrid( eltype, file, communicator, filename );
591 #endif // else if HAVE_ALUGRID
593 #endif // #ifndef DUNE_ALUGRID_DGF_HH