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 
31 namespace 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 
44  struct PolygonBlock
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