My Project
dgfparser.hh
1// -*- mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=2 sw=2 sts=2:
3#ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH
4#define DUNE_POLYHEDRALGRID_DGFPARSER_HH
5
6#include <algorithm>
7#include <numeric>
8#include <memory>
9#include <utility>
10
11#include <dune/common/typetraits.hh>
12#include <dune/common/version.hh>
13
14#include <dune/geometry/referenceelements.hh>
15
16#include <dune/grid/common/gridfactory.hh>
17#include <dune/grid/io/file/dgfparser/dgfparser.hh>
18
19#if DUNE_VERSION_NEWER(DUNE_GRID,2,7)
20#include <dune/grid/io/file/dgfparser/blocks/polyhedron.hh>
21#endif
22
23#include <opm/grid/polyhedralgrid/gridfactory.hh>
24
25#if HAVE_ECL_INPUT
26#include <opm/input/eclipse/Parser/ParseContext.hpp>
27#include <opm/input/eclipse/Parser/Parser.hpp>
28#include <opm/input/eclipse/Deck/Deck.hpp>
29#endif
30
31namespace Dune
32{
33
34 namespace dgf
35 {
36
37#if ! DUNE_VERSION_NEWER(DUNE_GRID,2,7)
38 namespace PolyhedralGrid
39 {
40
41 // PolygonBlock
42 // ------------
43
45 : public BasicBlock
46 {
47 PolygonBlock ( std::istream &in, int numVtx, int vtxOfs )
48 : BasicBlock( in, "Polygon" ), vtxBegin_( vtxOfs ), vtxEnd_( vtxOfs + numVtx )
49 {}
50
51 int get ( std::vector< std::vector< int > > &polygons )
52 {
53 reset();
54 std::vector< int > polygon;
55 while( getnextline() )
56 {
57 polygon.clear();
58 for( int vtxIdx; getnextentry( vtxIdx ); )
59 {
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_ );
63 }
64
65 polygons.push_back( polygon );
66 }
67 return polygons.size();
68 }
69
70 private:
71 int vtxBegin_, vtxEnd_;
72 };
73
74
75
76 // PolyhedronBlock
77 // ---------------
78
80 : public BasicBlock
81 {
82 explicit PolyhedronBlock ( std::istream &in, int numPolys )
83 : BasicBlock( in, "Polyhedron" ), numPolys_( numPolys )
84 {}
85
86 int get ( std::vector< std::vector< int > > &polyhedra )
87 {
88 reset();
89 std::vector< int > polyhedron;
90 int minPolyId = 1;
91 while( getnextline() )
92 {
93 polyhedron.clear();
94 for( int polyIdx; getnextentry( polyIdx ); )
95 {
96 if( (polyIdx < 0) || (polyIdx > numPolys_) )
97 DUNE_THROW( DGFException, "Error in " << *this << ": Invalid polygon index (" << polyIdx << " not int [0, " << numPolys_ << "])" );
98
99 minPolyId = std::min( minPolyId, polyIdx );
100 polyhedron.push_back( polyIdx );
101 }
102
103 polyhedra.push_back( polyhedron );
104 }
105
106 // substract minimal number to have 0 starting numbering
107 if( minPolyId > 0 )
108 {
109 const size_t polySize = polyhedra.size();
110 for( size_t i=0; i<polySize; ++i )
111 {
112 const size_t pSize = polyhedra[ i ].size();
113 for( size_t j=0; j<pSize; ++j )
114 {
115 polyhedra[ i ][ j ] -= minPolyId;
116 }
117 }
118 }
119 return polyhedra.size();
120 }
121
122 private:
123 const int numPolys_;
124 };
125
126 } // namespace PolyhedralGrid
127
128 using PolyhedralGrid :: PolygonBlock;
129 using PolyhedralGrid :: PolyhedronBlock;
130#endif
131
132 } // namespace dgf
133
134
135
136 // DGFGridFactory for PolyhedralGrid
137 // ---------------------------------
138
139 template< int dim, int dimworld, class coord_t >
140 struct DGFGridFactory< PolyhedralGrid< dim, dimworld, coord_t > >
141 {
143
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;
148
149 explicit DGFGridFactory ( std::istream &input, MPICommunicator = MPIHelper::getCommunicator() )
150 : gridPtr_(),
151 grid_( nullptr )
152 {
153 input.clear();
154 input.seekg( 0 );
155 if( !input )
156 DUNE_THROW( DGFException, "Error resetting input stream" );
157 generate( input );
158 }
159
160 explicit DGFGridFactory ( const std::string &filename, MPICommunicator /* comm */ = MPIHelper::getCommunicator() )
161 : gridPtr_(),
162 grid_( nullptr )
163 {
164 std::ifstream input( filename );
165 if( !input )
166 DUNE_THROW( DGFException, "Macrofile '" << filename << "' not found" );
167
168#if HAVE_ECL_INPUT
169 if( !DuneGridFormatParser::isDuneGridFormat( input ) )
170 {
171 Opm::Parser parser;
172 const auto deck = parser.parseString( filename );
173 std::vector<double> porv;
174
175 gridPtr_.reset( new Grid( deck, porv ) );
176 return ;
177 }
178 else
179#endif
180 {
181 generate( input );
182 }
183 }
184
185 Grid *grid () const
186 {
187 if( ! grid_ )
188 {
189 // set pointer to grid and avoid grid being deleted
190 grid_ = gridPtr_.release();
191 }
192 return grid_;
193 }
194
195 template< class Intersection >
196 bool wasInserted ( const Intersection& /*intersection*/ ) const
197 {
198 return false;
199 }
200
201 template< class Intersection >
202 int boundaryId ( const Intersection& ) const
203 {
204 return false;
205 }
206
207 bool haveBoundaryParameters () const { return false; }
208
209 template< int codim >
210 int numParameters () const
211 {
212 //return (codim == dimension ? numVtxParams_ : 0);;
213 return 0;
214 }
215
216 template< class Intersection >
217 const typename DGFBoundaryParameter::type &
218 boundaryParameter ( const Intersection& ) const
219 {
220 return DGFBoundaryParameter::defaultValue();;
221 }
222
223 template< class Entity >
224 std::vector< double > &parameter ( const Entity& )
225 {
226 static std::vector< double > dummy;
227 return dummy;
228 }
229
230 private:
231 int readVertices ( std::istream &input, std::vector< std::vector< double > > &vertices )
232 {
233 int dimWorld = Grid::dimensionworld ;
234 dgf::VertexBlock vtxBlock( input, dimWorld );
235 if( !vtxBlock.isactive() )
236 DUNE_THROW( DGFException, "Vertex block not found" );
237
238 vtxBlock.get( vertices, vtxParams_, numVtxParams_ );
239 return vtxBlock.offset();
240 }
241
242 std::vector< std::vector< int > > readPolygons ( std::istream &input, int numVtx, int vtxOfs )
243 {
244 dgf::PolygonBlock polygonBlock( input, numVtx, vtxOfs );
245 if( !polygonBlock.isactive() )
246 DUNE_THROW( DGFException, "Polygon block not found" );
247
248 std::vector< std::vector< int > > polygons;
249 polygonBlock.get( polygons );
250 return polygons;
251 }
252
253 std::vector< std::vector< int > > readPolyhedra ( std::istream &input, int numPolygons )
254 {
255 dgf::PolyhedronBlock polyhedronBlock( input, numPolygons );
256 std::vector< std::vector< int > > polyhedra;
257 if( polyhedronBlock.isactive() )
258 {
259 polyhedronBlock.get( polyhedra );
260 }
261 return polyhedra;
262 }
263
264 template< class Iterator >
265 void copy ( Iterator begin, Iterator end, double *dest )
266 {
267 for( ; begin != end; ++begin )
268 dest = std::copy( begin->begin(), begin->end(), dest );
269 }
270
271 template< class Iterator >
272 void copy ( Iterator begin, Iterator end, int *dest, int *offset )
273 {
274 int size = 0;
275 for( ; begin != end; ++begin )
276 {
277 *(offset++) = size;
278 size += begin->size();
279 dest = std::copy( begin->begin(), begin->end(), dest );
280 }
281 *offset = size;
282 }
283
284 void generate ( std::istream &input )
285 {
286 // check whether an interval block is active, otherwise read polyhedrons
287 dgf::IntervalBlock intervalBlock( input );
288 if( intervalBlock.isactive() )
289 {
290 if( intervalBlock.numIntervals() != 1 )
291 DUNE_THROW( DGFException, "Currently, CpGrid can only handle 1 interval block." );
292
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 );
296
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 ];
300
301 gridPtr_.reset( new Grid( interval.n, spacing ) );
302 return ;
303 }
304 else // polyhedral input
305 {
306 typedef std::vector< std::vector< double > > CoordinateVectorType;
307 CoordinateVectorType nodes;
308
309 typedef std::vector< std::vector< int > > IndexVectorType;
310 IndexVectorType faces;
311 IndexVectorType cells;
312
313 const int vtxOfs = readVertices( input, nodes );
314
315 faces = readPolygons ( input, nodes.size(), vtxOfs );
316 cells = readPolyhedra( input, faces.size() );
317
318 if( cells.empty() )
319 {
320 DUNE_THROW( DGFException, "Polyhedron block not found!" );
321 }
322
323 typedef GridFactory< Grid > GridFactoryType;
324 typedef typename GridFactoryType :: Coordinate Coordinate ;
325
326 GridFactoryType gridFactory;
327
328 const int nNodes = nodes.size();
329 Coordinate node( 0 );
330 for( int i=0; i<nNodes; ++i )
331 {
332 for( int d=0; d<Coordinate::dimension; ++d )
333 node[ d ] = nodes[ i ][ d ];
334
335 gridFactory.insertVertex( node );
336 }
337 //nodes.swap( CoordinateVectorType() );
338
339 // insert faces with type none/dim-1
340 GeometryType type;
341 type = Dune::GeometryTypes::none(Grid::dimension-1);
342 std::vector< unsigned int > numbers;
343
344 const int nFaces = faces.size();
345 for(int i = 0; i < nFaces; ++ i )
346 {
347 // copy values into appropriate data type
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 );
352 }
353
354 // free memory
355 //faces.swap( IndexVectorType() );
356
357 // insert cells with type none/dim
358 type = Dune::GeometryTypes::none(Grid::dimension);
359
360 const int nCells = cells.size();
361 for(int i = 0; i < nCells; ++ i )
362 {
363 // copy values into appropriate data type
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 );
368 }
369 //cells.swap( IndexVectorType() );
370
371#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7)
372 gridPtr_ = gridFactory.createGrid();
373#else
374 gridPtr_.reset( gridFactory.createGrid() );
375#endif
376
377 } // end else branch
378
379 // alternative conversion to polyhedral format that does not work yet.
380#if 0
381 else // convert dgf input to polyhedral
382 {
383 nodes.resize( dgf.nofvtx );
384 // copy vertices
385 std::copy( dgf.vtx.begin(), dgf.vtx.end(), nodes.begin() );
386
387 for( const auto& node : nodes )
388 {
389 for( size_t i=0; i<node.size(); ++i )
390 std::cout << node[i] << " ";
391 std::cout << std::endl;
392 }
393
394 const unsigned int nVx = dgf.elements[ 0 ].size();
395
396 typedef std::vector< int > face_t;
397 std::map< face_t, int > tmpFaces;
398
399 const int nFaces = (nVx == dim+1) ? dim+1 : 2*dim;
400
401 Dune::GeometryType type( (nVx == dim+1) ?
402 Impl :: SimplexTopology< dim > :: type :: id :
403 Impl :: CubeTopology < dim > :: type :: id,
404 dim );
405
406 const auto& refElem = Dune::referenceElement< typename Grid::ctype, dim >( type );
407
408 cells.resize( dgf.nofelements );
409
410 face_t face;
411 int faceNo = 0;
412 for( int n = 0; n < dgf.nofelements; ++n )
413 {
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 )
419 {
420 const int nFaceVx = refElem.size(f, 1, dim);
421 face.resize( nFaceVx );
422 for( int j=0; j<nFaceVx; ++j )
423 {
424 face[ j ] = elem[ refElem.subEntity(f, 1, j , dim) ];
425 }
426 std::sort( face.begin(), face.end() );
427 auto it = tmpFaces.find( face );
428 int myFaceNo = -1;
429 if( it == tmpFaces.end() )
430 {
431 myFaceNo = faceNo++;
432 tmpFaces.insert( std::make_pair( face, myFaceNo ) );
433 }
434 else
435 myFaceNo = it->second;
436
437 cell[ f ] = myFaceNo;
438 }
439
440 /*
441 for( const auto& c : cell )
442 {
443 std::cout << c << " ";
444 }
445 std::cout << std::endl;
446 */
447 }
448#endif
449 }
450
451 mutable std::unique_ptr< Grid > gridPtr_;
452 mutable Grid* grid_;
453 int numVtxParams_;
454 std::vector< std::vector< double > > vtxParams_;
455 };
456
457
458
459 // DGFGridInfo for PolyhedralGrid
460 // ------------------------------
461
462 template< int dim, int dimworld >
463 struct DGFGridInfo< PolyhedralGrid< dim, dimworld > >
464 {
465 static int refineStepsForHalf ()
466 {
467 return 0;
468 }
469
470 static double refineWeight ()
471 {
472 return 0;
473 }
474 };
475
476} // namespace Dune
477
478#endif // #ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH
identical grid wrapper
Definition: grid.hh:163
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Definition: dgfparser.hh:46