3 #ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH
4 #define DUNE_POLYHEDRALGRID_DGFPARSER_HH
11 #include <dune/common/typetraits.hh>
12 #include <dune/common/version.hh>
14 #include <dune/geometry/referenceelements.hh>
16 #include <dune/grid/common/gridfactory.hh>
17 #include <dune/grid/io/file/dgfparser/dgfparser.hh>
19 #if DUNE_VERSION_NEWER(DUNE_GRID,2,7)
20 #include <dune/grid/io/file/dgfparser/blocks/polyhedron.hh>
23 #include <opm/grid/polyhedralgrid/gridfactory.hh>
26 #include <opm/parser/eclipse/Parser/ParseContext.hpp>
27 #include <opm/parser/eclipse/Parser/Parser.hpp>
28 #include <opm/parser/eclipse/Deck/Deck.hpp>
37 #if ! DUNE_VERSION_NEWER(DUNE_GRID,2,7)
38 namespace PolyhedralGrid
47 PolygonBlock ( std::istream &in,
int numVtx,
int vtxOfs )
48 : BasicBlock( in,
"Polygon" ), vtxBegin_( vtxOfs ), vtxEnd_( vtxOfs + numVtx )
51 int get ( std::vector< std::vector< int > > &polygons )
54 std::vector< int > polygon;
55 while( getnextline() )
58 for(
int vtxIdx; getnextentry( vtxIdx ); )
60 if( (vtxBegin_ > vtxIdx) || (vtxIdx >= vtxEnd_) )
61 DUNE_THROW( DGFException,
"Error in " << *
this <<
": Invalid vertex index (" << vtxIdx <<
" not int [" << vtxBegin_ <<
", " << vtxEnd_ <<
"[)" );
62 polygon.push_back( vtxIdx - vtxBegin_ );
65 polygons.push_back( polygon );
67 return polygons.size();
71 int vtxBegin_, vtxEnd_;
83 : BasicBlock( in,
"Polyhedron" ), numPolys_( numPolys )
86 int get ( std::vector< std::vector< int > > &polyhedra )
89 std::vector< int > polyhedron;
91 while( getnextline() )
94 for(
int polyIdx; getnextentry( polyIdx ); )
96 if( (polyIdx < 0) || (polyIdx > numPolys_) )
97 DUNE_THROW( DGFException,
"Error in " << *
this <<
": Invalid polygon index (" << polyIdx <<
" not int [0, " << numPolys_ <<
"])" );
99 minPolyId = std::min( minPolyId, polyIdx );
100 polyhedron.push_back( polyIdx );
103 polyhedra.push_back( polyhedron );
109 const size_t polySize = polyhedra.size();
110 for(
size_t i=0; i<polySize; ++i )
112 const size_t pSize = polyhedra[ i ].size();
113 for(
size_t j=0; j<pSize; ++j )
115 polyhedra[ i ][ j ] -= minPolyId;
119 return polyhedra.size();
128 using PolyhedralGrid :: PolygonBlock;
129 using PolyhedralGrid :: PolyhedronBlock;
139 template<
int dim,
int dimworld,
class coord_t >
144 const static int dimension = Grid::dimension;
145 typedef MPIHelper::MPICommunicator MPICommunicator;
146 typedef typename Grid::template Codim<0>::Entity Element;
147 typedef typename Grid::template Codim<dimension>::Entity Vertex;
149 explicit DGFGridFactory ( std::istream &input, MPICommunicator = MPIHelper::getCommunicator() )
156 DUNE_THROW( DGFException,
"Error resetting input stream" );
160 explicit DGFGridFactory (
const std::string &filename, MPICommunicator = MPIHelper::getCommunicator() )
164 std::ifstream input( filename );
166 DUNE_THROW( DGFException,
"Macrofile '" << filename <<
"' not found" );
169 if( !DuneGridFormatParser::isDuneGridFormat( input ) )
172 const auto deck = parser.parseString( filename );
173 std::vector<double> porv;
175 gridPtr_.reset(
new Grid( deck, porv ) );
190 grid_ = gridPtr_.release();
195 template<
class Intersection >
196 bool wasInserted (
const Intersection& )
const
201 template<
class Intersection >
202 int boundaryId (
const Intersection& )
const
207 bool haveBoundaryParameters ()
const {
return false; }
209 template<
int codim >
210 int numParameters ()
const
216 template<
class Intersection >
217 const typename DGFBoundaryParameter::type &
218 boundaryParameter (
const Intersection& )
const
220 return DGFBoundaryParameter::defaultValue();;
223 template<
class Entity >
224 std::vector< double > ¶meter (
const Entity& )
226 static std::vector< double > dummy;
231 int readVertices ( std::istream &input, std::vector< std::vector< double > > &vertices )
233 int dimWorld = Grid::dimensionworld ;
234 dgf::VertexBlock vtxBlock( input, dimWorld );
235 if( !vtxBlock.isactive() )
236 DUNE_THROW( DGFException,
"Vertex block not found" );
238 vtxBlock.get( vertices, vtxParams_, numVtxParams_ );
239 return vtxBlock.offset();
242 std::vector< std::vector< int > > readPolygons ( std::istream &input,
int numVtx,
int vtxOfs )
244 dgf::PolygonBlock polygonBlock( input, numVtx, vtxOfs );
245 if( !polygonBlock.isactive() )
246 DUNE_THROW( DGFException,
"Polygon block not found" );
248 std::vector< std::vector< int > > polygons;
249 polygonBlock.get( polygons );
253 std::vector< std::vector< int > > readPolyhedra ( std::istream &input,
int numPolygons )
255 dgf::PolyhedronBlock polyhedronBlock( input, numPolygons );
256 std::vector< std::vector< int > > polyhedra;
257 if( polyhedronBlock.isactive() )
259 polyhedronBlock.get( polyhedra );
264 template<
class Iterator >
265 void copy ( Iterator begin, Iterator end,
double *dest )
267 for( ; begin != end; ++begin )
268 dest = std::copy( begin->begin(), begin->end(), dest );
271 template<
class Iterator >
272 void copy ( Iterator begin, Iterator end,
int *dest,
int *offset )
275 for( ; begin != end; ++begin )
278 size += begin->size();
279 dest = std::copy( begin->begin(), begin->end(), dest );
284 void generate ( std::istream &input )
287 dgf::IntervalBlock intervalBlock( input );
288 if( intervalBlock.isactive() )
290 if( intervalBlock.numIntervals() != 1 )
291 DUNE_THROW( DGFException,
"Currently, CpGrid can only handle 1 interval block." );
293 if( intervalBlock.dimw() != dimworld )
294 DUNE_THROW( DGFException,
"CpGrid cannot handle an interval of dimension " << intervalBlock.dimw() <<
"." );
295 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
297 std::vector< double > spacing( dimworld );
298 for(
int i=0; i<dimworld; ++i )
299 spacing[ i ] = (interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ]) / interval.n[ i ];
301 gridPtr_.reset(
new Grid( interval.n, spacing ) );
306 typedef std::vector< std::vector< double > > CoordinateVectorType;
307 CoordinateVectorType nodes;
309 typedef std::vector< std::vector< int > > IndexVectorType;
310 IndexVectorType faces;
311 IndexVectorType cells;
313 const int vtxOfs = readVertices( input, nodes );
315 faces = readPolygons ( input, nodes.size(), vtxOfs );
316 cells = readPolyhedra( input, faces.size() );
320 DUNE_THROW( DGFException,
"Polyhedron block not found!" );
323 typedef GridFactory< Grid > GridFactoryType;
324 typedef typename GridFactoryType :: Coordinate Coordinate ;
326 GridFactoryType gridFactory;
328 const int nNodes = nodes.size();
329 Coordinate node( 0 );
330 for(
int i=0; i<nNodes; ++i )
332 for(
int d=0; d<Coordinate::dimension; ++d )
333 node[ d ] = nodes[ i ][ d ];
335 gridFactory.insertVertex( node );
341 type = Dune::GeometryTypes::none(Grid::dimension-1);
342 std::vector< unsigned int > numbers;
344 const int nFaces = faces.size();
345 for(
int i = 0; i < nFaces; ++ i )
348 std::vector<int>& face = faces[ i ];
349 numbers.resize( face.size() );
350 std::copy( face.begin(), face.end(), numbers.begin() );
351 gridFactory.insertElement( type, numbers );
358 type = Dune::GeometryTypes::none(Grid::dimension);
360 const int nCells = cells.size();
361 for(
int i = 0; i < nCells; ++ i )
364 std::vector<int>& cell = cells[ i ];
365 numbers.resize( cell.size() );
366 std::copy( cell.begin(), cell.end(), numbers.begin() );
367 gridFactory.insertElement( type, numbers );
371 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7)
372 gridPtr_ = gridFactory.createGrid();
374 gridPtr_.reset( gridFactory.createGrid() );
383 nodes.resize( dgf.nofvtx );
385 std::copy( dgf.vtx.begin(), dgf.vtx.end(), nodes.begin() );
387 for(
const auto& node : nodes )
389 for(
size_t i=0; i<node.size(); ++i )
390 std::cout << node[i] <<
" ";
391 std::cout << std::endl;
394 const unsigned int nVx = dgf.elements[ 0 ].size();
396 typedef std::vector< int > face_t;
397 std::map< face_t, int > tmpFaces;
399 const int nFaces = (nVx == dim+1) ? dim+1 : 2*dim;
401 Dune::GeometryType type( (nVx == dim+1) ?
402 Impl :: SimplexTopology< dim > :: type :: id :
403 Impl :: CubeTopology < dim > :: type :: id,
406 const auto& refElem = Dune::referenceElement< typename Grid::ctype, dim >( type );
408 cells.resize( dgf.nofelements );
412 for(
int n = 0; n < dgf.nofelements; ++n )
414 const auto& elem = dgf.elements[ n ];
415 auto& cell = cells[ n ];
416 assert( elem.size() == nVx );
417 cell.resize( nFaces );
418 for(
int f=0; f<nFaces; ++f )
420 const int nFaceVx = refElem.size(f, 1, dim);
421 face.resize( nFaceVx );
422 for(
int j=0; j<nFaceVx; ++j )
424 face[ j ] = elem[ refElem.subEntity(f, 1, j , dim) ];
426 std::sort( face.begin(), face.end() );
427 auto it = tmpFaces.find( face );
429 if( it == tmpFaces.end() )
432 tmpFaces.insert( std::make_pair( face, myFaceNo ) );
435 myFaceNo = it->second;
437 cell[ f ] = myFaceNo;
451 mutable std::unique_ptr< Grid > gridPtr_;
454 std::vector< std::vector< double > > vtxParams_;
462 template<
int dim,
int dimworld >
465 static int refineStepsForHalf ()
470 static double refineWeight ()
identical grid wrapper
Definition: grid.hh:158
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Definition: dgfparser.hh:46
Definition: dgfparser.hh:81