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#include <dune/grid/io/file/dgfparser/blocks/polyhedron.hh>
24#include <opm/input/eclipse/Parser/ParseContext.hpp>
25#include <opm/input/eclipse/Parser/Parser.hpp>
26#include <opm/input/eclipse/Deck/Deck.hpp>
34 template<
int dim,
int dimworld,
class coord_t >
39 const static int dimension = Grid::dimension;
41 typedef typename Grid::template Codim<0>::Entity
Element;
42 typedef typename Grid::template Codim<dimension>::Entity
Vertex;
51 DUNE_THROW( DGFException,
"Error resetting input stream" );
59 std::ifstream input( filename );
61 DUNE_THROW( DGFException,
"Macrofile '" << filename <<
"' not found" );
64 if( !DuneGridFormatParser::isDuneGridFormat( input ) )
67 const auto deck = parser.parseString( filename );
68 std::vector<double> porv;
70 gridPtr_.reset(
new Grid( deck, porv ) );
85 grid_ = gridPtr_.release();
90 template<
class Intersection >
96 template<
class Intersection >
104 template<
int codim >
111 template<
class Intersection >
112 const typename DGFBoundaryParameter::type &
115 return DGFBoundaryParameter::defaultValue();;
118 template<
class Entity >
121 static std::vector< double > dummy;
126 int readVertices ( std::istream &input, std::vector< std::vector< double > > &vertices )
128 int dimWorld = Grid::dimensionworld ;
129 dgf::VertexBlock vtxBlock( input, dimWorld );
130 if( !vtxBlock.isactive() )
131 DUNE_THROW( DGFException,
"Vertex block not found" );
133 vtxBlock.get( vertices, vtxParams_, numVtxParams_ );
134 return vtxBlock.offset();
137 std::vector< std::vector< int > > readPolygons ( std::istream &input,
int numVtx,
int vtxOfs )
139 dgf::PolygonBlock polygonBlock( input, numVtx, vtxOfs );
140 if( !polygonBlock.isactive() )
141 DUNE_THROW( DGFException,
"Polygon block not found" );
143 std::vector< std::vector< int > > polygons;
144 polygonBlock.get( polygons );
148 std::vector< std::vector< int > > readPolyhedra ( std::istream &input,
int numPolygons )
150 dgf::PolyhedronBlock polyhedronBlock( input, numPolygons );
151 std::vector< std::vector< int > > polyhedra;
152 if( polyhedronBlock.isactive() )
154 polyhedronBlock.get( polyhedra );
159 template<
class Iterator >
160 void copy ( Iterator begin, Iterator end,
double *dest )
162 for( ; begin != end; ++begin )
163 dest = std::copy( begin->begin(), begin->end(), dest );
166 template<
class Iterator >
167 void copy ( Iterator begin, Iterator end,
int *dest,
int *offset )
170 for( ; begin != end; ++begin )
173 size += begin->size();
174 dest = std::copy( begin->begin(), begin->end(), dest );
179 void generate ( std::istream &input )
182 dgf::IntervalBlock intervalBlock( input );
183 if( intervalBlock.isactive() )
185 if( intervalBlock.numIntervals() != 1 )
186 DUNE_THROW( DGFException,
"Currently, CpGrid can only handle 1 interval block." );
188 if( intervalBlock.dimw() != dimworld )
189 DUNE_THROW( DGFException,
"CpGrid cannot handle an interval of dimension " << intervalBlock.dimw() <<
"." );
190 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
192 std::vector< double > spacing( dimworld );
193 for(
int i=0; i<dimworld; ++i )
194 spacing[ i ] = (interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ]) / interval.n[ i ];
196 gridPtr_.reset(
new Grid( interval.n, spacing ) );
201 typedef std::vector< std::vector< double > > CoordinateVectorType;
202 CoordinateVectorType nodes;
204 typedef std::vector< std::vector< int > > IndexVectorType;
205 IndexVectorType faces;
206 IndexVectorType cells;
208 const int vtxOfs = readVertices( input, nodes );
210 faces = readPolygons ( input, nodes.size(), vtxOfs );
211 cells = readPolyhedra( input, faces.size() );
215 DUNE_THROW( DGFException,
"Polyhedron block not found!" );
218 typedef GridFactory< Grid > GridFactoryType;
219 typedef typename GridFactoryType :: Coordinate Coordinate ;
221 GridFactoryType gridFactory;
223 const int nNodes = nodes.size();
224 Coordinate node( 0 );
225 for(
int i=0; i<nNodes; ++i )
227 for(
int d=0; d<Coordinate::dimension; ++d )
228 node[ d ] = nodes[ i ][ d ];
230 gridFactory.insertVertex( node );
236 type = Dune::GeometryTypes::none(Grid::dimension-1);
237 std::vector< unsigned int > numbers;
239 const int nFaces = faces.size();
240 for(
int i = 0; i < nFaces; ++ i )
243 std::vector<int>& face = faces[ i ];
244 numbers.resize( face.size() );
245 std::copy( face.begin(), face.end(), numbers.begin() );
246 gridFactory.insertElement( type, numbers );
253 type = Dune::GeometryTypes::none(Grid::dimension);
255 const int nCells = cells.size();
256 for(
int i = 0; i < nCells; ++ i )
259 std::vector<int>& cell = cells[ i ];
260 numbers.resize( cell.size() );
261 std::copy( cell.begin(), cell.end(), numbers.begin() );
262 gridFactory.insertElement( type, numbers );
266 gridPtr_ = gridFactory.createGrid();
273 nodes.resize( dgf.nofvtx );
275 std::copy( dgf.vtx.begin(), dgf.vtx.end(), nodes.begin() );
277 for(
const auto& node : nodes )
279 for(
size_t i=0; i<node.size(); ++i )
280 std::cout << node[i] <<
" ";
281 std::cout << std::endl;
284 const unsigned int nVx = dgf.elements[ 0 ].size();
286 typedef std::vector< int > face_t;
287 std::map< face_t, int > tmpFaces;
289 const int nFaces = (nVx == dim+1) ? dim+1 : 2*dim;
291 Dune::GeometryType type( (nVx == dim+1) ?
292 Impl :: SimplexTopology< dim > :: type :: id :
293 Impl :: CubeTopology < dim > :: type :: id,
296 const auto& refElem = Dune::referenceElement< typename Grid::ctype, dim >( type );
298 cells.resize( dgf.nofelements );
302 for(
int n = 0; n < dgf.nofelements; ++n )
304 const auto& elem = dgf.elements[ n ];
305 auto& cell = cells[ n ];
306 assert( elem.size() == nVx );
307 cell.resize( nFaces );
308 for(
int f=0; f<nFaces; ++f )
310 const int nFaceVx = refElem.size(f, 1, dim);
311 face.resize( nFaceVx );
312 for(
int j=0; j<nFaceVx; ++j )
314 face[ j ] = elem[ refElem.subEntity(f, 1, j , dim) ];
316 std::sort( face.begin(), face.end() );
317 auto it = tmpFaces.find( face );
319 if( it == tmpFaces.end() )
322 tmpFaces.insert( std::make_pair( face, myFaceNo ) );
325 myFaceNo = it->second;
327 cell[ f ] = myFaceNo;
341 mutable std::unique_ptr< Grid > gridPtr_;
344 std::vector< std::vector< double > > vtxParams_;
352 template<
int dim,
int dimworld >
identical grid wrapper
Definition: grid.hh:159
The namespace Dune is the main namespace for all Dune code.
Definition: common/CartesianIndexMapper.hpp:10
int numParameters() const
Definition: polyhedralgrid/dgfparser.hh:105
Grid * grid() const
Definition: polyhedralgrid/dgfparser.hh:80
bool haveBoundaryParameters() const
Definition: polyhedralgrid/dgfparser.hh:102
Grid::template Codim< 0 >::Entity Element
Definition: polyhedralgrid/dgfparser.hh:41
PolyhedralGrid< dim, dimworld, coord_t > Grid
Definition: polyhedralgrid/dgfparser.hh:37
int boundaryId(const Intersection &) const
Definition: polyhedralgrid/dgfparser.hh:97
std::vector< double > & parameter(const Entity &)
Definition: polyhedralgrid/dgfparser.hh:119
bool wasInserted(const Intersection &) const
Definition: polyhedralgrid/dgfparser.hh:91
DGFGridFactory(const std::string &filename, MPICommunicator=MPIHelper::getCommunicator())
Definition: polyhedralgrid/dgfparser.hh:55
const DGFBoundaryParameter::type & boundaryParameter(const Intersection &) const
Definition: polyhedralgrid/dgfparser.hh:113
DGFGridFactory(std::istream &input, MPICommunicator=MPIHelper::getCommunicator())
Definition: polyhedralgrid/dgfparser.hh:44
Grid::template Codim< dimension >::Entity Vertex
Definition: polyhedralgrid/dgfparser.hh:42
MPIHelper::MPICommunicator MPICommunicator
Definition: polyhedralgrid/dgfparser.hh:40
static double refineWeight()
Definition: polyhedralgrid/dgfparser.hh:360
static int refineStepsForHalf()
Definition: polyhedralgrid/dgfparser.hh:355