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 const typename Intersection::Entity & entity = intersection.inside();
156 const int face = intersection.indexInInside();
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 )
164 const int k = refElem.subEntity( face, 1, i,
dimension );
165 bound[ i ] =
factory_.insertionIndex( entity.template subEntity< dimension >( k ) );
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;
173 return (intersection.boundary() ? 1 : 0);
176 template<
class GG,
class II >
177 const typename DGFBoundaryParameter::type &
180 typedef Dune::Intersection< GG, II > Intersection;
181 const typename Intersection::Entity & entity = intersection.inside();
183 const int face = intersection.indexInInside();
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 )
191 const int k = refElem.subEntity( face, 1, i,
dimension );
192 bound[ i ] =
factory_.insertionIndex( entity.template subEntity< dimension >( k ) );
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;
200 return DGFBoundaryParameter::defaultValue();
203 template<
int codim >
207 return dgf_.nofelparams;
209 return dgf_.nofvtxparams;
217 return dgf_.haveBndParameters;
222 if( numParameters< 0 >() <= 0 )
224 DUNE_THROW( InvalidStateException,
225 "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
227 return dgf_.elParams[
factory_.insertionIndex( element ) ];
232 if( numParameters< dimension >() <= 0 )
234 DUNE_THROW( InvalidStateException,
235 "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
237 return dgf_.vtxParams[
factory_.insertionIndex( vertex ) ];
244 const std::string &filename );
249 const char *filename,
252 if( !std::is_same< MPICommunicatorType, No_Comm >::value )
255 std :: stringstream tmps;
256 tmps << filename <<
"." <<
rank;
257 const std :: string &tmp = tmps.str();
261 return new Grid( tmp.c_str(), communicator );
268 return new Grid( filename , communicator );
272 DUNE_THROW( GridError,
"Unable to create " << gridname <<
" from '"
273 << filename <<
"'." );
282 return new Grid( communicator );
286 std :: ifstream testfile( fileName );
296 MPI_Comm_rank( MPICOMM, &
rank );
304 MPI_Comm_size( MPICOMM, &
size );
313 template <
int dim,
int dimw, ALUGr
idElementType eltype, ALUGr
idRefinementType refinementtype,
class Comm >
314 struct DGFGridFactory<
ALUGrid< dim, dimw, eltype, refinementtype, Comm > > :
315 public DGFBaseFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
321 using BaseType :: grid_;
322 using BaseType :: callDirectly;
331 DUNE_THROW( DGFException,
"Error resetting input stream." );
332 generate( input, comm );
339 std::ifstream input( filename.c_str() );
340 bool fileFound = input.is_open() ;
342 fileFound = generate( input, comm, filename );
346 std::stringstream gridname;
347 gridname <<
"ALUGrid< " << dim <<
", " << dimw <<
", eltype, ref, comm >";
348 grid_ = callDirectly( gridname.str(), this->rank( comm ), filename.c_str(), comm );
353 bool generate( std::istream &file, MPICommunicatorType comm,
const std::string &filename =
"" );
361 :
public GridParameterBlock
364 : GridParameterBlock( in ),
367 if( findtoken(
"tolerance" ) )
370 if( getnextentry(x) )
376 dwarn <<
"GridParameterBlock: found keyword `tolerance' but no value, "
377 <<
"defaulting to `" <<
tolerance_ <<
"'!" << std::endl;
381 DUNE_THROW( DGFException,
"Nonpositive tolerance specified!" );
387 dwarn <<
"GridParameterBlock: Parameter 'tolerance' not specified, "
388 <<
"defaulting to `" <<
tolerance_ <<
"'!" << std::endl;
405 const std::string &filename )
407 typedef G DGFGridType ;
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;
417 const bool isDGF = dgf_.isDuneGridFormat( file );
423 #if ALU3DGRID_PARALLEL
424 MPI_Comm_rank( communicator, &rank );
427 dgf::GridParameterBlock parameter( file );
429 typedef FieldVector< typename DGFGridType :: ctype, dimworld > CoordinateType ;
431 ALUParallelBlock aluParallelBlock( file );
432 const bool readFromParallelDGF = aluParallelBlock.isactive();
433 bool parallelFileExists =
false;
435 std::string newfilename;
436 if (readFromParallelDGF)
439 for (
int p=0;p<=rank && ok;++p)
440 ok = aluParallelBlock.next(newfilename);
443 parallelFileExists =
true;
444 std::ifstream newfile(newfilename.c_str());
447 std::cout <<
"prozess " << rank <<
" failed to open file " << newfilename <<
" ... abort" << std::endl;
448 DUNE_THROW( InvalidStateException,
"parallel DGF file could not opend" );
451 return generateALUGrid(eltype,newfile,communicator,filename);
455 GlobalVertexIndexBlock vertexIndex( file );
456 const bool globalVertexIndexFound = vertexIndex.isactive();
457 if( rank == 0 || globalVertexIndexFound )
459 if( !dgf_.readDuneGrid( file, dimgrid, dimworld ) )
460 DUNE_THROW( InvalidStateException,
"DGF file not recognized on second call." );
465 dgf_.setOrientation( 2, 3 );
467 dgf_.setRefinement( 0, 1);
470 if( parallelFileExists && !globalVertexIndexFound )
471 DUNE_THROW( DGFException,
"Parallel DGF file requires GLOBALVERTEXINDEX block." );
473 for(
int n = 0; n < dgf_.nofvtx; ++n )
476 for(
int i = 0; i < dimworld; ++i )
477 pos[ i ] = dgf_.vtx[ n ][ i ];
478 if ( !globalVertexIndexFound )
479 factory_.insertVertex( pos );
483 bool ok = vertexIndex.next(globalIndex);
485 DUNE_THROW( DGFException,
"Not enough values in GlobalVertexIndex block" );
486 factory_.insertVertex( pos, globalIndex );
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 );
496 const int nFaces = (eltype ==
simplex) ? dimgrid+1 : 2*dimgrid;
497 for(
int n = 0; n < dgf_.nofelements; ++n )
499 factory_.insertElement( elementType, dgf_.elements[ n ] );
500 for(
int face = 0; face <nFaces; ++face )
502 typedef DuneGridFormatParser::facemap_t::key_type Key;
503 typedef DuneGridFormatParser::facemap_t::iterator Iterator;
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 );
512 dgf::ProjectionBlock projectionBlock( file, dimworld );
513 const DuneBoundaryProjection< dimworld > *projection
514 = projectionBlock.defaultProjection< dimworld >();
516 if( projection != 0 )
517 factory_.insertBoundaryProjection( *projection );
519 const size_t numBoundaryProjections = projectionBlock.numBoundaryProjections();
520 GeometryType type( faceTopoId, dimgrid-1 );
521 for(
size_t i = 0; i < numBoundaryProjections; ++i )
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 );
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 )
536 const Transformation &trafo = trafoBlock.transformation( k );
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 );
543 typename GridFactory::WorldVector shift;
544 for(
int i = 0; i < dimworld; ++i )
545 shift[ i ] = trafo.shift[ i ];
547 factory_.insertFaceTransformation( matrix, shift );
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 );
556 if( !parameter.dumpFileName().empty() )
557 grid_ = factory_.createGrid( addMissingBoundariesGlobal,
false, parameter.dumpFileName() );
559 grid_ = factory_.createGrid( addMissingBoundariesGlobal,
true, filename );
563 template <
int dim,
int dimw, ALUGr
idElementType eltype, ALUGr
idRefinementType refinementtype,
class Comm>
564 inline bool DGFGridFactory< ALUGrid< dim, dimw, eltype, refinementtype, Comm > >
567 return BaseType :: generateALUGrid( eltype, file, communicator, filename );
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
static double refineWeight()
Definition: dgf.hh:109
static int refineStepsForHalf()
Definition: dgf.hh:108
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