diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index a65852d748..80d056d1da 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -91,7 +91,7 @@ namespace Dune typedef Dune::Geometry< dimension-codim, dimensionworld, const Grid, PolyhedralGridGeometry > Geometry; typedef PolyhedralGridLocalGeometry< dimension-codim, dimensionworld, const Grid> LocalGeometryImpl; - typedef Dune::Geometry< dimension-codim, dimension, const Grid, PolyhedralGridLocalGeometry > LocalGeometry; + typedef Dune::Geometry< dimension-codim, dimensionworld, const Grid, PolyhedralGridLocalGeometry > LocalGeometry; typedef PolyhedralGridEntity< codim, dimension, const Grid > EntityImpl; typedef Dune::Entity< codim, dimension, const Grid, PolyhedralGridEntity > Entity; @@ -179,7 +179,8 @@ namespace Dune static UnstructuredGridPtr allocateGrid ( std::size_t nCells, std::size_t nFaces, std::size_t nFaceNodes, std::size_t nCellFaces, std::size_t nNodes ) { - UnstructuredGridType *grid = allocate_grid( dim, nCells, nFaces, nFaceNodes, nCellFaces, nNodes ); + // Note that we here assign a grid of dimension dimworld in order to obtain global coordinates in the correct dimension + UnstructuredGridType *grid = allocate_grid( dimworld, nCells, nFaces, nFaceNodes, nCellFaces, nNodes ); if( !grid ) DUNE_THROW( GridError, "Unable to allocate grid" ); return UnstructuredGridPtr( grid ); @@ -972,26 +973,13 @@ namespace Dune { const int codim = EntitySeed :: codimension; const int index = seed.index(); - switch (codim) - { - case 0: - { - return cellVertices_[ index ].size(); - } - case 1: - { - //return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; - return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ]; - } - case dim: - { - return 1; - } - default: - { - return 0; - } - } + if (codim==0) + return cellVertices_[ index ].size(); + if (codim==1) + return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ]; + if (codim==dim) + return 1; + return 0; } template @@ -999,33 +987,27 @@ namespace Dune corner ( const EntitySeed& seed, const int i ) const { const int codim = EntitySeed :: codimension; - switch (codim) + if (codim==0) { - case 0: - { - const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ]; - return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); - } - case 1: - { - // for faces we need to swap vertices in 3d since in UnstructuredGrid - // those are ordered counter clockwise, for 2d this does not matter - // TODO: Improve this for performance reasons - const int crners = corners( seed ); - const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i; - const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + crner ]; - return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex ); - } - case dim: - { - const int coordIndex = GlobalCoordinate :: dimension * seed.index(); - return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); - } - default: - { - return GlobalCoordinate( 0 ); - } + const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ]; + return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); + } + if (codim==1) + { + // for faces we need to swap vertices in 3d since in UnstructuredGrid + // those are ordered counter clockwise, for 2d this does not matter + // TODO: Improve this for performance reasons + const int crners = corners( seed ); + const int crner = (crners == 4 && EntitySeed :: dimension == 3 && i > 1 ) ? 5 - i : i; + const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + crner ]; + return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex ); } + if (codim==dim) + { + const int coordIndex = GlobalCoordinate :: dimension * seed.index(); + return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); + } + return GlobalCoordinate( 0 ); } template @@ -1034,25 +1016,19 @@ namespace Dune const int index = seed.index(); if( seed.codimension == 0 ) { - switch (codim) - { - case 0: - return 1; - case 1: - return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; - case dim: - return cellVertices_[ index ].size(); - } + if (codim==0) + return 1; + if (codim==1) + return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; + if (codim==dim) + return cellVertices_[ index ].size(); } else if( seed.codimension == 1 ) { - switch (codim) - { - case 1: - return 1; - case dim: - return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ]; - } + if (codim==1) + return 1; + if (codim==dim) + return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ]; } else if ( seed.codimension == dim ) { @@ -1532,7 +1508,6 @@ namespace Dune // check face normals { - typedef Dune::FieldVector< double, dim > Coordinate; const int faces = grid_.number_of_faces; for( int face = 0 ; face < faces; ++face ) { @@ -1544,41 +1519,41 @@ namespace Dune if( grid_.face_areas[ face ] < 0 ) std::abort(); - Coordinate centerDiff( 0 ); + GlobalCoordinate centerDiff( 0 ); if( b >= 0 ) { - for( int d=0; d= 0 ) { - for( int d=0; d 2 ) + for (int codim = 0; codim <= dim; ++codim) { - for( const auto& elemType : geomTypes_[ 0 ] ) + if( hasSimplex ) { - if( elemType.isSimplex() ) - tmp = Dune::GeometryTypes::simplex(2); - else if ( elemType.isCube() ) - tmp = Dune::GeometryTypes::cube(2); - else - tmp = Dune::GeometryTypes::none(2); - geomTypes_[ 1 ].push_back( tmp ); + tmp = Dune::GeometryTypes::simplex(dim - codim); + geomTypes_[ codim ].push_back( tmp ); + } + else if ( hasCube ) + { + tmp = Dune::GeometryTypes::cube(dim - codim); + geomTypes_[ codim ].push_back( tmp ); + } + else if (hasPolyhedron) + { + tmp = Dune::GeometryTypes::none(dim - codim); + geomTypes_[ codim ].push_back( tmp ); + } + else + { + OPM_THROW(std::runtime_error, "Grid error, unkown geometry type."); } } @@ -1691,17 +1652,6 @@ namespace Dune } } } - - /* - for( int i=0; i<= dim ; ++ i ) - { - for( const auto& geomType : geomTypes_[ i ] ) - { - std::cout << "Codim " << i << " type = " << geomType << std::endl; - } - } - */ - //print( std::cout, grid_ ); } void print( std::ostream& out, const UnstructuredGridType& grid ) const diff --git a/opm/grid/polyhedralgrid/intersection.hh b/opm/grid/polyhedralgrid/intersection.hh index b519873422..0d4aa325be 100644 --- a/opm/grid/polyhedralgrid/intersection.hh +++ b/opm/grid/polyhedralgrid/intersection.hh @@ -50,6 +50,12 @@ namespace Dune intersectionIdx_( -1 ) {} + PolyhedralGridIntersection ( ) + : data_( ), + seed_(), + intersectionIdx_( -1 ) + {} + PolyhedralGridIntersection ( ExtraData data, const EntitySeed& seed, const int intersectionIdx ) : data_( data ), seed_( seed ), diff --git a/tests/test_polyhedralgrid.cpp b/tests/test_polyhedralgrid.cpp index f6178007b9..2f8dd9d408 100644 --- a/tests/test_polyhedralgrid.cpp +++ b/tests/test_polyhedralgrid.cpp @@ -340,5 +340,118 @@ int main(int argc, char** argv ) Dune::GridPtr< Grid > gridPtr( dgfFile ); gridcheck( *gridPtr ); } + + { + // Test 2D grid embedded in a 3D domain + const std::string gridFileName = "polyhedral_grid_test.txt"; + std::stringstream gridFile; + std::ofstream out(gridFileName); + + gridFile << "3 1 3 3 6 3 0 0" << std::endl; + gridFile << "0 0 0" << std::endl; + gridFile << "0.0 0.0 1.0 0.0 1.0 0.0 0.0 1.0 1.0" << std::endl; + gridFile << "0 2 4 6" << std::endl; + gridFile << "0 1 0 2 1 2" << std::endl; + gridFile << "0 -1 0 -1 0 -1" << std::endl; + gridFile << "1.4142135623730951 1.0 1.0" << std::endl; + gridFile << "0.0 0.5 0.5 0.0 0.5 1.0 0.0 1.0 0.5" << std::endl; + gridFile << "-1.1102230246251565e-16 -1.0 -1.0 1.1102230246251565e-16 0.0 1.0 0.0 1.0 0.0" << std::endl; + gridFile << "0 3" << std::endl; + gridFile << "0 1 2" << std::endl; + gridFile << "0.5" << std::endl; + gridFile << "0.0 0.6666666666666666 0.6666666666666666" << std::endl; + + out << gridFile.str(); + out.close(); + + const char* c_str = gridFileName.c_str(); + typedef Dune::PolyhedralGrid< 2, 3, float > Grid; + UnstructuredGrid* grid = read_grid(c_str); + std::remove(c_str); + + if (!grid) { + std::string msg = "RuntimeError: UnstructuredGrid could not read grid file"; + throw std::runtime_error(msg); + } + // check different coordinate field type here + Grid::UnstructuredGridPtr ugPtr( grid ); + Dune::GridPtr< Grid > gridPtr( new Grid(*ugPtr) ); + gridcheck( *gridPtr ); + } + + { + // Test 1D grid embedded in a 3D domain + const std::string gridFileName = "polyhedral_grid_test.txt"; + std::stringstream gridFile; + std::ofstream out(gridFileName); + gridFile << "3 2 3 3 3 4 0 0" << std::endl; + gridFile << "0 0 0" << std::endl; + gridFile << "0.0 0.0 0.0 0.5 0.5 1.0 1.0 1.0 2.0" << std::endl; + gridFile << "0 1 2 3" << std::endl; + gridFile << "0 1 2" << std::endl; + gridFile << "-1 0 0 1 1 -1" << std::endl; + gridFile << "1.0 1.0 1.0" << std::endl; + gridFile << "0.0 0.0 0.0 0.5 0.5 1.0 1.0 1.0 2.0" << std::endl; + gridFile << "0.4082482904638631 0.4082482904638631 0.8164965809277261 0.4082482904638631 0.4082482904638631 0.8164965809277261 0.4082482904638631 0.4082482904638631 0.8164965809277261" << std::endl; + gridFile << "0 2 4" << std::endl; + gridFile << "0 1 1 2" << std::endl; + gridFile << "1.224744871391589 1.224744871391589" << std::endl; + gridFile << "0.25 0.25 0.5 0.75 0.75 1.5" << std::endl; + + out << gridFile.str(); + out.close(); + + const char* c_str = gridFileName.c_str(); + UnstructuredGrid* grid = read_grid(c_str); + std::remove(c_str); + + if (!grid) { + std::string msg = "RuntimeError: UnstructuredGrid could not read grid file"; + throw std::runtime_error(msg); + } + // check different coordinate field type here + typedef Dune::PolyhedralGrid< 1, 3, float > Grid; + Grid::UnstructuredGridPtr ugPtr( grid ); + Dune::GridPtr< Grid > gridPtr( new Grid(*ugPtr) ); + gridcheck( *gridPtr ); + } + + { + // Test 1D grid embedded in a 2D domain + const std::string gridFileName = "polyhedral_grid_test.txt"; + std::stringstream gridFile; + std::ofstream out(gridFileName); + gridFile << "2 2 3 3 3 4 0 0" << std::endl; + gridFile << "0 0" << std::endl; + gridFile << "0.0 -0.0 0.5 -0.5 1.0 -1.0" << std::endl; + gridFile << "0 1 2 3" << std::endl; + gridFile << "0 1 2" << std::endl; + gridFile << "-1 0 0 1 1 -1" << std::endl; + gridFile << "1.0 1.0 1.0" << std::endl; + gridFile << "0.0 -0.0 0.5 -0.5 1.0 -1.0" << std::endl; + gridFile << "0.7071067811865475 -0.7071067811865475 0.7071067811865475 -0.7071067811865475 0.7071067811865475 -0.7071067811865475" << std::endl; + gridFile << "0 2 4" << std::endl; + gridFile << "0 1 1 2" << std::endl; + gridFile << "0.7071067811865476 0.7071067811865476" << std::endl; + gridFile << "0.25 -0.25 0.75 -0.75" << std::endl; + + out << gridFile.str(); + out.close(); + + const char* c_str = gridFileName.c_str(); + UnstructuredGrid* grid = read_grid(c_str); + std::remove(c_str); + + if (!grid) { + std::string msg = "RuntimeError: UnstructuredGrid could not read grid file"; + throw std::runtime_error(msg); + } + // check different coordinate field type here + typedef Dune::PolyhedralGrid< 1, 2, float > Grid; + Grid::UnstructuredGridPtr ugPtr( grid ); + Dune::GridPtr< Grid > gridPtr( new Grid(*ugPtr) ); + gridcheck( *gridPtr ); + } + return 0; }