156 :
public GridDefaultImplementation
157 < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > >
162 typedef GridDefaultImplementation
173 destroy_grid( grdPtr );
178 typedef std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > UnstructuredGridPtr;
180 static UnstructuredGridPtr
181 allocateGrid ( std::size_t nCells, std::size_t nFaces, std::size_t nFaceNodes, std::size_t nCellFaces, std::size_t nNodes )
184 UnstructuredGridType *grid = allocate_grid( dimworld, nCells, nFaces, nFaceNodes, nCellFaces, nNodes );
186 DUNE_THROW( GridError,
"Unable to allocate grid" );
187 return UnstructuredGridPtr( grid );
191 computeGeometry ( UnstructuredGridPtr& ug )
197 compute_geometry( ugPtr );
201 typedef PolyhedralGridFamily< dim, dimworld, coord_t > GridFamily;
208 typedef typename GridFamily::Traits
Traits;
216 template<
int codim >
237 template< PartitionIteratorType pitype >
240 typedef typename GridFamily::Traits::template Partition< pitype >::LevelGridView LevelGridView;
241 typedef typename GridFamily::Traits::template Partition< pitype >::LeafGridView LeafGridView;
246 typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
308 typedef typename Traits::ctype
ctype;
314 typedef typename Traits :: GlobalCoordinate GlobalCoordinate;
328 const std::vector<double>& poreVolumes = std::vector<double> ())
329 : gridPtr_( createGrid( inputGrid, poreVolumes ) ),
331 comm_( MPIHelper::getCommunicator() ),
332 leafIndexSet_( *this ),
333 globalIdSet_( *this ),
334 localIdSet_( *this ),
347 const std::vector< double >& dx )
348 : gridPtr_( createGrid( n, dx ) ),
350 comm_( MPIHelper::getCommunicator()),
351 leafIndexSet_( *this ),
352 globalIdSet_( *this ),
353 localIdSet_( *this ),
366 : gridPtr_( std::move( gridPtr ) ),
368 comm_( MPIHelper::getCommunicator() ),
369 leafIndexSet_( *this ),
370 globalIdSet_( *this ),
371 localIdSet_( *this ),
387 comm_( MPIHelper::getCommunicator() ),
388 leafIndexSet_( *this ),
389 globalIdSet_( *this ),
390 localIdSet_( *this ),
400 operator const UnstructuredGridType& ()
const {
return grid_; }
427 int size (
int ,
int codim )
const
429 return size( codim );
444 else if ( codim == 1 )
448 else if ( codim == dim )
454 std::cerr <<
"Warning: codimension " << codim <<
" not available in PolyhedralGrid" << std::endl;
467 int size (
int , GeometryType type )
const
469 return size( dim - type.dim() );
476 int size ( GeometryType type )
const
478 return size( dim - type.dim() );
489 return nBndSegments_;
493 template<
int codim >
496 return leafbegin< codim, All_Partition >();
499 template<
int codim >
502 return leafend< codim, All_Partition >();
505 template<
int codim, PartitionIteratorType pitype >
506 typename Codim< codim >::template Partition< pitype >::LeafIterator
509 typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
510 return Impl( extraData(),
true );
513 template<
int codim, PartitionIteratorType pitype >
514 typename Codim< codim >::template Partition< pitype >::LeafIterator
517 typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
518 return Impl( extraData(),
false );
521 template<
int codim >
524 return leafbegin< codim, All_Partition >();
527 template<
int codim >
530 return leafend< codim, All_Partition >();
533 template<
int codim, PartitionIteratorType pitype >
534 typename Codim< codim >::template Partition< pitype >::LevelIterator
535 lbegin (
const int )
const
537 return leafbegin< codim, pitype > ();
540 template<
int codim, PartitionIteratorType pitype >
541 typename Codim< codim >::template Partition< pitype >::LevelIterator
542 lend (
const int )
const
544 return leafend< codim, pitype > ();
559 return leafIndexSet();
564 return leafIndexSet_;
567 void globalRefine (
int )
597 template<
class DataHandle >
626 return (codim == 0 ) ? 1 : 0;
662 template<
class DataHandle>
665 CommunicationDirection ,
668 OPM_THROW(std::runtime_error,
"communicate not implemented for polyhedreal grid!");
683 template<
class DataHandle>
686 CommunicationDirection )
const
688 OPM_THROW(std::runtime_error,
"communicate not implemented for polyhedreal grid!");
694 OPM_THROW(std::runtime_error,
"switch to global view not implemented for polyhedreal grid!");
700 OPM_THROW(std::runtime_error,
"switch to distributed view not implemented for polyhedreal grid!");
711 const CommunicationType &
comm ()
const
747 template<
class DataHandle,
class Data >
767 template<
class DofManager >
774 template< PartitionIteratorType pitype >
777 typedef typename Partition< pitype >::LevelGridView View;
778 typedef typename View::GridViewImp ViewImp;
779 return View( ViewImp( *
this ) );
783 template< PartitionIteratorType pitype >
786 typedef typename Traits::template Partition< pitype >::LeafGridView View;
787 typedef typename View::GridViewImp ViewImp;
788 return View( ViewImp( *
this ) );
794 typedef typename LevelGridView::GridViewImp ViewImp;
801 typedef typename LeafGridView::GridViewImp ViewImp;
802 return LeafGridView( ViewImp( *
this ) );
806 template<
class EntitySeed >
813 return EntityPointer( EntityPointerImpl( EntityImpl( extraData(), seed ) ) );
817 template<
class EntitySeed >
822 return EntityImpl( extraData(), seed );
844 const std::array<int, 3>& logicalCartesianSize()
const
849 const int* globalCell()
const
855 const int* globalCellPtr()
const
860 void getIJK(
const int c, std::array<int,3>& ijk)
const
862 int gc = globalCell()[c];
863 ijk[0] = gc % logicalCartesianSize()[0]; gc /= logicalCartesianSize()[0];
864 ijk[1] = gc % logicalCartesianSize()[1];
865 ijk[2] = gc / logicalCartesianSize()[1];
882 template<
class DataHandle>
885 OPM_THROW(std::runtime_error,
"ScatterData not implemented for polyhedral grid!");
890 UnstructuredGridType* createGrid(
const Opm::EclipseGrid& inputGrid,
const std::vector< double >& poreVolumes )
const
894 g.
dims[0] = inputGrid.getNX();
895 g.dims[1] = inputGrid.getNY();
896 g.dims[2] = inputGrid.getNZ();
898 std::vector<double> coord = inputGrid.getCOORD( );
899 std::vector<double> zcorn = inputGrid.getZCORN( );
900 std::vector<int> actnum = inputGrid.getACTNUM( );
902 g.coord = coord.data();
903 g.zcorn = zcorn.data();
904 g.actnum = actnum.data();
906 const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0;
908 if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != Opm::MinpvMode::Inactive))
911 const std::vector<double>& minpvv = inputGrid.getMinpvVector();
916 const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
917 std::vector<double> thickness(cartGridSize);
918 for (
size_t i = 0; i < cartGridSize; ++i) {
919 thickness[i] = inputGrid.getCellThickness(i);
921 mp.process(thickness, z_tolerance, inputGrid.getPinchMaxEmptyGap() , poreVolumes,
922 minpvv, actnum, opmfil, zcorn.data());
936 UnstructuredGridType* cgrid = create_grid_cornerpoint(&g, z_tolerance);
938 OPM_THROW(std::runtime_error,
"Failed to construct grid.");
944 UnstructuredGridType* createGrid(
const std::vector< int >& n,
const std::vector< double >& dx )
const
946 UnstructuredGridType* cgrid = nullptr ;
947 assert(
int(n.size()) == dim );
950 cgrid = create_grid_cart2d( n[ 0 ], n[ 1 ], dx[ 0 ], dx[ 1 ] );
954 cgrid = create_grid_hexa3d( n[ 0 ], n[ 1 ], n[ 2 ], dx[ 0 ], dx[ 1 ], dx[ 2 ] );
959 OPM_THROW(std::runtime_error,
"Failed to construct grid.");
965 typedef typename Traits :: ExtraData ExtraData;
966 ExtraData extraData ()
const {
return this; }
968 template <
class EntitySeed>
969 int corners(
const EntitySeed& seed )
const
971 const int codim = EntitySeed :: codimension;
972 const int index = seed.index();
974 return cellVertices_[ index ].size();
982 template <
class EntitySeed>
984 corner (
const EntitySeed& seed,
const int i )
const
986 const int codim = EntitySeed :: codimension;
989 const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ];
997 const int crners = corners( seed );
998 const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i;
1000 return copyToGlobalCoordinate( grid_.
node_coordinates + GlobalCoordinate :: dimension * faceVertex );
1004 const int coordIndex = GlobalCoordinate :: dimension * seed.index();
1007 return GlobalCoordinate( 0 );
1010 template <
class EntitySeed>
1011 int subEntities(
const EntitySeed& seed,
const int codim )
const
1013 const int index = seed.index();
1014 if( seed.codimension == 0 )
1021 return cellVertices_[ index ].size();
1023 else if( seed.codimension == 1 )
1030 else if ( seed.codimension == dim )
1038 template <
int codim,
class EntitySeedArg >
1039 typename Codim<codim>::EntitySeed
1040 subEntitySeed(
const EntitySeedArg& baseSeed,
const int i )
const
1042 assert( codim >= EntitySeedArg::codimension );
1043 assert( i>= 0 && i<subEntities( baseSeed, codim ) );
1044 typedef typename Codim<codim>::EntitySeed EntitySeed;
1047 if( codim == EntitySeedArg::codimension )
1049 return EntitySeed( baseSeed.index() );
1052 if( EntitySeedArg::codimension == 0 )
1058 else if ( codim == dim )
1060 return EntitySeed( cellVertices_[ baseSeed.index() ][ i ] );
1063 else if ( EntitySeedArg::codimension == 1 && codim == dim )
1068 DUNE_THROW(NotImplemented,
"codimension not available");
1069 return EntitySeed();
1072 template <
int codim>
1073 typename Codim<codim>::EntitySeed
1074 subEntitySeed(
const typename Codim<1>::EntitySeed& faceSeed,
const int i )
const
1076 assert( i>= 0 && i<subEntities( faceSeed, codim ) );
1077 typedef typename Codim<codim>::EntitySeed EntitySeed;
1080 return EntitySeed( faceSeed.index() );
1082 else if ( codim == dim )
1088 DUNE_THROW(NotImplemented,
"codimension not available");
1092 bool hasBoundaryIntersections(
const typename Codim<0>::EntitySeed& seed )
const
1094 const int faces = subEntities( seed, 1 );
1095 for(
int f=0; f<faces; ++f )
1097 const auto faceSeed = this->
template subEntitySeed<1>( seed, f );
1098 if( isBoundaryFace( faceSeed ) )
1104 bool isBoundaryFace(
const int face )
const
1107 const int facePos = 2 * face;
1111 bool isBoundaryFace(
const typename Codim<1>::EntitySeed& faceSeed )
const
1113 assert( faceSeed.isValid() );
1114 return isBoundaryFace( faceSeed.index() );
1117 int boundarySegmentIndex(
const typename Codim<0>::EntitySeed& seed,
const int face )
const
1119 const auto faceSeed = this->
template subEntitySeed<1>( seed, face );
1120 assert( faceSeed.isValid() );
1121 const int facePos = 2 * faceSeed.index();
1128 const std::vector< GeometryType > &geomTypes (
const unsigned int codim )
const
1130 static std::vector< GeometryType > emptyDummy;
1131 if (codim < geomTypes_.size())
1133 return geomTypes_[codim];
1139 template <
class Seed >
1140 GeometryType geometryType(
const Seed& seed )
const
1142 if( Seed::codimension == 0 )
1144 assert(!geomTypes( Seed::codimension ).empty());
1145 return geomTypes( Seed::codimension )[ 0 ];
1150 if( dim == 3 && Seed::codimension == 1 )
1153 const int nVx = corners( seed );
1155 face = Dune::GeometryTypes::cube(2);
1157 face = Dune::GeometryTypes::simplex(2);
1159 face = Dune::GeometryTypes::none(2);
1164 assert(!geomTypes( Seed::codimension ).empty());
1165 return geomTypes( Seed::codimension )[ 0 ];
1169 int indexInInside(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1171 return ( grid_.
cell_facetag ) ? cartesianIndexInInside( seed, i ) : i;
1174 int cartesianIndexInInside(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1176 assert( i>= 0 && i<subEntities( seed, 1 ) );
1180 typename Codim<0>::EntitySeed
1181 neighbor(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1183 const int face = this->
template subEntitySeed<1>( seed, i ).index();
1185 if( nb == seed.index() )
1190 typedef typename Codim<0>::EntitySeed EntitySeed;
1191 return EntitySeed( nb );
1195 indexInOutside(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1200 const int in_inside = cartesianIndexInInside( seed, i );
1201 return in_inside + ((in_inside % 2) ? -1 : 1);
1205 typedef typename Codim<0>::EntitySeed EntitySeed;
1206 EntitySeed nb = neighbor( seed, i );
1207 const int faces = subEntities( seed, 1 );
1208 for(
int face = 0; face<faces; ++ face )
1210 if( neighbor( nb, face ).equals(seed) )
1212 return indexInInside( nb, face );
1215 DUNE_THROW(InvalidStateException,
"inverse intersection not found");
1220 template <
class EntitySeed>
1222 outerNormal(
const EntitySeed& seed,
const int i )
const
1224 const int face = this->
template subEntitySeed<1>( seed, i ).index();
1225 const int normalIdx = face * GlobalCoordinate :: dimension ;
1226 GlobalCoordinate normal = copyToGlobalCoordinate( grid_.
face_normals + normalIdx );
1228 if( nb != seed.index() )
1235 template <
class EntitySeed>
1237 unitOuterNormal(
const EntitySeed& seed,
const int i )
const
1239 const int face = this->
template subEntitySeed<1>( seed, i ).index();
1240 if( seed.index() == grid_.
face_cells[ 2*face ] )
1242 return unitOuterNormals_[ face ];
1246 GlobalCoordinate normal = unitOuterNormals_[ face ];
1252 template <
class EntitySeed>
1253 GlobalCoordinate centroids(
const EntitySeed& seed )
const
1255 if( ! seed.isValid() )
1256 return GlobalCoordinate( 0 );
1258 const int index = GlobalCoordinate :: dimension * seed.index();
1259 const int codim = EntitySeed::codimension;
1260 assert( index >= 0 && index <
size( codim ) * GlobalCoordinate :: dimension );
1266 else if ( codim == 1 )
1270 else if( codim == dim )
1276 DUNE_THROW(InvalidStateException,
"codimension not implemented");
1277 return GlobalCoordinate( 0 );
1281 GlobalCoordinate copyToGlobalCoordinate(
const double* coords )
const
1283 GlobalCoordinate coordinate;
1284 for(
int i=0; i<GlobalCoordinate::dimension; ++i )
1286 coordinate[ i ] = coords[ i ];
1291 template <
class EntitySeed>
1292 double volumes(
const EntitySeed& seed )
const
1294 static const int codim = EntitySeed::codimension;
1295 if( codim == dim || ! seed.isValid() )
1301 assert( seed.isValid() );
1307 else if ( codim == 1 )
1313 DUNE_THROW(InvalidStateException,
"codimension not implemented");
1323 for(
int i=0; i<3; ++i )
1325 cartDims_[ i ] = grid_.
cartdims[ i ];
1329 const int numCells =
size( 0 );
1331 cellVertices_.resize( numCells );
1336 typedef std::array<int, 3> KeyType;
1337 std::map< const KeyType, const int > vertexFaceTags;
1338 const int vertexFacePattern [8][3] = {
1349 for(
int i=0; i<8; ++i )
1351 KeyType key; key.fill( 4 );
1352 for(
int j=0; j<dim; ++j )
1354 key[ j ] = vertexFacePattern[ i ][ j ];
1357 vertexFaceTags.insert( std::make_pair( key, i ) );
1360 for (
int c = 0; c < numCells; ++c)
1370 typedef std::map<int,int> vertexmap_t;
1371 typedef typename vertexmap_t :: iterator iterator;
1373 std::vector< vertexmap_t > cell_pts( dim*2 );
1382 const int node = grid_.
face_nodes[ nodepos ];
1383 iterator it = cell_pts[ faceTag ].find( node );
1384 if( it == cell_pts[ faceTag ].end() )
1386 cell_pts[ faceTag ].insert( std::make_pair( node, 1 ) );
1396 typedef std::map< int, std::set<int> > vertexlist_t;
1397 vertexlist_t vertexList;
1399 for(
int faceTag = 0; faceTag<dim*2; ++faceTag )
1401 for( iterator it = cell_pts[ faceTag ].begin(),
1402 end = cell_pts[ faceTag ].end(); it != end; ++it )
1406 if( (*it).second == 1 )
1408 vertexList[ (*it).first ].insert( faceTag );
1413 assert(
int(vertexList.size()) == ( dim == 2 ? 4 : 8) );
1415 cellVertices_[ c ].resize( vertexList.size() );
1416 for(
auto it = vertexList.begin(), end = vertexList.end(); it != end; ++it )
1418 assert( (*it).second.size() == dim );
1419 KeyType key; key.fill( 4 );
1421 std::copy( (*it).second.begin(), (*it).second.end(), key.begin() );
1422 auto vx = vertexFaceTags.find( key );
1423 assert( vx != vertexFaceTags.end() );
1424 if( vx != vertexFaceTags.end() )
1426 if( (*vx).second >=
int(cellVertices_[ c ].size()) )
1427 cellVertices_[ c ].resize( (*vx).second+1 );
1429 cellVertices_[ c ][ (*vx).second ] = (*it).first ;
1435 geomTypes_.resize(dim + 1);
1437 for (
int codim = 0; codim <= dim; ++codim)
1439 tmp = Dune::GeometryTypes::cube(dim - codim);
1440 geomTypes_[codim].push_back(tmp);
1446 int minVx = std::numeric_limits<int>::max();
1448 for (
int c = 0; c < numCells; ++c)
1450 std::set<int> cell_pts;
1456 cell_pts.insert(fnbeg, fnend);
1459 cellVertices_[ c ].resize( cell_pts.size() );
1460 std::copy(cell_pts.begin(), cell_pts.end(), cellVertices_[ c ].begin() );
1461 maxVx = std::max( maxVx,
int( cell_pts.size() ) );
1462 minVx = std::min( minVx,
int( cell_pts.size() ) );
1465 if( minVx == maxVx && maxVx == 4 )
1467 for (
int c = 0; c < numCells; ++c)
1469 assert( cellVertices_[ c ].
size() == 4 );
1470 GlobalCoordinate center( 0 );
1471 GlobalCoordinate p[ dim+1 ];
1472 for(
int i=0; i<dim+1; ++i )
1474 const int vertex = cellVertices_[ c ][ i ];
1476 for(
int d=0; d<dim; ++d )
1483 for(
int d=0; d<dim; ++d )
1488 Dune::GeometryType simplex;
1489 simplex = Dune::GeometryTypes::simplex(dim);
1491 typedef Dune::AffineGeometry< ctype, dim, dimworld> AffineGeometryType;
1492 AffineGeometryType geometry( simplex, p );
1500 for(
int face = 0 ; face < faces; ++face )
1503 const int b = grid_.
face_cells[ 2*face + 1 ];
1505 assert( a >=0 || b >=0 );
1510 GlobalCoordinate centerDiff( 0 );
1513 for(
int d=0; d<dimworld; ++d )
1520 for(
int d=0; d<dimworld; ++d )
1528 for(
int d=0; d<dimworld; ++d )
1535 for(
int d=0; d<dimworld; ++d )
1541 GlobalCoordinate normal( 0 );
1542 for(
int d=0; d<dimworld; ++d )
1547 if( centerDiff.two_norm() < 1e-10 )
1551 if( centerDiff * normal < 0 )
1559 bool allSimplex = true ;
1560 bool allCube = true ;
1562 for (
int c = 0; c < numCells; ++c)
1564 const int nVx = cellVertices_[ c ].size();
1576 geomTypes_.resize(dim + 1);
1578 for (
int codim = 0; codim <= dim; ++codim)
1582 tmp = Dune::GeometryTypes::simplex(dim - codim);
1583 geomTypes_[ codim ].push_back( tmp );
1587 tmp = Dune::GeometryTypes::cube(dim - codim);
1588 geomTypes_[ codim ].push_back( tmp );
1592 tmp = Dune::GeometryTypes::none(dim - codim);
1593 geomTypes_[ codim ].push_back( tmp );
1603 const int normalIdx = face * GlobalCoordinate :: dimension ;
1604 GlobalCoordinate normal = copyToGlobalCoordinate( grid_.
face_normals + normalIdx );
1605 normal /= normal.two_norm();
1606 unitOuterNormals_[ face ] = normal;
1608 if( isBoundaryFace( face ) )
1612 const int facePos = 2 * face ;
1617 grid_.
face_cells[ facePos ] = -nBndSegments_;
1619 else if ( grid_.
face_cells[ facePos+1 ] < 0 )
1621 grid_.
face_cells[ facePos+1 ] = -nBndSegments_;
1627 void print( std::ostream& out,
const UnstructuredGridType& grid )
const
1629 const int numCells = grid.number_of_cells;
1630 for(
int c=0; c<numCells; ++c )
1632 out <<
"cell " << c <<
" : faces = " << std::endl;
1633 for (
int hf=grid.cell_facepos[ c ]; hf < grid.cell_facepos[c+1]; ++hf)
1638 out << f <<
" vx = " ;
1639 while( fnbeg != fnend )
1641 out << *fnbeg <<
" ";
1648 const auto& vx = cellVertices_[ c ];
1649 out <<
"cell " << c <<
" : vertices = ";
1650 for(
size_t i=0; i<vx.size(); ++i )
1651 out << vx[ i ] <<
" ";
1658 UnstructuredGridPtr gridPtr_;
1659 const UnstructuredGridType& grid_;
1661 CommunicationType comm_;
1662 std::array< int, 3 > cartDims_;
1663 std::vector< std::vector< GeometryType > > geomTypes_;
1664 std::vector< std::vector< int > > cellVertices_;
1666 std::vector< GlobalCoordinate > unitOuterNormals_;
1672 size_t nBndSegments_;