3 #ifndef DUNE_POLYHEDRALGRID_GRID_HH
4 #define DUNE_POLYHEDRALGRID_GRID_HH
10 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
13 #include <dune/common/version.hh>
14 #include <dune/common/parallel/mpihelper.hh>
17 #include <dune/grid/common/grid.hh>
19 #if DUNE_VERSION_GTE(DUNE_COMMON, 2, 7)
20 #include <dune/common/parallel/communication.hh>
22 #include <dune/common/parallel/collectivecommunication.hh>
26 #include <opm/grid/polyhedralgrid/capabilities.hh>
27 #include <opm/grid/polyhedralgrid/declaration.hh>
28 #include <opm/grid/polyhedralgrid/entity.hh>
29 #include <opm/grid/polyhedralgrid/entityseed.hh>
30 #include <opm/grid/polyhedralgrid/geometry.hh>
31 #include <opm/grid/polyhedralgrid/gridview.hh>
32 #include <opm/grid/polyhedralgrid/idset.hh>
35 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
37 #include <opm/grid/utility/ErrorMacros.hpp>
42 #include <opm/grid/GridManager.hpp>
44 #include <opm/grid/MinpvProcessor.hpp>
53 template<
int dim,
int dimworld,
typename coord_t >
60 typedef coord_t ctype;
68 static const int dimension = dim;
69 static const int dimensionworld = dimworld;
71 typedef Dune::FieldVector< ctype, dimensionworld > GlobalCoordinate ;
78 typedef Dune::Intersection< const Grid, LeafIntersectionImpl > LeafIntersection;
79 typedef Dune::Intersection< const Grid, LevelIntersectionImpl > LevelIntersection;
81 typedef Dune::IntersectionIterator< const Grid, LeafIntersectionIteratorImpl, LeafIntersectionImpl > LeafIntersectionIterator;
82 typedef Dune::IntersectionIterator< const Grid, LevelIntersectionIteratorImpl, LevelIntersectionImpl > LevelIntersectionIterator;
85 typedef Dune::EntityIterator< 0, const Grid, HierarchicIteratorImpl > HierarchicIterator;
97 typedef Dune::Entity< codim, dimension, const Grid, PolyhedralGridEntity > Entity;
100 typedef Entity EntityPointer;
105 template< PartitionIteratorType pitype >
109 typedef Dune::EntityIterator< codim, const Grid, LeafIteratorImpl > LeafIterator;
111 typedef LeafIterator LevelIterator;
114 typedef typename Partition< All_Partition >::LeafIterator LeafIterator;
115 typedef typename Partition< All_Partition >::LevelIterator LevelIterator;
124 typedef Dune::MPIHelper::MPICommunicator MPICommunicator;
125 typedef Dune::CollectiveCommunication<MPICommunicator> CollectiveCommunication;
127 template< PartitionIteratorType pitype >
130 typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, ctype, pitype > > LeafGridView;
131 typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, ctype, pitype > > LevelGridView;
134 typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
135 typedef typename Partition< All_Partition >::LevelGridView LevelGridView;
152 template <
int dim,
int dimworld,
typename coord_t >
155 :
public GridDefaultImplementation
156 < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > >
161 typedef GridDefaultImplementation
177 typedef std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > UnstructuredGridPtr;
179 static UnstructuredGridPtr
180 allocateGrid ( std::size_t nCells, std::size_t nFaces, std::size_t nFaceNodes, std::size_t nCellFaces, std::size_t nNodes )
185 DUNE_THROW( GridError,
"Unable to allocate grid" );
186 return UnstructuredGridPtr( grid );
190 computeGeometry ( UnstructuredGridPtr& ug )
200 typedef PolyhedralGridFamily< dim, dimworld, coord_t > GridFamily;
207 typedef typename GridFamily::Traits
Traits;
215 template<
int codim >
236 template< PartitionIteratorType pitype >
239 typedef typename GridFamily::Traits::template Partition< pitype >::LevelGridView LevelGridView;
240 typedef typename GridFamily::Traits::template Partition< pitype >::LeafGridView LeafGridView;
245 typedef typename Partition< All_Partition >::LeafGridView LeafGridView;
307 typedef typename Traits::ctype
ctype;
312 typedef typename Traits :: GlobalCoordinate GlobalCoordinate;
326 const std::vector<double>& poreVolumes = std::vector<double> ())
327 : gridPtr_( createGrid( inputGrid, poreVolumes ) ),
329 comm_( MPIHelper::getCommunicator() ),
330 leafIndexSet_( *this ),
331 globalIdSet_( *this ),
332 localIdSet_( *this ),
345 const std::vector< double >& dx )
346 : gridPtr_( createGrid( n, dx ) ),
348 comm_( MPIHelper::getCommunicator()),
349 leafIndexSet_( *this ),
350 globalIdSet_( *this ),
351 localIdSet_( *this ),
364 : gridPtr_( std::move( gridPtr ) ),
366 comm_( MPIHelper::getCommunicator() ),
367 leafIndexSet_( *this ),
368 globalIdSet_( *this ),
369 localIdSet_( *this ),
385 comm_( MPIHelper::getCommunicator() ),
386 leafIndexSet_( *this ),
387 globalIdSet_( *this ),
388 localIdSet_( *this ),
398 operator const UnstructuredGridType& ()
const {
return grid_; }
425 int size (
int ,
int codim )
const
427 return size( codim );
442 else if ( codim == 1 )
446 else if ( codim == dim )
452 std::cerr <<
"Warning: codimension " << codim <<
" not available in PolyhedralGrid" << std::endl;
465 int size (
int , GeometryType type )
const
467 return size( dim - type.dim() );
474 int size ( GeometryType type )
const
476 return size( dim - type.dim() );
487 return nBndSegments_;
491 template<
int codim >
494 return leafbegin< codim, All_Partition >();
497 template<
int codim >
500 return leafend< codim, All_Partition >();
503 template<
int codim, PartitionIteratorType pitype >
504 typename Codim< codim >::template Partition< pitype >::LeafIterator
507 typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
508 return Impl( extraData(),
true );
511 template<
int codim, PartitionIteratorType pitype >
512 typename Codim< codim >::template Partition< pitype >::LeafIterator
515 typedef typename Traits::template Codim< codim >::template Partition< pitype >::LeafIteratorImpl Impl;
516 return Impl( extraData(),
false );
519 template<
int codim >
522 return leafbegin< codim, All_Partition >();
525 template<
int codim >
528 return leafend< codim, All_Partition >();
531 template<
int codim, PartitionIteratorType pitype >
532 typename Codim< codim >::template Partition< pitype >::LevelIterator
533 lbegin (
const int )
const
535 return leafbegin< codim, pitype > ();
538 template<
int codim, PartitionIteratorType pitype >
539 typename Codim< codim >::template Partition< pitype >::LevelIterator
540 lend (
const int )
const
542 return leafend< codim, pitype > ();
557 return leafIndexSet();
562 return leafIndexSet_;
565 void globalRefine (
int )
595 template<
class DataHandle >
624 return (codim == 0 ) ? 1 : 0;
660 template<
class DataHandle>
663 CommunicationDirection ,
666 OPM_THROW(std::runtime_error,
"communicate not implemented for polyhedreal grid!");
681 template<
class DataHandle>
684 CommunicationDirection )
const
686 OPM_THROW(std::runtime_error,
"communicate not implemented for polyhedreal grid!");
692 OPM_THROW(std::runtime_error,
"switch to global view not implemented for polyhedreal grid!");
698 OPM_THROW(std::runtime_error,
"switch to distributed view not implemented for polyhedreal grid!");
745 template<
class DataHandle,
class Data >
765 template<
class DofManager >
772 template< PartitionIteratorType pitype >
775 typedef typename Partition< pitype >::LevelGridView View;
776 typedef typename View::GridViewImp ViewImp;
777 return View( ViewImp( *
this ) );
781 template< PartitionIteratorType pitype >
784 typedef typename Traits::template Partition< pitype >::LeafGridView View;
785 typedef typename View::GridViewImp ViewImp;
786 return View( ViewImp( *
this ) );
792 typedef typename LevelGridView::GridViewImp ViewImp;
799 typedef typename LeafGridView::GridViewImp ViewImp;
800 return LeafGridView( ViewImp( *
this ) );
804 template<
class EntitySeed >
811 return EntityPointer( EntityPointerImpl( EntityImpl( extraData(), seed ) ) );
815 template<
class EntitySeed >
820 return EntityImpl( extraData(), seed );
842 const std::array<int, 3>& logicalCartesianSize()
const
847 const int* globalCell()
const
853 const int* globalCellPtr()
const
858 void getIJK(
const int c, std::array<int,3>& ijk)
const
860 int gc = globalCell()[c];
861 ijk[0] = gc % logicalCartesianSize()[0]; gc /= logicalCartesianSize()[0];
862 ijk[1] = gc % logicalCartesianSize()[1];
863 ijk[2] = gc / logicalCartesianSize()[1];
871 template<
class DataHandle>
883 OPM_THROW(std::runtime_error,
"ScatterData not implemented for polyhedral grid!");
888 UnstructuredGridType* createGrid(
const Opm::EclipseGrid& inputGrid,
const std::vector< double >& poreVolumes )
const
892 g.
dims[0] = inputGrid.getNX();
893 g.dims[1] = inputGrid.getNY();
894 g.dims[2] = inputGrid.getNZ();
896 std::vector<double>
coord = inputGrid.getCOORD( );
897 std::vector<double>
zcorn = inputGrid.getZCORN( );
898 std::vector<int>
actnum = inputGrid.getACTNUM( );
900 g.coord =
coord.data();
901 g.zcorn =
zcorn.data();
904 if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive))
907 const std::vector<double>& minpvv = inputGrid.getMinpvVector();
912 const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
913 std::vector<double> thickness(cartGridSize);
914 for (
size_t i = 0; i < cartGridSize; ++i) {
915 thickness[i] = inputGrid.getCellThickness(i);
917 const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0;
918 mp.process(thickness, z_tolerance, poreVolumes, minpvv,
actnum, opmfil,
zcorn.data());
932 const double z_tolerance = inputGrid.isPinchActive() ?
933 inputGrid.getPinchThresholdThickness() : 0.0;
936 OPM_THROW(std::runtime_error,
"Failed to construct grid.");
942 UnstructuredGridType* createGrid(
const std::vector< int >& n,
const std::vector< double >& dx )
const
944 UnstructuredGridType* cgrid = nullptr ;
945 assert(
int(n.size()) == dim );
957 OPM_THROW(std::runtime_error,
"Failed to construct grid.");
963 #if DUNE_VERSION_LT_REV(DUNE_GRID, 2, 7, 1)
964 using Base::getRealImplementation;
966 typedef typename Traits :: ExtraData ExtraData;
967 ExtraData extraData ()
const {
return this; }
969 template <
class EntitySeed>
970 int corners(
const EntitySeed& seed )
const
972 const int codim = EntitySeed :: codimension;
973 const int index = seed.index();
975 return cellVertices_[ index ].size();
983 template <
class EntitySeed>
985 corner (
const EntitySeed& seed,
const int i )
const
987 const int codim = EntitySeed :: codimension;
990 const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ];
998 const int crners = corners( seed );
999 const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i;
1001 return copyToGlobalCoordinate( grid_.
node_coordinates + GlobalCoordinate :: dimension * faceVertex );
1005 const int coordIndex = GlobalCoordinate :: dimension * seed.index();
1008 return GlobalCoordinate( 0 );
1011 template <
class EntitySeed>
1012 int subEntities(
const EntitySeed& seed,
const int codim )
const
1014 const int index = seed.index();
1015 if( seed.codimension == 0 )
1022 return cellVertices_[ index ].size();
1024 else if( seed.codimension == 1 )
1031 else if ( seed.codimension == dim )
1039 template <
int codim,
class EntitySeedArg >
1040 typename Codim<codim>::EntitySeed
1041 subEntitySeed(
const EntitySeedArg& baseSeed,
const int i )
const
1043 assert( codim >= EntitySeedArg::codimension );
1044 assert( i>= 0 && i<subEntities( baseSeed, codim ) );
1045 typedef typename Codim<codim>::EntitySeed EntitySeed;
1048 if( codim == EntitySeedArg::codimension )
1050 return EntitySeed( baseSeed.index() );
1053 if( EntitySeedArg::codimension == 0 )
1059 else if ( codim == dim )
1061 return EntitySeed( cellVertices_[ baseSeed.index() ][ i ] );
1064 else if ( EntitySeedArg::codimension == 1 && codim == dim )
1069 DUNE_THROW(NotImplemented,
"codimension not available");
1070 return EntitySeed();
1073 template <
int codim>
1074 typename Codim<codim>::EntitySeed
1075 subEntitySeed(
const typename Codim<1>::EntitySeed& faceSeed,
const int i )
const
1077 assert( i>= 0 && i<subEntities( faceSeed, codim ) );
1078 typedef typename Codim<codim>::EntitySeed EntitySeed;
1081 return EntitySeed( faceSeed.index() );
1083 else if ( codim == dim )
1089 DUNE_THROW(NotImplemented,
"codimension not available");
1093 bool hasBoundaryIntersections(
const typename Codim<0>::EntitySeed& seed )
const
1095 const int faces = subEntities( seed, 1 );
1096 for(
int f=0; f<faces; ++f )
1098 const auto faceSeed = this->
template subEntitySeed<1>( seed, f );
1099 if( isBoundaryFace( faceSeed ) )
1105 bool isBoundaryFace(
const int face )
const
1108 const int facePos = 2 * face;
1112 bool isBoundaryFace(
const typename Codim<1>::EntitySeed& faceSeed )
const
1114 assert( faceSeed.isValid() );
1115 return isBoundaryFace( faceSeed.index() );
1118 int boundarySegmentIndex(
const typename Codim<0>::EntitySeed& seed,
const int face )
const
1120 const auto faceSeed = this->
template subEntitySeed<1>( seed, face );
1121 assert( faceSeed.isValid() );
1122 const int facePos = 2 * faceSeed.index();
1129 const std::vector< GeometryType > &geomTypes (
const unsigned int codim )
const
1131 static std::vector< GeometryType > emptyDummy;
1132 if (codim < geomTypes_.size())
1134 return geomTypes_[codim];
1140 template <
class Seed >
1141 GeometryType geometryType(
const Seed& seed )
const
1143 if( Seed::codimension == 0 )
1145 if( cellGeomTypes_.empty() )
1147 assert( geomTypes( Seed::codimension ).
size() == 1 );
1148 return geomTypes( Seed::codimension )[ 0 ];
1152 assert(
static_cast<size_t>(seed.index()) < cellGeomTypes_.size() );
1153 return cellGeomTypes_[ seed.index() ];
1159 if( dim == 3 && Seed::codimension == 1 )
1162 const int nVx = corners( seed );
1164 face = Dune::GeometryTypes::cube(2);
1166 face = Dune::GeometryTypes::simplex(2);
1168 face = Dune::GeometryTypes::none(2);
1173 return geomTypes( Seed::codimension )[ 0 ];
1177 int indexInInside(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1179 return ( grid_.
cell_facetag ) ? cartesianIndexInInside( seed, i ) : i;
1182 int cartesianIndexInInside(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1184 assert( i>= 0 && i<subEntities( seed, 1 ) );
1188 typename Codim<0>::EntitySeed
1189 neighbor(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1191 const int face = this->
template subEntitySeed<1>( seed, i ).index();
1193 if( nb == seed.index() )
1198 typedef typename Codim<0>::EntitySeed EntitySeed;
1199 return EntitySeed( nb );
1203 indexInOutside(
const typename Codim<0>::EntitySeed& seed,
const int i )
const
1208 const int in_inside = cartesianIndexInInside( seed, i );
1209 return in_inside + ((in_inside % 2) ? -1 : 1);
1213 typedef typename Codim<0>::EntitySeed EntitySeed;
1214 EntitySeed nb = neighbor( seed, i );
1215 const int faces = subEntities( seed, 1 );
1216 for(
int face = 0; face<faces; ++ face )
1218 if( neighbor( nb, face ).equals(seed) )
1220 return indexInInside( nb, face );
1223 DUNE_THROW(InvalidStateException,
"inverse intersection not found");
1228 template <
class EntitySeed>
1230 outerNormal(
const EntitySeed& seed,
const int i )
const
1232 const int face = this->
template subEntitySeed<1>( seed, i ).index();
1233 const int normalIdx = face * GlobalCoordinate :: dimension ;
1234 GlobalCoordinate normal = copyToGlobalCoordinate( grid_.
face_normals + normalIdx );
1236 if( nb != seed.index() )
1243 template <
class EntitySeed>
1245 unitOuterNormal(
const EntitySeed& seed,
const int i )
const
1247 const int face = this->
template subEntitySeed<1>( seed, i ).index();
1248 if( seed.index() == grid_.
face_cells[ 2*face ] )
1250 return unitOuterNormals_[ face ];
1254 GlobalCoordinate normal = unitOuterNormals_[ face ];
1260 template <
class EntitySeed>
1261 GlobalCoordinate centroids(
const EntitySeed& seed )
const
1263 if( ! seed.isValid() )
1264 return GlobalCoordinate( 0 );
1266 const int index = GlobalCoordinate :: dimension * seed.index();
1267 const int codim = EntitySeed::codimension;
1268 assert( index >= 0 && index <
size( codim ) * GlobalCoordinate :: dimension );
1274 else if ( codim == 1 )
1278 else if( codim == dim )
1284 DUNE_THROW(InvalidStateException,
"codimension not implemented");
1285 return GlobalCoordinate( 0 );
1289 GlobalCoordinate copyToGlobalCoordinate(
const double* coords )
const
1291 GlobalCoordinate coordinate;
1292 for(
int i=0; i<GlobalCoordinate::dimension; ++i )
1294 coordinate[ i ] = coords[ i ];
1299 template <
class EntitySeed>
1300 double volumes(
const EntitySeed& seed )
const
1302 static const int codim = EntitySeed::codimension;
1303 if( codim == dim || ! seed.isValid() )
1309 assert( seed.isValid() );
1315 else if ( codim == 1 )
1321 DUNE_THROW(InvalidStateException,
"codimension not implemented");
1331 for(
int i=0; i<3; ++i )
1333 cartDims_[ i ] = grid_.
cartdims[ i ];
1336 cellGeomTypes_.clear();
1339 const int numCells =
size( 0 );
1341 cellVertices_.resize( numCells );
1346 typedef std::array<int, 3> KeyType;
1347 std::map< const KeyType, const int > vertexFaceTags;
1348 const int vertexFacePattern [8][3] = {
1359 for(
int i=0; i<8; ++i )
1361 KeyType key; key.fill( 4 );
1362 for(
int j=0; j<dim; ++j )
1364 key[ j ] = vertexFacePattern[ i ][ j ];
1367 vertexFaceTags.insert( std::make_pair( key, i ) );
1370 for (
int c = 0; c < numCells; ++c)
1380 typedef std::map<int,int> vertexmap_t;
1381 typedef typename vertexmap_t :: iterator iterator;
1383 std::vector< vertexmap_t > cell_pts( dim*2 );
1392 const int node = grid_.
face_nodes[ nodepos ];
1393 iterator it = cell_pts[ faceTag ].find( node );
1394 if( it == cell_pts[ faceTag ].end() )
1396 cell_pts[ faceTag ].insert( std::make_pair( node, 1 ) );
1406 typedef std::map< int, std::set<int> > vertexlist_t;
1407 vertexlist_t vertexList;
1409 for(
int faceTag = 0; faceTag<dim*2; ++faceTag )
1411 for( iterator it = cell_pts[ faceTag ].begin(),
1412 end = cell_pts[ faceTag ].end(); it != end; ++it )
1416 if( (*it).second == 1 )
1418 vertexList[ (*it).first ].insert( faceTag );
1423 assert(
int(vertexList.size()) == ( dim == 2 ? 4 : 8) );
1425 cellVertices_[ c ].resize( vertexList.size() );
1426 for(
auto it = vertexList.begin(), end = vertexList.end(); it != end; ++it )
1428 assert( (*it).second.size() == dim );
1429 KeyType key; key.fill( 4 );
1431 std::copy( (*it).second.begin(), (*it).second.end(), key.begin() );
1432 auto vx = vertexFaceTags.find( key );
1433 assert( vx != vertexFaceTags.end() );
1434 if( vx != vertexFaceTags.end() )
1436 if( (*vx).second >=
int(cellVertices_[ c ].size()) )
1437 cellVertices_[ c ].resize( (*vx).second+1 );
1439 cellVertices_[ c ][ (*vx).second ] = (*it).first ;
1445 geomTypes_.resize(dim + 1);
1447 for (
int codim = 0; codim <= dim; ++codim)
1449 tmp = Dune::GeometryTypes::cube(dim - codim);
1450 geomTypes_[codim].push_back(tmp);
1456 int minVx = std::numeric_limits<int>::max();
1458 for (
int c = 0; c < numCells; ++c)
1460 std::set<int> cell_pts;
1466 cell_pts.insert(fnbeg, fnend);
1469 cellVertices_[ c ].resize( cell_pts.size() );
1470 std::copy(cell_pts.begin(), cell_pts.end(), cellVertices_[ c ].begin() );
1471 maxVx = std::max( maxVx,
int( cell_pts.size() ) );
1472 minVx = std::min( minVx,
int( cell_pts.size() ) );
1475 if( minVx == maxVx && maxVx == 4 )
1477 for (
int c = 0; c < numCells; ++c)
1479 assert( cellVertices_[ c ].
size() == 4 );
1480 GlobalCoordinate center( 0 );
1481 GlobalCoordinate p[ dim+1 ];
1482 for(
int i=0; i<dim+1; ++i )
1484 const int vertex = cellVertices_[ c ][ i ];
1486 for(
int d=0; d<dim; ++d )
1493 for(
int d=0; d<dim; ++d )
1498 Dune::GeometryType simplex;
1499 simplex = Dune::GeometryTypes::simplex(dim);
1501 typedef Dune::AffineGeometry< ctype, dim, dimworld> AffineGeometryType;
1502 AffineGeometryType geometry( simplex, p );
1510 for(
int face = 0 ; face < faces; ++face )
1513 const int b = grid_.
face_cells[ 2*face + 1 ];
1515 assert( a >=0 || b >=0 );
1520 GlobalCoordinate centerDiff( 0 );
1523 for(
int d=0; d<dimworld; ++d )
1530 for(
int d=0; d<dimworld; ++d )
1538 for(
int d=0; d<dimworld; ++d )
1545 for(
int d=0; d<dimworld; ++d )
1551 GlobalCoordinate normal( 0 );
1552 for(
int d=0; d<dimworld; ++d )
1557 if( centerDiff.two_norm() < 1e-10 )
1561 if( centerDiff * normal < 0 )
1573 tmp = Dune::GeometryTypes::none(dim);
1574 cellGeomTypes_.resize( numCells );
1575 std::fill( cellGeomTypes_.begin(), cellGeomTypes_.end(), tmp );
1577 bool hasSimplex = false ;
1578 bool hasCube = false ;
1579 bool hasPolyhedron =
false;
1581 for (
int c = 0; c < numCells; ++c)
1583 const int nVx = cellVertices_[ c ].size();
1586 cellGeomTypes_[ c ] = Dune::GeometryTypes::simplex(dim);
1591 cellGeomTypes_[ c ] = Dune::GeometryTypes::cube(dim);
1597 hasPolyhedron =
true;
1601 geomTypes_.resize(dim + 1);
1602 for (
int codim = 0; codim <= dim; ++codim)
1606 tmp = Dune::GeometryTypes::simplex(dim - codim);
1607 geomTypes_[ codim ].push_back( tmp );
1611 tmp = Dune::GeometryTypes::cube(dim - codim);
1612 geomTypes_[ codim ].push_back( tmp );
1614 else if (hasPolyhedron)
1616 tmp = Dune::GeometryTypes::none(dim - codim);
1617 geomTypes_[ codim ].push_back( tmp );
1621 OPM_THROW(std::runtime_error,
"Grid error, unkown geometry type.");
1631 const int normalIdx = face * GlobalCoordinate :: dimension ;
1632 GlobalCoordinate normal = copyToGlobalCoordinate( grid_.
face_normals + normalIdx );
1633 normal /= normal.two_norm();
1634 unitOuterNormals_[ face ] = normal;
1636 if( isBoundaryFace( face ) )
1640 const int facePos = 2 * face ;
1645 grid_.
face_cells[ facePos ] = -nBndSegments_;
1647 else if ( grid_.
face_cells[ facePos+1 ] < 0 )
1649 grid_.
face_cells[ facePos+1 ] = -nBndSegments_;
1655 void print( std::ostream& out,
const UnstructuredGridType& grid )
const
1658 for(
int c=0; c<numCells; ++c )
1660 out <<
"cell " << c <<
" : faces = " << std::endl;
1666 out << f <<
" vx = " ;
1667 while( fnbeg != fnend )
1669 out << *fnbeg <<
" ";
1676 const auto& vx = cellVertices_[ c ];
1677 out <<
"cell " << c <<
" : vertices = ";
1678 for(
size_t i=0; i<vx.size(); ++i )
1679 out << vx[ i ] <<
" ";
1686 UnstructuredGridPtr gridPtr_;
1687 const UnstructuredGridType& grid_;
1690 std::array< int, 3 > cartDims_;
1691 std::vector< std::vector< GeometryType > > geomTypes_;
1692 std::vector< std::vector< int > > cellVertices_;
1694 std::vector< GlobalCoordinate > unitOuterNormals_;
1697 std::vector< GeometryType > cellGeomTypes_;
1703 size_t nBndSegments_;
1715 template<
int dim,
int dimworld,
typename coord_t >
1716 template<
int codim >
1766 template< PartitionIteratorType pitype >
1770 ::template Partition< pitype >::LeafIterator
1773 ::template Partition< pitype >::LevelIterator
1800 #include <opm/grid/polyhedralgrid/persistentcontainer.hh>
1801 #include <opm/grid/polyhedralgrid/cartesianindexmapper.hh>
1802 #include <opm/grid/polyhedralgrid/gridhelpers.hh>
Main OPM-Core grid data structure along with helper functions for construction, destruction and readi...
void destroy_grid(struct UnstructuredGrid *g)
Destroy and deallocate an UnstructuredGrid and all its data.
Definition: UnstructuredGrid.c:32
struct UnstructuredGrid * allocate_grid(size_t ndims, size_t ncells, size_t nfaces, size_t nfacenodes, size_t ncellfaces, size_t nnodes)
Allocate and initialise an UnstructuredGrid where pointers are set to location with correct size.
Definition: UnstructuredGrid.c:87
Routines to construct fully formed grid structures from a simple Cartesian (i.e., tensor product) des...
struct UnstructuredGrid * create_grid_cart2d(int nx, int ny, double dx, double dy)
Form geometrically Cartesian grid in two space dimensions with equally sized cells.
Definition: cart_grid.c:95
struct UnstructuredGrid * create_grid_hexa3d(int nx, int ny, int nz, double dx, double dy, double dz)
Form geometrically Cartesian grid in three space dimensions with equally sized cells.
Definition: cart_grid.c:59
Definition: entityseed.hh:16
Definition: entity.hh:152
Definition: geometry.hh:230
Definition: indexset.hh:24
Definition: intersectioniterator.hh:16
Definition: intersection.hh:20
Definition: iterator.hh:21
Definition: geometry.hh:249
identical grid wrapper
Definition: grid.hh:158
void postAdapt()
Definition: grid.hh:602
Traits::ctype ctype
type of vector coordinates (e.g., double)
Definition: grid.hh:307
int ghostSize(int, int codim) const
obtain size of ghost region for a grid level
Definition: grid.hh:642
bool loadBalance(CommDataHandleIF< DataHandle, Data > &)
rebalance the load each process has to handle
Definition: grid.hh:746
const CollectiveCommunication & comm() const
obtain CollectiveCommunication object
Definition: grid.hh:709
Traits::template Codim< EntitySeed::codimension >::EntityPointer entityPointer(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:806
Traits::template Codim< EntitySeed::codimension >::Entity entity(const EntitySeed &seed) const
obtain EntityPointer from EntitySeed.
Definition: grid.hh:817
bool preAdapt()
Definition: grid.hh:580
int ghostSize(int codim) const
obtain size of ghost region for the leaf grid
Definition: grid.hh:622
bool loadBalance(DofManager &)
rebalance the load each process has to handle
Definition: grid.hh:766
Traits::GlobalIdSet GlobalIdSet
type of global id set
Definition: grid.hh:282
int maxLevel() const
obtain maximal grid level
Definition: grid.hh:412
Traits::LevelIndexSet LevelIndexSet
type of level index set
Definition: grid.hh:270
LevelGridView levelGridView(int) const
View for a grid level for All_Partition.
Definition: grid.hh:790
PolyhedralGrid(const std::vector< int > &n, const std::vector< double > &dx)
constructor
Definition: grid.hh:344
Partition< pitype >::LeafGridView leafGridView() const
View for the leaf grid.
Definition: grid.hh:782
Traits::HierarchicIterator HierarchicIterator
iterator over the grid hierarchy
Definition: grid.hh:224
void scatterData([[maybe_unused]] DataHandle &handle) const
Moves data from the global (all data on process) view to the distributed view.
Definition: grid.hh:881
int size(int, int codim) const
obtain number of entites on a level
Definition: grid.hh:425
Partition< pitype >::LevelGridView levelGridView(int) const
View for a grid level.
Definition: grid.hh:773
Partition< All_Partition >::LevelGridView LevelGridView
View types for All_Partition.
Definition: grid.hh:244
void update()
update grid caches
Definition: grid.hh:836
void switchToDistributedView()
Switch to the distributed view.
Definition: grid.hh:696
int overlapSize(int, int) const
obtain size of overlap region for a grid level
Definition: grid.hh:632
int size(int codim) const
obtain number of leaf entities
Definition: grid.hh:436
Traits::CollectiveCommunication CollectiveCommunication
communicator with all other processes having some part of the grid
Definition: grid.hh:310
Traits::LocalIdSet LocalIdSet
type of local id set
Definition: grid.hh:299
int overlapSize(int) const
obtain size of overlap region for the leaf grid
Definition: grid.hh:613
Traits::LevelIntersectionIterator LevelIntersectionIterator
iterator over intersections with other entities on the same level
Definition: grid.hh:228
GridFamily::Traits Traits
type of the grid traits
Definition: grid.hh:207
PolyhedralGrid(UnstructuredGridPtr &&gridPtr)
constructor
Definition: grid.hh:363
Traits::LeafIntersectionIterator LeafIntersectionIterator
iterator over intersections with other entities on the leaf level
Definition: grid.hh:226
bool adapt(DataHandle &)
Definition: grid.hh:596
int size(int, GeometryType type) const
obtain number of entites on a level
Definition: grid.hh:465
void communicate(DataHandle &, InterfaceType, CommunicationDirection) const
communicate information on leaf entities
Definition: grid.hh:682
size_t numBoundarySegments() const
obtain number of leaf entities
Definition: grid.hh:485
void switchToGlobalView()
Switch to the global view.
Definition: grid.hh:690
bool loadBalance()
rebalance the load each process has to handle
Definition: grid.hh:725
PolyhedralGrid(const UnstructuredGridType &grid)
constructor
Definition: grid.hh:382
LeafGridView leafGridView() const
View for the leaf grid for All_Partition.
Definition: grid.hh:797
int size(GeometryType type) const
returns the number of boundary segments within the macro grid
Definition: grid.hh:474
void communicate(DataHandle &, InterfaceType, CommunicationDirection, int) const
communicate information on a grid level
Definition: grid.hh:661
bool adapt()
Definition: grid.hh:586
Traits::LeafIndexSet LeafIndexSet
type of leaf index set
Definition: grid.hh:260
Transform a corner-point grid ZCORN field to account for MINPV processing.
Definition: MinpvProcessor.hpp:34
Routines to form a complete UnstructuredGrid from a corner-point specification.
struct UnstructuredGrid * create_grid_cornerpoint(const struct grdecl *in, double tol)
Construct grid representation from corner-point specification of a particular geological model.
Definition: cornerpoint_grid.c:164
void compute_geometry(struct UnstructuredGrid *g)
Compute derived geometric primitives in a grid.
Definition: cornerpoint_grid.c:137
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Low-level corner-point processing routines and supporting data structures.
traits structure containing types for a codimension
Definition: grid.hh:1719
Partition< All_Partition >::LevelIterator LevelIterator
type of leaf iterator
Definition: grid.hh:1793
Traits::template Codim< codim >::Entity Entity
type of entity
Definition: grid.hh:1728
Traits::template Codim< codim >::LocalGeometry LocalGeometry
type of local geometry
Definition: grid.hh:1759
Traits::template Codim< codim >::EntityPointer EntityPointer
type of entity pointer
Definition: grid.hh:1734
Partition< All_Partition >::LeafIterator LeafIterator
type of level iterator
Definition: grid.hh:1784
Traits::template Codim< codim >::Geometry Geometry
type of world geometry
Definition: grid.hh:1749
Types for GridView.
Definition: grid.hh:238
Data structure for an unstructured grid, unstructured meaning that any cell may have an arbitrary num...
Definition: UnstructuredGrid.h:99
int * face_nodes
Contains for each face, the indices of its adjacent nodes.
Definition: UnstructuredGrid.h:121
int * face_nodepos
For a face f, face_nodepos[f] contains the starting index for f's nodes in the face_nodes array.
Definition: UnstructuredGrid.h:127
int number_of_cells
The number of cells in the grid.
Definition: UnstructuredGrid.h:109
double * face_areas
Exact or approximate face areas.
Definition: UnstructuredGrid.h:173
int * cell_faces
Contains for each cell, the indices of its adjacent faces.
Definition: UnstructuredGrid.h:146
int number_of_faces
The number of faces in the grid.
Definition: UnstructuredGrid.h:111
double * cell_centroids
Exact or approximate cell centroids, stored consecutively for each cell.
Definition: UnstructuredGrid.h:192
int * cell_facepos
For a cell c, cell_facepos[c] contains the starting index for c's faces in the cell_faces array.
Definition: UnstructuredGrid.h:152
int * cell_facetag
If non-null, this array contains a number for cell-face adjacency indicating the face's position with...
Definition: UnstructuredGrid.h:244
double * cell_volumes
Exact or approximate cell volumes.
Definition: UnstructuredGrid.h:197
int * face_cells
For a face f, face_cells[2*f] and face_cells[2*f + 1] contain the cell indices of the cells adjacent ...
Definition: UnstructuredGrid.h:138
double * face_centroids
Exact or approximate face centroids, stored consecutively for each face.
Definition: UnstructuredGrid.h:168
int number_of_nodes
The number of nodes in the grid.
Definition: UnstructuredGrid.h:113
int cartdims[3]
Contains the size of the logical cartesian structure (if any) of the grid.
Definition: UnstructuredGrid.h:227
int * global_cell
If non-null, this array contains the logical cartesian indices (in a lexicographic ordering) of each ...
Definition: UnstructuredGrid.h:214
double * face_normals
Exact or approximate face normals, stored consecutively for each face.
Definition: UnstructuredGrid.h:184
double * node_coordinates
Node coordinates, stored consecutively for each node.
Definition: UnstructuredGrid.h:160
Raw corner-point specification of a particular geological model.
Definition: preprocess.h:56
const double * coord
Pillar end-points.
Definition: preprocess.h:58
int dims[3]
Cartesian box dimensions.
Definition: preprocess.h:57
const double * zcorn
Corner-point depths.
Definition: preprocess.h:59
const int * actnum
Explicit "active" map.
Definition: preprocess.h:60