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) ?
507 GeometryType::simplex :
508 GeometryType::cube, dimgrid );
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) ?
537 GeometryType::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 );