diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 295be4f2c3..19be201c3e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -58,6 +58,10 @@ list (APPEND MAIN_SOURCE_FILES opm/grid/utility/cartesianToCompressed.cpp opm/grid/utility/StopWatch.cpp opm/grid/utility/WachspressCoord.cpp + opm/grid/verteq/topsurf.cpp + opm/grid/verteq/nav.cpp + opm/grid/verteq/utility/exc.cpp + opm/grid/verteq/utility/runlen.cpp ) if (opm-common_FOUND) @@ -173,6 +177,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/grid/common/WellConnections.hpp opm/grid/common/ZoltanGraphFunctions.hpp opm/grid/common/ZoltanPartition.hpp + opm/grid/polyhedralgrid.hh opm/grid/polyhedralgrid/capabilities.hh opm/grid/polyhedralgrid/cartesianindexmapper.hh opm/grid/polyhedralgrid/declaration.hh @@ -184,7 +189,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/grid/polyhedralgrid/gridhelpers.hh opm/grid/polyhedralgrid/grid.hh opm/grid/polyhedralgrid/gridview.hh - opm/grid/polyhedralgrid.hh opm/grid/polyhedralgrid/idset.hh opm/grid/polyhedralgrid/indexset.hh opm/grid/polyhedralgrid/intersection.hh @@ -221,4 +225,9 @@ list (APPEND PUBLIC_HEADER_FILES opm/grid/utility/OpmParserIncludes.hpp opm/grid/utility/platform_dependent/disable_warnings.h opm/grid/utility/platform_dependent/reenable_warnings.h + opm/grid/verteq/topsurf.hpp + opm/grid/verteq/nav.hpp + opm/grid/verteq/utility/visibility.h + opm/grid/verteq/utility/exc.hpp + opm/grid/verteq/utility/runlen.hpp ) diff --git a/opm/grid/UnstructuredGrid.c b/opm/grid/UnstructuredGrid.c index 271376795c..6ede525c86 100644 --- a/opm/grid/UnstructuredGrid.c +++ b/opm/grid/UnstructuredGrid.c @@ -166,6 +166,67 @@ allocate_grid(size_t ndims , return G; } +#define SWAP(a,b) { int tmp = a; a = b; b =tmp; } + +void print_grid(const struct UnstructuredGrid *grd ) +{ + int dim = grd->dimensions; + for(int i=0; inumber_of_nodes; ++i ) + { + printf("Vx (%d) = ( ",i); + for( int d=0; dnode_coordinates[i*dim + d ]); + } + printf(")\n"); + } + + for( int f=0; fnumber_of_faces; ++f ) + { + //for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) + printf(" Face (%d) (",f); + for( int vx=grd->face_nodepos[ f ]; vxface_nodepos[ f+1 ]; ++vx ) + { + printf(" %d ",grd->face_nodes[ vx ]); + } + printf(")"); + int bnd = grd->face_cells[ 2*f ] < 0 || grd->face_cells[ 2*f+1 ] < 0; + if( bnd ) + printf(" boundary"); + printf("\n"); + } + + /* + if( dim == 2 ) + { + for( int c=0; cnumber_of_cells; ++c ) + { + int f=grd->cell_facepos[ c ]; + SWAP( grd->cell_faces[ f+1 ], grd->cell_faces[ f+2 ] ); + SWAP( grd->cell_facetag[ f+1 ], grd->cell_facetag[ f+2 ] ); + } + } + */ + + for( int c=0; cnumber_of_cells; ++c ) + { + printf("Cell %d = (", c); + //for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) + for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) + { + printf(" %d ",grd->cell_faces[ f ]); + + /* + for( int vx=grd->face_nodepos[ f ]; vxface_nodepos[ f+1 ]; ++vx ) + { + printf(" %d ",grd->face_nodes[ vx ]); + } + printf(")\n"); + */ + } + printf(")\n"); + } +} #define GRID_NMETA 6 #define GRID_NDIMS 0 diff --git a/opm/grid/UnstructuredGrid.h b/opm/grid/UnstructuredGrid.h index c2f8a543ae..48cdf5ede9 100644 --- a/opm/grid/UnstructuredGrid.h +++ b/opm/grid/UnstructuredGrid.h @@ -299,6 +299,8 @@ allocate_grid(size_t ndims , size_t ncellfaces, size_t nnodes ); +void print_grid(const struct UnstructuredGrid *); + /** Will allocate storage internally in the grid object to hold a copy diff --git a/opm/grid/common/VerteqColumnUtility.hpp b/opm/grid/common/VerteqColumnUtility.hpp new file mode 100644 index 0000000000..c0472ff799 --- /dev/null +++ b/opm/grid/common/VerteqColumnUtility.hpp @@ -0,0 +1,56 @@ +#ifndef OPM_VERTEQCOLUMNUTILITY_HEADER +#define OPM_VERTEQCOLUMNUTILITY_HEADER + +#include +#include + +#include +#include + + +namespace Dune +{ + + class ColumnCell + { + typedef Opm::TopSurf TopSurfaceGridType; + const TopSurfaceGridType& topSurf_; + const int index_; + + public: + ColumnCell( const TopSurfaceGridType& topSurf, const int index ) + : topSurf_( topSurf ), + index_( index ) + {} + + int index () const { return index_; } + + double dz () const { return topSurf_.dz[ index_ ]; } + double h () const { return topSurf_.h [ index_ ]; } + + int maxVertRes() const { return topSurf_.max_vert_res; } + + // this can be used to access an entity in the original 3D grid + int fineCellIndex() const { return topSurf_.fine_col[ index_ ]; } + }; + + + /** \brief Interface class to access the logical Cartesian grid as used in industry + standard simulator decks. + */ + template< class Grid > + class VerteqColumnUtility + { + public: + /** \brief dimension of the grid */ + static const int dimension = Grid :: dimension ; + + /** \brief constructor taking grid */ + explicit VerteqColumnUtility( const Grid& ) + { + DUNE_THROW(InvalidStateException,"CartesianIndexMapper not specialized for given grid"); + } + }; + +} // end namespace Opm +#endif diff --git a/opm/grid/polyhedralgrid/geometry.hh b/opm/grid/polyhedralgrid/geometry.hh index 2ac43c4194..ea81902a01 100644 --- a/opm/grid/polyhedralgrid/geometry.hh +++ b/opm/grid/polyhedralgrid/geometry.hh @@ -191,6 +191,7 @@ namespace Dune return storage_.volume(); } + JacobianTransposed jacobianTransposed ( const LocalCoordinate & local ) const { if( geometryImpl_ ) diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index a65852d748..53e6623800 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -30,6 +30,7 @@ #include #include #include +#include // Re-enable warnings. #include @@ -43,6 +44,10 @@ #include #include +#include + +#include + namespace Dune { @@ -163,6 +168,7 @@ namespace Dune public: typedef UnstructuredGrid UnstructuredGridType; + typedef Opm::TopSurf TopSurfaceGridType; protected: struct UnstructuredGridDeleter @@ -174,6 +180,7 @@ namespace Dune }; public: + typedef PolyhedralMesh< dim, dimworld, coord_t, int > PolyhedralMeshType; typedef std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > UnstructuredGridPtr; static UnstructuredGridPtr @@ -344,6 +351,7 @@ namespace Dune const std::vector< double >& dx ) : gridPtr_( createGrid( n, dx ) ), grid_( *gridPtr_ ), + topSurfaceGrid_( nullptr ), comm_( MPIHelper::getCommunicator()), leafIndexSet_( *this ), globalIdSet_( *this ), @@ -362,6 +370,8 @@ namespace Dune explicit PolyhedralGrid ( UnstructuredGridPtr &&gridPtr ) : gridPtr_( std::move( gridPtr ) ), grid_( *gridPtr_ ), + topSurfaceGrid_( nullptr ), + //polyhedralMesh_( grid_ ), comm_( MPIHelper::getCommunicator() ), leafIndexSet_( *this ), globalIdSet_( *this ), @@ -381,6 +391,7 @@ namespace Dune explicit PolyhedralGrid ( const UnstructuredGridType& grid ) : gridPtr_(), grid_( grid ), + topSurfaceGrid_( nullptr ), comm_( MPIHelper::getCommunicator() ), leafIndexSet_( *this ), globalIdSet_( *this ), @@ -390,6 +401,27 @@ namespace Dune init(); } + /** \brief constructor + * + * The references to ug are stored in the grid. + * Therefore, they must remain valid until the grid is destroyed. + * + * \param[in] ug UnstructuredGrid reference + */ + explicit PolyhedralGrid ( const TopSurfaceGridType& topSurf ) + : gridPtr_(), + grid_( static_cast< UnstructuredGridType > ( topSurf ) ), + topSurfaceGrid_( &topSurf ), + comm_( *this ), + leafIndexSet_( *this ), + globalIdSet_( *this ), + localIdSet_( *this ), + nBndSegments_( 0 ) + { + // std::cout << "Creating TopSurfaceGrid" << topSurfaceGrid_ << std::endl; + init(); + } + /** \} */ /** \name Casting operators @@ -398,6 +430,8 @@ namespace Dune /** \} */ + const TopSurfaceGridType* topSurfaceGrid() const { return topSurfaceGrid_; } + /** \name Size Methods * \{ */ @@ -434,6 +468,7 @@ namespace Dune */ int size ( int codim ) const { + //return polyhedralMesh_.size( codim ); if( codim == 0 ) { return grid_.number_of_cells; @@ -1353,6 +1388,8 @@ namespace Dune protected: void init () { + //std::cout << "PolyhedralGrid init" << std::endl; + // copy Cartesian dimensions for( int i=0; i<3; ++i ) { @@ -1465,6 +1502,11 @@ namespace Dune cellVertices_[ c ][ (*vx).second ] = (*it).first ; } } + + for( int c=0; c cartDims_; std::vector< std::vector< GeometryType > > geomTypes_; @@ -1852,5 +1895,6 @@ namespace Dune #include #include #include +#include #endif // #ifndef DUNE_POLYHEDRALGRID_GRID_HH diff --git a/opm/grid/polyhedralgrid/gridfactory.hh b/opm/grid/polyhedralgrid/gridfactory.hh index 80fef4429a..936e20ffcb 100644 --- a/opm/grid/polyhedralgrid/gridfactory.hh +++ b/opm/grid/polyhedralgrid/gridfactory.hh @@ -18,7 +18,6 @@ namespace Dune // GridFactory for PolyhedralGrid // --------------------------------- - template< int dim, int dimworld, class coord_t > class GridFactory< PolyhedralGrid< dim, dimworld, coord_t > > : public GridFactoryInterface< PolyhedralGrid< dim, dimworld, coord_t > > diff --git a/opm/grid/polyhedralgrid/idset.hh b/opm/grid/polyhedralgrid/idset.hh index f001cee153..5bc91e7a87 100644 --- a/opm/grid/polyhedralgrid/idset.hh +++ b/opm/grid/polyhedralgrid/idset.hh @@ -48,6 +48,15 @@ namespace Dune } } +#if ! DUNE_VERSION_NEWER(DUNE_GRID,2,4) + //! id meethod for entity and specific codim + template< int codim > + IdType id ( const typename Traits::template Codim< codim >::EntityPointer &entityPointer ) const + { + return id( *entityPointer ); + } +#endif + //! id method of all entities template< class Entity > IdType id ( const Entity &entity ) const diff --git a/opm/grid/polyhedralgrid/intersection.hh b/opm/grid/polyhedralgrid/intersection.hh index b519873422..1ab538ce75 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/opm/grid/polyhedralgrid/polyhedralmesh.hh b/opm/grid/polyhedralgrid/polyhedralmesh.hh new file mode 100644 index 0000000000..7bcc7f12d8 --- /dev/null +++ b/opm/grid/polyhedralgrid/polyhedralmesh.hh @@ -0,0 +1,253 @@ +#ifndef DUNE_POLYHEDRALGRID_MESH_HH +#define DUNE_POLYHEDRALGRID_MESH_HH + +#include + +namespace Dune +{ + template + class PolyhedralMesh + { + static constexpr int volumeData = 0; + static constexpr int centroidData = 1; + static constexpr int normalData = centroidData + dimworld; + + static constexpr int elementData = centroidData + dimworld; + static constexpr int faceData = normalData + dimworld; + public: + typedef T ctype; + typedef I Index; + + typedef Dune::FieldVector< ctype, dimworld > GlobalCoordinate ; + + PolyhedralMesh () + : centroids_( dim+1 ), + faceNormals_(), + subEntity_( 1 ), + subEntityPos_( 1, std::vector(1, Index(0)) ), + topology_( dim+1 ), + topologyPos_( dim+1, std::vector(1, Index(0)) ), + volumes_( dim-1 ), + geomTypes_() + {} + + template + PolyhedralMesh (const UnstructuredGrid& ug) + : centroids_( dim+1 ), + faceNormals_(), + subEntity_( 1 ), + subEntityPos_( 1, std::vector(1, Index(0)) ), + topology_( dim+1 ), + topologyPos_( dim+1, std::vector(1, Index(0)) ), + volumes_( dim-1 ), + geomTypes_() + { + centroids_[ dim ].resize( ug.number_of_nodes ); + centroids_[ 0 ].resize( ug.number_of_cells ); + volumes_[ 0 ].resize( ug.number_of_cells ); + + centroids_[ 1 ].resize( ug.number_of_faces ); + faceNormals_.resize( ug.number_of_faces ); + volumes_[ 1 ].resize( ug.number_of_faces ); + GeometryType tmp; + tmp.makeNone( dim ); + geomTypes_.resize( ug.number_of_cells, tmp ); + + int id = 0; + for( int i=0; i::max(); + + const int numCells = ug.number_of_cells; + topology_[ 0 ].reserve( numCells * 8 ); + topologyPos_[ 0 ].resize( numCells + 1 ); + topologyPos_[ 0 ][ 0 ] = 0; + for (int c = 0; c < numCells; ++c) + { + topologyPos_[ 0 ][ c ] = topology_[ 0 ].size(); + std::set cell_pts; + for (int hf=ug.cell_facepos[ c ]; hf < ug.cell_facepos[c+1]; ++hf) + { + int f = ug.cell_faces[ hf ]; + const int* fnbeg = ug.face_nodes + ug.face_nodepos[f]; + const int* fnend = ug.face_nodes + ug.face_nodepos[f+1]; + cell_pts.insert(fnbeg, fnend); + } + + for( const auto& vertex : cell_pts ) + { + topology_[ 0 ].push_back( vertex ); + } + maxVx = std::max( maxVx, int( cell_pts.size() ) ); + minVx = std::min( minVx, int( cell_pts.size() ) ); + } + topologyPos_[ 0 ][ numCells ] = topology_[ 0 ].size(); + } + + /** \brief return number of cells in the mesh */ + Index size( const int codim ) const { return centroids_[ codim ].size(); } + + GlobalCoordinate& coordinate( const Index idx ) { return coordinates()[ idx ]; } + const GlobalCoordinate& coordinate( const Index idx ) const { return coordinates()[ idx ]; } + + std::pair< Index, Index* > entity( const Index en, const int codim ) + { + return std::make_pair( topologyPos_[ codim ][ en+1 ] - topologyPos_[ codim ][ en ], topology_[ codim ].data() + topologyPos_[ codim ][ en ] ); + } + + /** \brief return number of sub entities, i.e. corners of a face or element + * \param entity entity index, such as face number or element number + * \param codim codimension of the entity + * */ + Index subEntities( const Index entity, const int codim ) const + { + return topologyPos_[ codim ][ entity+1 ] - topologyPos_[ codim ][ entity ]; + } + + /** \brief return sub entities, i.e. corners of a face or element + * \param entity entity index, such as face number or element number + * \param i i-th vertex requested + * \param codim codimension of the entity + * */ + Index subEntity( const Index entity, const int i, const int codim ) const + { + return topology_[ codim ][ topologyPos_[ codim ][ entity ] + i ]; + } + + std::pair< Index, Index* > element( const Index en ) + { + return entity( en, 0 ); + } + + ctype volume( const Index& entity, const int codim ) const + { + return volumes_[ codim ][ entity ]; + } + + GlobalCoordinate& center( const Index& entity, const int codim ) + { + return centroids_[ codim ][ entity ]; + } + + const GlobalCoordinate& center( const Index& entity, const int codim ) const + { + return centroids_[ codim ][ entity ]; + } + + GlobalCoordinate& faceNormal( const Index& face) + { + return faceNormals_[ face ]; + } + + const GlobalCoordinate& faceNormal( const Index& face) const + { + return faceNormals_[ face ]; + } + + void insertVertex( const GlobalCoordinate& vertex ) + { + coordinates().push_back( vertex ); + } + + void insertElement( const GeometryType& type, const std::vector& items ) + { + // could be face or element + const int codim = dim - type.dim(); + if( codim == 1 ) // face + { + insertItems( items, topology_[ codim ], topologyPos_[ codim ] ); + } + else if( codim == 0 ) + { + // store geometry type of element + geomTypes_.push_back( type ); + // for polyhedral cells items are the faces + if( type.isNone() ) + { + insertItems( items, subEntity_[ codim ], subEntityPos_[ codim ] ); + } + else + { + // otherwise items are vertices + insertItems( items, topology_[ codim ], topologyPos_[ codim ] ); + } + } + + // TODO: compute volume and normals etc. + } + + protected: + std::vector< GlobalCoordinate >& coordinates() { return centroids_[ dim ]; } + const std::vector< GlobalCoordinate >& coordinates() const { return centroids_[ dim ]; } + + void insertItems( const std::vector& items, + std::vector< Index >& entities, + std::vector< Index >& entityPos ) + { + const int iSize = items.size(); + entities.reserve( entities.size() + items.size() ); + assert( entityPos.back() == items.size() ); + for( int i=0; i > centroids_; + std::vector< GlobalCoordinate > faceNormals_; + + std::vector< std::vector< Index > > subEntity_; + std::vector< std::vector< Index > > subEntityPos_; + + std::vector< std::vector< Index > > topology_; + std::vector< std::vector< Index > > topologyPos_; + + std::vector< Index > faceNeighbors_; + + std::vector< std::vector< ctype > > volumes_; // mostly codim 0 and 1 + std::vector< GeometryType > geomTypes_; + }; + +}// end namespace Dune + +#endif diff --git a/opm/grid/polyhedralgrid/verteqcolumnutility.hh b/opm/grid/polyhedralgrid/verteqcolumnutility.hh new file mode 100644 index 0000000000..f3ea18dda0 --- /dev/null +++ b/opm/grid/polyhedralgrid/verteqcolumnutility.hh @@ -0,0 +1,105 @@ +#ifndef OPM_POLYHEDRALGRID_VERTEQCOLUMNUTILITY_HEADER +#define OPM_POLYHEDRALGRID_VERTEQCOLUMNUTILITY_HEADER + +#include +#include + +#include +#include + +#include + +namespace Dune +{ + /** \brief Interface class to access the columns in a verteq grid + */ + template< int dim, int dimworld, typename coord_t > + class VerteqColumnUtility< PolyhedralGrid< dim, dimworld, coord_t > > + { + typedef PolyhedralGrid< dim, dimworld, coord_t > Grid; + typedef typename Grid :: TopSurfaceGridType TopSurfaceGridType; + + const TopSurfaceGridType* topSurfaceGrid_; + + // iterator through columns, + // ForwardIteratorFacade provides the standard iterator methods + class Iterator : public ForwardIteratorFacade< Iterator, ColumnCell, ColumnCell > + { + // see VerteqColumnUtility.hpp + typedef ColumnCell ColumnCellType; + + const TopSurfaceGridType* topSurfaceGrid_; + int idx_; + + public: + Iterator( const TopSurfaceGridType* topSurfaceGrid, const int idx ) + : topSurfaceGrid_( topSurfaceGrid ), + idx_( idx ) + {} + + ColumnCellType dereference() const + { + // -1 means non-valid iterator + assert( idx_ >= 0 ); + // if iterator is valid topSurf must be available. + assert( topSurfaceGrid_ ); + return ColumnCellType( *topSurfaceGrid_, idx_ ); + } + + void increment () + { + // -1 means non-valid iterator + assert( idx_ >= 0 ); + ++idx_; + } + + bool equals ( const Iterator& other ) const + { + return idx_ == other.idx_; + } + }; + + public: + typedef typename Grid :: template Codim< 0 > :: Entity Entity; + + typedef Iterator IteratorType; + + /** \brief dimension of the grid */ + static const int dimension = Grid :: dimension ; + + /** \brief constructor taking grid */ + explicit VerteqColumnUtility( const Grid& grid ) + : topSurfaceGrid_( grid.topSurfaceGrid() ) + { + } + + IteratorType begin( const Entity& entity ) const + { + if( topSurfaceGrid_ ) + { + const int colIdx = topSurfaceGrid_->col_cellpos[ entity.impl().index() ]; + std::cout << "begin " << colIdx << std::endl; + return IteratorType( topSurfaceGrid_, colIdx ); + } + else + return end( entity ); + } + + IteratorType end( const Entity& entity ) const + { + if( topSurfaceGrid_ ) + { + const int colIdx = topSurfaceGrid_->col_cellpos[ entity.impl().index()+1 ]; + std::cout << "end " << colIdx << std::endl; + return IteratorType( topSurfaceGrid_, colIdx ); + } + else + { + std::cout << "end no top surf " << std::endl; + return IteratorType( nullptr, -1 ); + } + } + }; + +} // end namespace Dune +#endif diff --git a/opm/grid/verteq/nav.cpp b/opm/grid/verteq/nav.cpp new file mode 100644 index 0000000000..8efc315cc5 --- /dev/null +++ b/opm/grid/verteq/nav.cpp @@ -0,0 +1,92 @@ +// Copyright (C) 2013 Uni Research AS +// This file is licensed under the GNU General Public License v3.0 +#include +#include +#include +using namespace std; + +const Dir Dir::DEC (0); +const Dir Dir::INC (1); + +const Dim2D Dim2D::X (0); +const Dim2D Dim2D::Y (1); +const Dim3D Dim3D::Z (2); + +template Side +Side ::from_tag (int tag) { + // direction is the minor bits in the enumeration + const div_t bits = div (tag, Dir::COUNT); + const int dir_val = bits.rem; + const int dim_val = bits.quot; + return Side (Dim (dim_val), Dir (dir_val)); +} + +// template instantiation to satisfy linker +template Side Side ::from_tag (int); +template Side Side ::from_tag (int); + +// these needs to be initialized here instead of in the header +// because vector::resize takes a reference to the data and not +// a value as a parameter (to avoid copying) +const int Cart2D::NO_ELEM = -1; +const int Cart2D::NO_FACE = -1; +const int Cart2D::NO_NODE = -1; + +// enumeration of all possible sides +template <> +const Side Side ::ALL[] = { + Side (Dim2D::X, Dir::DEC), // I- + Side (Dim2D::X, Dir::INC), // I+ + Side (Dim2D::Y, Dir::DEC), // J- + Side (Dim2D::Y, Dir::INC) // J+ +}; + +template <> +const Side Side ::ALL[] = { + Side (Dim2D::X, Dir::DEC), // I- + Side (Dim2D::X, Dir::INC), // I+ + Side (Dim2D::Y, Dir::DEC), // J- + Side (Dim2D::Y, Dir::INC), // J+ + Side (Dim3D::Z, Dir::DEC), // K- + Side (Dim3D::Z, Dir::INC) // K+ +}; + +// sides that exists as standalone constants +const Side3D UP (Dim3D::Z, Dir::DEC); +const Side3D DOWN (Dim3D::Z, Dir::INC); + +// print carriers (for debugging) +ostream& operator << (ostream& os, const Coord2D& c) { + return (os << '(' << c.i () << ',' << c.j () << ')'); // e.g. "(3,2)" +} + +ostream& operator << (ostream& os, const Coord3D& c) { + // e.g. "(3,2,5)" + return (os << '(' << c.i () << ',' << c.j () << ',' << c.k () << ')'); +} + +static const char DIR_NAMES[] = {'-', '+'}; +static const char DIM_NAMES[] = {'I', 'J', 'K'}; + +ostream& operator << (ostream& os, const Dir& d) { + return (os << DIR_NAMES[d.val]); // e.g. '-' +} + +ostream& operator << (ostream& os, const Dim2D& d) { + return (os << DIM_NAMES[d.val]); // e.g. 'I' +} + +ostream& operator << (ostream& os, const Side2D& s) { + return (os << s.dim () << s.dir ()); // e.g. "I-" +} + +ostream& operator << (ostream& os, const Side3D& s) { + return (os << s.dim () << s.dir ()); // e.g. "I-" +} + +ostream& operator << (ostream& os, const Corn3D& c) { + // e.g. "(I-,J+)" + return (os << '(' << DIM_NAMES[Dim3D::X.val] << c.i () << ',' + << DIM_NAMES[Dim3D::Y.val] << c.j () << ',' + << DIM_NAMES[Dim3D::Z.val] << c.k () << ')'); +} diff --git a/opm/grid/verteq/nav.hpp b/opm/grid/verteq/nav.hpp new file mode 100644 index 0000000000..61583127dc --- /dev/null +++ b/opm/grid/verteq/nav.hpp @@ -0,0 +1,425 @@ +#ifndef OPM_VERTEQ_NAV_HPP_INCLUDED +#define OPM_VERTEQ_NAV_HPP_INCLUDED + +// Copyright (C) 2013 Uni Research AS +// This file is licensed under the GNU General Public License v3.0 + +#ifndef OPM_GRID_HEADER_INCLUDED +#include +#endif + +#include // div_t +#include // ostream + +/** + * There are three types of indices used in this module: + * + * (1) Cartesian coordinate + * (2) Cartesian index + * (3) Global identity + * + * The Cartesian coordinate is an (i,j)-tuple into the cornerpoint grid + * structure. + * + * The Cartesian index, is a flat integer which has is + * determined solely by the structure of the grid, regardless of whether + * there are any active elements. It can be calculated from the coordinates + * if the extent of the grid is known. We use this to efficiently store + * data for each possible position in the grid without using a predefined + * size or dynamic structure for each column. + * + * The global identity is the index which is assigned to it in the grid + * structure, after inactive elements are discarded. This is the 'identity' + * of the cell in the rest of the simulator. Cells that aren't active are + * not assigned an identity. + * + * The value types defined here provide a way to address location in + * the grid in a type-safe manner to let the compiler help us keep track + * of the real meaning of indices. + * + * The navigation classes provide a way to define an enumeration of the + * grid without resorting to inline integer arithmetic inside the other + * functions. An optimizing compiler should be able to generate equally + * fast code as hand-coded index calculations, when using these classes. + */ + +/** + * Index tuple in two-dimensional cornerpoint grid + * + * This structure represents the carrier of Cartesian coordinates. They + * should be thought of as an integral type, along the lines of complex + * numbers. + */ +struct Coord2D { + Coord2D (int i_, int j_) : m_i (i_), m_j (j_) { } + + int i() const { return m_i; } + int j() const { return m_j; } + + /** + * Compare two coordinates + */ + bool operator == (const Coord2D& rhs) const { + return (m_i == rhs.m_i) && (m_j == rhs.m_j); + } + +protected: + const int m_i; + const int m_j; + + friend std::ostream& operator << (std::ostream& s, const Coord2D& c); +}; + +/// Index tuple in three-dimensional cornerpoint grid. +struct Coord3D : public Coord2D { + Coord3D (int i_, int j_, int k_) + : Coord2D (i_, j_) + , m_k (k_) { + } + + int k() const { return m_k; } + +protected: + const int m_k; + + friend std::ostream& operator << (std::ostream& s, const Coord3D& c); +}; + +// forward declaration +template struct Side; + +/// Type-safe enumeration of axis directions. +struct Dir { + /// Towards the end of the axis with lesser numbers. + static const Dir DEC; // = 0 + + /// Towards the end of the axis with greater numbers. + static const Dir INC; // = 1 + + /// Number of possible directions + static const int COUNT = 2; + + /// Integer representation suitable for indexing in array + const int val; + + Dir (const Dir& rhs) : val (rhs.val) {} + bool operator == (const Dir& rhs) const { return val == rhs.val; } + + /// Opposite direction to this one + Dir opposite () const { return Dir (-(val - 1)); } + +protected: + /// Private constructor to avoid initialization outside domain + Dir (int i) : val (i) { } + + template friend struct Side; + + friend std::ostream& operator << (std::ostream& os, const Dir& d); +}; + +struct Dim1D { + static const int COUNT = 1; +}; + +/// Type-safe enumeration of axis dimensions +struct Dim2D : public Dim1D { + // two spatial directions + static const Dim2D X; // = 0 + static const Dim2D Y; // = 1 + + // number of dimensions + static const int COUNT = 2; + + const int val; + + Dim2D (const Dim2D& rhs) : val (rhs.val) { } + bool operator == (const Dim2D& rhs) const { return val == rhs.val; } + + /// Orthogonal dimension to this one + Dim2D orthogonal () const { return Dim2D (-(val - 1)); } + +protected: + Dim2D (int i) : val (i) { } + + friend struct Side ; + + friend std::ostream& operator << (std::ostream& os, const Dim2D& d); +}; + +/// Type-safe enumeration of axis dimensions in 3D +struct Dim3D : public Dim2D { + // added dimension in 3D + static const Dim3D Z; // = 2 + + // number of dimensions (shadows Dim2D) + static const int COUNT = 3; + + Dim3D (const Dim3D& rhs) : Dim2D (rhs.val) { } + bool operator == (const Dim2D& rhs) const { return val == rhs.val; } + + // allow X and Y to work in 3D too + Dim3D (const Dim2D& rhs) : Dim2D (rhs) {} + +protected: + Dim3D (int i) : Dim2D (i) { } + + friend struct Side ; +}; + +/** + * Value type that addresses sides in a n-dimensional grid cell. + * A side is identified by a dimensions and a direction in that + * dimension. The side will be located in that direction in that + * dimension from the center of the cell. + */ +template +struct Side { + Side (Dim dim_, Dir dir_) : m_dim (dim_), m_dir (dir_) { } + Side (const Side& rhs) : m_dim (rhs.m_dim), m_dir (rhs.m_dir) { } + + /** + * Numeric tag of an enumeration of the sides. The sides are enumerated + * in the same order the dimensions are specified in the grid; I, J then + * K, wih the decreasing direction first, then the increasing direction. + * + * @see UnstructuredGrid.cell_facetag + */ + int facetag () const { + return dim().val * Dir::COUNT + dir().val; + } + + /** + * Construct a side value from the facetag stored in the grid structure + */ + static Side from_tag (int tag); + + Dim dim() const { return m_dim; } + Dir dir() const { return m_dir; } + + /** + * Number of possible sides in an element + */ + static const int COUNT = Dim::COUNT * Dir::COUNT; + + /** + * Iterator for all possible sides + */ + static const Side* begin () { return &ALL[0]; } + static const Side* end () { return &ALL[COUNT]; } + + /** + * Comparison of two sides + */ + bool operator == (const Side & rhs) const { + return (m_dim == rhs.m_dim) && (m_dir == rhs.m_dir); + } + +protected: + const Dim m_dim; + const Dir m_dir; + + // fixed enumeration of all sides + static const Side ALL []; + + template + friend std::ostream& operator << (std::ostream& os, const Side& s); +}; + +// specializations for the dimensions we work with +typedef Side Side2D; +typedef Side Side3D; + +// forward declaration of stream operator present in library +std::ostream& operator << (std::ostream& os, const Side2D& s); +std::ostream& operator << (std::ostream& os, const Side3D& s); + +// standalone constants for sides that we use; we call them 'up' and +// 'down' so that U and D are mnemonics, in contrast to 'top' and 'bottom' +// where the 'b' would conflict with 'back'. +extern const Side3D UP; +extern const Side3D DOWN; + +/** + * Value type that addresses corners in a two-dimensional grid cell. + * A corner is identified by directions in both dimensions. + */ +struct Corn2D { + Corn2D (Dir i_, Dir j_) : m_i (i_), m_j (j_) { } + Corn2D (const Corn2D& rhs) : m_i (rhs.m_i), m_j (rhs.m_j) { } + + Dir i() const { return m_i; } + Dir j() const { return m_j; } + +protected: + const Dir m_i; + const Dir m_j; +}; + +/** + * Three-dimensional corner type. It inherits from the two-dimensional + * one since we can project a corner onto the two-dimensional surface. + */ +struct Corn3D : public Corn2D { + Corn3D (Dir i_, Dir j_, Dir k_) : Corn2D (i_, j_), m_k (k_) { } + Corn3D (const Corn3D& rhs) : Corn2D (rhs), m_k (rhs.m_k) { } + + /** + * Initialize a new corner where one dimension has been (optionally) + * pivoted in another direction. We use this to correct classifier + * information about a corner; by enumerating all the vertices of an + * element through the sides (pivoting a corner to the dimension in + * which its containing face is), we can figure out in which corner + * they belong in. + */ + Corn3D pivot (Dim3D dim, Dir dir) { + return Corn3D (dim == Dim3D::X ? dir : m_i, + dim == Dim3D::Y ? dir : m_j, + dim == Dim3D::Z ? dir : m_k); + } + + Dir k() const { return m_k; } + + /** + * Compare two corners + */ + bool operator == (const Corn3D& rhs) const { + return (m_i == rhs.m_i) && (m_j == rhs.m_j) && (m_k == rhs.m_k); + } + +protected: + const Dir m_k; + + friend std::ostream& operator << (std::ostream& os, const Corn3D& c); +}; + +/** + * Navigate a Cartesian grid in a structured way so that clearly defined + * mapping between the enumeration index and the coordinate. + */ +struct Cart2D { + // number of cells in each direction + const int ni; + const int nj; + + // initialize from the size of the grid + Cart2D (int ni_, int nj_) : ni (ni_), nj (nj_) { } + + // use these two value types for structured coordinates and + // flattened indices, respectively + typedef Coord2D coord_t; + typedef int elem_t; + typedef int node_t; + typedef int face_t; + + /// Value used to indicate that a reference is not to a valid element + static const int NO_ELEM; // = -1 + static const int NO_FACE; // = -1 + static const int NO_NODE; // = -1 + + /// Number of (possible) elements in the grid + int num_elems () const { + return ni * nj; + } + + /// Cartesian (flattened) index for a coordinate + elem_t cart_ndx (const coord_t& coord) const { + return coord.j() * ni + coord.i(); + } + + /// Cartesian coordinate for a (flattened) index + coord_t coord (const elem_t& cart_ndx) const { + const std::div_t strip = std::div (cart_ndx, ni); + const int i = strip.rem; + const int j = strip.quot; + return coord_t (i, j); + } + + /** + * As each element has points on both sides (in both dimensions), there + * is an extra row and column of points compared to elements. + */ + int num_nodes () const { + return (ni + 1) * (nj + 1); + } + + node_t node_ndx (const coord_t& coord, const Corn2D& corn) { + return (coord.j() + corn.j().val) * (ni + 1) + (coord.i() + corn.i().val); + } + + /** + * Each column has one more faces than there are rows, and each row + * has one more face than there are columns, due to the boundary. + */ + int num_faces () const { + // there are ni+1 faces oriented in j-direction in each column, + // for a total of (ni+1)*nj, and nj+1 faces in i-direction in + // each row, for a total of (nj+1)*ni. + return (ni + 1) * nj + ni * (nj + 1); + } + + /** + * Translate element coordinate plus relative side of this into a + * Cartesian index of the face. + */ + face_t face_ndx (const coord_t& coord, const Side2D& side) { + // flags that code 1 (include) or 0 (don't) for each of the tests + const int dirv = side.dir().val; + const int idim = side.dim() == Dim2D::X ? 1 : 0; + const int idir = idim ? dirv : 0; + const int jdir = idim ? 0 : dirv; + // there is a left side and a lower side face for each element in a + // column, plus one face at the top of each column; 2*ni+1. the j- + // coordinate plus one extra if we are selecting the right face, tells + // us which column which is to the left of or "above" the face. + // if the side is in the i-direction to the center of the element, then + // the face is aligned with the j axis, so we skip idim*ni i-aligned + // faces before it. + // finally, determine the row by using the i-coordinate, and i-direction + // to determine whether it is the lower or upper face (if applicable). + return (coord.j() + jdir) * (2 * ni + 1) + (idim * ni) + (coord.i() + idir); + } +}; + +/** + * Navigate a three-dimensional grid. + * + * In this module, we only need this to get the structured index of + * each three-dimensional element, in order to know into which column + * we should assign this element. However, we keep the design in the + * same manner as the two-dimensional case. + */ +struct Cart3D { + // number of cells in each direction + const int ni; + const int nj; + const int nk; + + /// Initialize POD from an existing (3D) grid + Cart3D (const UnstructuredGrid& g) + : ni (g.cartdims [0]) + , nj (g.cartdims [1]) + , nk (g.cartdims [2]) { } + + /// Project grid into a surface + Cart2D project () const { + return Cart2D (ni, nj); + } + + // use these two value types for structured coordinates and + // flattened indices, respectively + typedef Coord3D coord_t; + typedef int elem_t; + + /// Deconstruct Cartesian index into coordinates + coord_t coord (const elem_t& cart_ndx) const { + // the i-index moves fastest, as this is Fortran-indexing + const div_t strip = div (cart_ndx, ni); + const int i = strip.rem; + const div_t plane = div (strip.quot, nj); + const int j = plane.rem; + const int k = plane.quot; + return coord_t (i, j, k); + } +}; + +#endif // OPM_VERTEQ_NAV_HPP_INCLUDED diff --git a/opm/grid/verteq/topsurf.cpp b/opm/grid/verteq/topsurf.cpp new file mode 100644 index 0000000000..3a371e1525 --- /dev/null +++ b/opm/grid/verteq/topsurf.cpp @@ -0,0 +1,826 @@ +// Copyright (C) 2013 Uni Research AS +// This file is licensed under the GNU General Public License v3.0 + +#include +#include +#include +#include // rlw_int +#include // compute_geometry +#include // ios_all_saver +#include // min, max +#include // INT_MIN, INT_MAX +#include // div +#include // ostream +#include +#include // unique_ptr +#include // accumulate, iota +#include +#include // pair + +using namespace boost; +using namespace Opm; +using namespace std; + +/// Helper routine to print of map +template +void dump_map (ostream& os, const map& m) { + for (typename map::const_iterator it = m.begin(); it != m.end(); ++it) { + boost::io::ios_all_saver state (os); + os << it->first << ":\t"; + state.restore (); + os << it->second << endl; + } +} + +/** + * @brief Process to extract the top surface from a structured grid. + * + * This object encapsulates a procedure with variables shared amongst + * several sub-procedures (like in Pascal). These objects are not + * supposed to linger on afterwards. + */ +struct TopSurfBuilder { + // source grid from which we get the input data + const UnstructuredGrid& fine_grid; + + // target grid we are constructing + TopSurf& ts; + + // number of grid dimensions in each direction of the plane + Cart3D three_d; + + // dimensions needed to create a two-dimensional projection + // of the top surface + Cart2D two_d; + + // map from a two-dimensional Cartesian coordinate to the final + // id of an active element in the grid, or NO_ELEM if nothing is assigned + // this vector is first valid after create_elements() have been done + vector elms; + + // map from a two-dimensional Cartesian node coordinate to the final + // id of active nodes in the grid, or NO_NODE if nothing is assigned. + // this vector is first valid after create_nodes() have been done + vector nodes; + + // map from a two-dimensional Cartesion face coordinate to the final + // id of active faces in the grid, or NO_FACE if nothing is assigned. + // this vector is first valid after create_faces() have been done. + vector faces; + + // logical Cartesian indices for items in the fine grid. we need our + // own copy of this since not all grids provide it + vector fine_global; + + TopSurfBuilder (const UnstructuredGrid& from, TopSurf& into) + // link to the fine grid for the duration of the construction + : fine_grid (from) + + // allocate memory for the grid. it is initially empty + , ts (into) + + // extract dimensions from the source grid + , three_d (fine_grid) + , two_d (three_d.project ()) + , fine_global (from.number_of_cells, 0) { + + // check that the fine grid contains structured information; + // this is essential to mapping cells to columns + const int prod = std::accumulate(&fine_grid.cartdims[0], + &fine_grid.cartdims[fine_grid.dimensions], + 1, + std::multiplies()); + if (!prod) { + throw OPM_EXC ("Find grid is not (logically) structured"); + } + + // some cartesian grids (most notably those generated with + // create_grid_cart{2,3}d) have no global cell + if (!fine_grid.global_cell) { + std::iota (fine_global.begin(), fine_global.end(), 0); + } + else { + std::copy (fine_grid.global_cell, + fine_grid.global_cell + fine_grid.number_of_cells, + fine_global.begin ()); + } + + // create frame of the new top surface + create_dimensions (); + + // identify active columns in the grid + create_elements (); + + // identify active points in the grid + create_nodes (); + + // identify active faces in the grid + create_faces (); + + // cache fine block and column metrics + create_heights (); + } + + // various stages of the build process, supposed to be called in + // this order. (I have separated them into separate procedures to + // make it more obvious what parts that needs to be shared between + // them) +private: + void create_dimensions () { + // we are going to create two-dimensional grid + ts.dimensions = 2; + ts.cartdims[0] = two_d.ni; + ts.cartdims[1] = two_d.nj; + ts.cartdims[2] = 1; + } + + void create_elements() { + // statistics of the deepest and highest active k-index of + // each column in the grid. to know each index into the column, + // we only need to know the deepest k and the count; the highest + // is sampled to do consistency checks afterwards + const int num_cols = two_d.num_elems (); + + // assume initially that there are no active elements in each column + vector act_cnt (num_cols, 0); + ts.number_of_cells = 0; + + // initialize these to values that are surely out of range, so that + // the first invocation of min or max always set the value. we use + // this to detect whether anything was written later on. since the + // numbering of the grid starts at the top, then the deepest cell + // has the *largest* k-index, thus we need a value smaller than all + vector deep_k (num_cols, INT_MIN); + vector high_k (num_cols, INT_MAX); + + // loop once through the fine grid to gather statistics of the + // size of the surface so we know what to allocate + for (int fine_elem = 0; fine_elem != fine_grid.number_of_cells; ++fine_elem) { + // get the cartesian index for this cell; this is the cell + // number in a grid that also includes the inactive cells + const Cart3D::elem_t cart_ndx = fine_global [fine_elem]; + + // deconstruct the cartesian index into (i,j,k) constituents; + // the i-index moves fastest, as this is Fortran-indexing + const Coord3D ijk = three_d.coord (cart_ndx); + + // figure out which column this item belongs to (in 2D) + const Cart2D::elem_t col = two_d.cart_ndx (ijk); + + // update the statistics for this column; 'deepest' is the largest + // k-index seen so far, 'highest' is the smallest (ehm) + deep_k[col] = max (deep_k[col], ijk.k()); + high_k[col] = min (high_k[col], ijk.k()); + + // we have seen an element in this column; it becomes active. only + // columns with active cells will get active elements in the surface + // grid. if this cell wasn't marked as active before we can check it + // off now + if (!act_cnt[col]++) { + ts.number_of_cells++; + } + } + + // check that we have a continuous range of elements in each column; + // this must be the case to assume that the entire column can be merged + for (int col = 0; col < num_cols; ++col) { + if (act_cnt[col]) { + if (high_k[col] + act_cnt[col] - 1 != deep_k[col]) { + const Coord2D coord = two_d.coord (col); + throw OPM_EXC ("Non-continuous column at (%d, %d)", coord.i(), coord.j()); + } + } + } + + // allocate memory needed to hold the elements in the grid structure + // if we throw an exception at some point, the destructor of the TopSurf + // will take care of deallocating this memory for us + ts.global_cell = new int [ts.number_of_cells]; + ts.col_cellpos = new int [ts.number_of_cells+1]; + + // we haven't filled any columns yet, so this is a sensible init value + ts.max_vert_res = 0; + + // there is no elements before the first column, so this number is + // always zero. it is written to the array to maintain a straight code + // path in calculations. + ts.col_cellpos[0] = 0; + + // now we know enough to start assigning ids to active columns.if + // memory is a real shortage, we could reuse the act_cnt array for this. + elms.resize (num_cols, Cart2D::NO_ELEM); + + // loop through the grid and assign an id for all columns that have + // active elements + int elem_id = 0; + for (int col = 0; col < num_cols; ++col) { + if (act_cnt[col]) { + elms[col] = elem_id; + + // dual pointer that maps the other way; what is the structured + // index of this element (flattened into an integer) + ts.global_cell[elem_id] = col; + + // note the number of elements there now are before the next column; + // in addition to all the previous ones, our elements are now added + ts.col_cellpos[elem_id+1] = ts.col_cellpos[elem_id] + act_cnt[col]; + + // update the largest number of these seen so far + ts.max_vert_res = max (ts.max_vert_res, act_cnt[col]); + + // only increment this if we found an active column, elem_id <= col + elem_id++; + } + } + + // now write indices from the fine grid into the column map of the surface + // we end up with a list of element that are in each column + ts.col_cells = new int [fine_grid.number_of_cells]; + ts.fine_col = new int [fine_grid.number_of_cells]; + for (int cell = 0; cell < fine_grid.number_of_cells; ++cell) { + // get the Cartesian index for this element + const Cart3D::elem_t cart_ndx = fine_global[cell]; + const Coord3D ijk = three_d.coord (cart_ndx); + + // get the id of the column in which this element now belongs + const Cart2D::elem_t col = two_d.cart_ndx (ijk); + const int elem_id = elms[col]; + + // start of the list of elements for this particular column + const int segment = ts.col_cellpos[elem_id]; + + // since there is supposed to be a continuous range of elements in + // each column, we can calculate the relative position in the list + // based on the k part of the coordinate. + const int offset = ijk.k() - high_k[col]; + + // write the fine grid cell number in the column list; since we + // have calculated the position based on depth, the list will be + // sorted downwards up when we are done + ts.col_cells[segment + offset] = cell; + + // reverse mapping; allows us to quickly figure out the corresponding + // column of a location (for instance for a well) + ts.fine_col[cell] = elem_id; + } + + // these members will be filled by computeGeometry, but needs valid + // memory to work with + ts.cell_volumes = new double [ts.number_of_cells]; + ts.cell_centroids = new double [ts.dimensions * ts.number_of_cells]; + } + + void create_nodes () { + // construct a dual Cartesian grid consisting of the points + const int num_nodes = two_d.num_nodes (); + + // vectors which will hold the coordinates for each active point. + // at first we sum all the points, then we divide by the count to + // get the average. as long as the count is zero, there is no + // registered active point at this location + vector x (num_nodes, 0.); + vector y (num_nodes, 0.); + vector cnt (num_nodes, 0); + + // number of nodes needed in the top surface + int active_nodes = 0; + + // tag of the top side in a cell; we're looking for this + const int top_tag = Side3D (Dim3D::Z, Dir::DEC).facetag (); + + // initial corner value. this could really be anything, since + // we expect all the fields to be overwritten. + const Corn3D blank (Dir::DEC, Dir::DEC, Dir::DEC); + + // this map holds the classification of each node locally for the + // element being currently processed. we reuse the map to avoid + // needless memory allocations. + typedef map cls_t; + cls_t classifier; + + // loop through all active cells in the top surface + for (int col = 0; col < ts.number_of_cells; ++col) { + // get the highest element in this column; since we have them + // sorted by k-index this should be the first item in the + // extended column info + const Cart3D::elem_t top_cell_glob_id = ts.col_cells [ts.col_cellpos[col]]; + + // start afresh whenever we start working on a new element + classifier.clear (); + int top_face_glob_id = Cart2D::NO_FACE; + + // loop through all the faces of the top element + for (int face_pos = fine_grid.cell_facepos[top_cell_glob_id]; + face_pos != fine_grid.cell_facepos[top_cell_glob_id+1]; + ++face_pos) { + + // get the (normal) dimension and direction of this face + const int this_tag = fine_grid.cell_facetag[face_pos]; + Side3D s = Side3D::from_tag (this_tag); + + // identifier of the face, which is the index in the next arary + const int face_glob_id = fine_grid.cell_faces[face_pos]; + + // remember it if we've found the top face + if (this_tag == top_tag) { + if (top_face_glob_id != Cart2D::NO_FACE) { + throw OPM_EXC ("More than one top face in element %d", top_cell_glob_id); + } + top_face_glob_id = face_glob_id; + } + + // loop through all nodes in this face, adding them to the + // classifier. when we are through with all the faces, we have + // found in which corner a node is, defined by a direction in + // each of the three dimensions + for (int node_pos = fine_grid.face_nodepos[face_glob_id]; + node_pos != fine_grid.face_nodepos[face_glob_id+1]; + ++node_pos) { + const int node_glob_id = fine_grid.face_nodes[node_pos]; + + // locate pointer to data record ("iterator" in stl parlance) + // for this node, if it is already there. otherwise, just start + // out with some blank data (which eventually will get overwritten) + cls_t::iterator ptr = classifier.find (node_glob_id); + Corn3D prev (ptr == classifier.end () ? blank : ptr->second); + + // update the dimension in which this face is pointing + if (ptr != classifier.end ()) { + classifier.erase (ptr); + } + const Corn3D upd_corn = prev.pivot (s.dim(), s.dir()); + classifier.insert (make_pair (node_glob_id, upd_corn)); + } + + // after this loop, we have a map of each node local to the element, + // classified into in which corner it is located (it cannot be in + // both directions in the same dimension -- then it would have to + // belong to two opposite faces, unless the grid is degenerated) + } + /* + cerr << "elem: " << three_d.coord(top_cell_glob_id) << ':' << endl; + dump_map (cerr, classifier); + */ + + // cannot handle degenerate grids without top face properly + if (top_face_glob_id == Cart2D::NO_FACE) { + throw OPM_EXC ("No top face in cell %d", top_cell_glob_id); + } + + // get the Cartesian ij coordinate of this cell + const Cart2D::elem_t top_cell_cart_ndx = ts.global_cell [top_cell_glob_id]; + const Coord2D ij = two_d.coord (top_cell_cart_ndx); + + // loop through all the nodes of the top face, and write their position + // into the corresponding two-d node. this has to be done separately + // after we have classified *all* the nodes of the element, in order for + // the corner values to be set correctly, i.e. we cannot merge this into + // the loop above. + for (int node_pos = fine_grid.face_nodepos[top_face_glob_id]; + node_pos != fine_grid.face_nodepos[top_face_glob_id+1]; + ++node_pos) { + const int node_glob_id = fine_grid.face_nodes[node_pos]; + + // get which corner this node has; this returns a three-dimensional + // corner, but by using the base class part of it we automatically + // project it to a flat surface + cls_t::iterator ptr = classifier.find (node_glob_id); + const Corn3D corn (ptr->second); + + // get the structured index for this particular corner + const Cart2D::node_t cart_node = two_d.node_ndx(ij, corn); + + // add these coordinates to the average position for this junction + // if we activate a corner, then add it to the total count + x[cart_node] += fine_grid.node_coordinates[Dim3D::COUNT*node_glob_id+0]; + y[cart_node] += fine_grid.node_coordinates[Dim3D::COUNT*node_glob_id+1]; + if (!cnt[cart_node]++) { + ++active_nodes; + } + } + } + + // after this loop we the accumulated coordinates for each of the + // corners that are part of active elements (the nodes that are + // needed in the top surface) + + // assign identifiers and find average coordinate for each point + ts.number_of_nodes = active_nodes; + ts.node_coordinates = new double [active_nodes * Dim2D::COUNT]; + nodes.resize (num_nodes, Cart2D::NO_NODE); + int next_node_id = 0; + for (int cart_node = 0; cart_node != num_nodes; ++cart_node) { + if (cnt[cart_node]) { + nodes[cart_node] = next_node_id; + const int start = Dim2D::COUNT * next_node_id; + ts.node_coordinates[start+0] = x[cart_node] / cnt[cart_node]; + ts.node_coordinates[start+1] = y[cart_node] / cnt[cart_node]; + ++next_node_id; + } + } + + // dump node topology to console + /* + for (int cart_node_ndx = 0; cart_node_ndx != num_nodes; ++cart_node_ndx) { + const int glob_node_id = nodes[cart_node_ndx]; + if (glob_node_id != Cart2D::NO_NODE) { + cerr << "node " << nodes[cart_node_ndx] << ": (" + << ts.node_coordinates[2*glob_node_id+0] << ',' + << ts.node_coordinates[2*glob_node_id+1] << ')' << endl; + } + } + */ + + // TODO: check for degeneracy by comparing each node's coordinates + // with those in the opposite direction in both dimensions (separately) + } + + void create_faces () { + // number of possible (but not necessarily active) faces + const int num_faces = two_d.num_faces (); + + // assign identifiers into this array. start out with the value + // NO_FACE which means that unless we write in an id, the face + // is not active + faces.resize (num_faces, Cart2D::NO_FACE); + + // a face will be referenced from two elements (apart from boundary), + // denoted the "primary" and "secondary" neighbours. the nodes in a + // face needs to be specified so that the normal to the directed LINE_NODES + // points towards the element center of the *primary* neighbour. + + // we use the convention that every face that are in the J-direction + // are directed from decreasing direction to increasing direction, + // whereas every face that are in the I-direction are directed the + // opposite way, from increasing to decreasing. this way, an element + // is primary neighbour for a face if the face is on side which is + // classified as decreasing, relative to the center, and secondary + // if the face is in the increasing direction of whatever axis. + + // I+ (I+,J+) (I+,J-) + // o <---- o o <---- o + // | | | | + // J+ | | J- | | + // v v v v + // o <---- o o <---- o + // I- (I-,J+) (I-,J-) + + // TODO: Possible to instead create an S-shaped traversal with + // different directions every odd/even column, so the faces ends + // up naturally in clockwise directions around the element? + + // allocate two arrays that will hold the global ids of the two + // elements on each side, or NO_ELEM if there is no active element + // on that side. if both sides are NO_ELEM then the face can be + // optimized away from the resulting grid. + vector pri_elem (num_faces, Cart2D::NO_ELEM); + vector sec_elem (num_faces, Cart2D::NO_ELEM); + + vector src (num_faces, Cart2D::NO_NODE); + vector dst (num_faces, Cart2D::NO_NODE); + + // keep count on how many faces that actually have at least one + // active neighbouring element + int active_faces = 0; + + // loop through all the active elements in the grid, and write + // their identifier in the neighbour array for their faces + for (int j = 0; j != two_d.nj; ++j) { + for (int i = 0; i != two_d.ni; ++i) { + // where are we in the grid? + const Coord2D coord (i, j); + const int col = two_d.cart_ndx (coord); + + // check if this element is active; is there an id assignment? + const int elem_glob_id = elms[col]; + if (elem_glob_id != Cart2D::NO_ELEM) { + + // loop through all sides of this element; by assigning identities + // to the faces in this manner, the faces around each element are + // relatively local to eachother in the array. + for (const Side2D* s = Side2D::begin(); s != Side2D::end(); ++s) { + // cartesian index of this face, i.e. index just depending on the + // extent of the grid, not whether face is active or not + const int cart_face = two_d.face_ndx (coord, *s); + + // select primary or secondary neighbour collection based on + // the direction of the face relative to the center of the + // element, see discussion above. if the vector from the center + // to the face points in the INC direction (i.e. this is an INC + // face), then that vector is aligned with the right-normal of + // the face, and this node is the primary. + const bool is_primary = s->dir () == Dir::INC; + vector & neighbour = is_primary ? pri_elem : sec_elem; + vector & other = is_primary ? sec_elem : pri_elem; + + // put this identifier in there + if (neighbour[cart_face] != Cart2D::NO_ELEM) { + throw OPM_EXC ("Duplicate neighbour assignment in column (%d,%d)", + coord.i(), coord.j()); + } + neighbour[cart_face] = elem_glob_id; + + // if this was the first time we assigned a neighbour to the + // face, then count it as active and assign an identity + if (other[cart_face] == Cart2D::NO_ELEM) { + faces[cart_face] = active_faces++; + + // faces that are in the I-dimension (I- and I+) starts at the DEC + // direction in the J-dimension, where as for the faces in the J- + // dimension it is opposite + const Dir src_dir = s->dim () == Dim2D::X ? Dir::DEC : Dir::INC; + const Dir dst_dir = src_dir.opposite (); + + // use the direction of the side for its dimension, as both the corners + // of the side will be here, and use the two other directions for the + // remaining dimension. the trick is to know in which corner the face + // should start, see above. + const Corn2D src_corn (s->dim () == Dim2D::X ? s->dir () : src_dir, + s->dim () == Dim2D::X ? src_dir : s->dir ()); + const Corn2D dst_corn (s->dim () == Dim2D::X ? s->dir () : dst_dir, + s->dim () == Dim2D::X ? dst_dir : s->dir ()); + + // get the identity of the two corners, and take note of these + const int src_cart_ndx = two_d.node_ndx (coord, src_corn); + const int dst_cart_ndx = two_d.node_ndx (coord, dst_corn); + src[cart_face] = nodes[src_cart_ndx]; + dst[cart_face] = nodes[dst_cart_ndx]; + } + } + } + } + } + // after this loop we have a map of all possible faces, and know + // how many of these are active. + + // dump face topology to console + /* + for (int cart_face_ndx = 0; cart_face_ndx != num_faces; ++cart_face_ndx) { + const int face_glob_id = faces[cart_face_ndx]; + if (face_glob_id != Cart2D::NO_FACE) { + cerr << "face " << face_glob_id << ": " << endl; + cerr << "\t" << "elements: " << pri_elem[cart_face_ndx] << "," + << sec_elem[cart_face_ndx] << endl; + cerr << "\t" << "nodes: " << src[cart_face_ndx] << "," + << dst[cart_face_ndx] << endl; + } + } + */ + + // each cell has 4 sides in 2D; we assume no degenerate sides + const int QUAD_SIDES = Dim2D::COUNT * Dir::COUNT; + const int num_sides = QUAD_SIDES * active_faces; + + // nodes in faces; each face in 2D is only 1D, so simplices always + // have only two nodes (a LINE_NODES, with corners in each direction) + const int LINE_NODES = Dim1D::COUNT * Dir::COUNT; + const int num_corns = LINE_NODES * active_faces; + + // number of element neighbours for each face. this is always 2, + // the reason for not using the number is to make it searchable + const int NEIGHBOURS = 2; + + // allocate memory; this will be freed in the TopSurf destructor + ts.face_nodes = new int [num_corns]; + ts.face_nodepos = new int [active_faces + 1]; + ts.face_cells = new int [NEIGHBOURS * active_faces]; + ts.cell_faces = new int [num_sides]; + ts.cell_facepos = new int [ts.number_of_cells + 1]; + ts.cell_facetag = new int [num_sides]; + ts.number_of_faces = active_faces; + + // we need to allocate memory for computeGeometry to fill + ts.face_centroids = new double [Dim2D::COUNT * active_faces]; + ts.face_areas = new double [active_faces]; + ts.face_normals = new double [Dim2D::COUNT * active_faces]; + + // write the internal data structures to UnstructuredGrid representation + + // face <-> node topology + for (int cart_face = 0; cart_face != num_faces; ++cart_face) { + const int face_glob_id = faces[cart_face]; + if (face_glob_id != Cart2D::NO_FACE) { + // since each face has exactly two coordinates, we can easily + // calculate the position based only on the face number + const int start_pos = LINE_NODES * face_glob_id; + ts.face_nodepos[face_glob_id] = start_pos; + ts.face_nodes[start_pos + 0] = src[cart_face]; + ts.face_nodes[start_pos + 1] = dst[cart_face]; + + // TODO: If a vertical fault displaces two column so that there + // is no longer connection between them, they will be reconnected + // here. This condition can be detected by comparing the top and + // bottom surface. + + // neighbours should already be stored in the right orientation + ts.face_cells[NEIGHBOURS * face_glob_id + 0] = pri_elem[cart_face]; + ts.face_cells[NEIGHBOURS * face_glob_id + 1] = sec_elem[cart_face]; + } + } + ts.face_nodepos[ts.number_of_faces] = LINE_NODES * ts.number_of_faces; + + // cell <-> face topology + for (int j = 0; j != two_d.nj; ++j) { + for (int i = 0; i != two_d.ni; ++i) { + // get various indices for this element + const Coord2D coord (i, j); + const int cart_elem = two_d.cart_ndx (coord); + const int elem_glob_id = elms[cart_elem]; + + // each element is assumed to be a quad, so we can calculate the + // number of accumulated sides based on the absolute id + const int start_pos = QUAD_SIDES * elem_glob_id; + ts.cell_facepos[elem_glob_id] = start_pos; + + // write all faces for this element + for (const Side2D* s = Side2D::begin(); s != Side2D::end(); ++s) { + // get the global id of this face + const int face_cart_ndx = two_d.face_ndx (coord, *s); + const int face_glob_id = faces[face_cart_ndx]; + + // the face tag can also serve as an offset into a regular element + const int ofs = s->facetag (); + ts.cell_faces[start_pos + ofs] = face_glob_id; + ts.cell_facetag[start_pos + ofs] = ofs; + } + } + } + ts.cell_facepos[ts.number_of_cells] = QUAD_SIDES * ts.number_of_cells; + } + + /** + * Specific face number of a given side of an element. + * + * @param glob_elem_id Element index in the fine grid. + * @param s Side to locate + * @return Index of the face of the element which is this side + * + * @see Opm::UP, Opm::DOWN + */ + int find_face (int glob_elem_id, const Side3D& s) { + // this is the tag we are looking for + const int target_tag = s.facetag (); + + // this is the matrix we are looking in + const rlw_int cell_facetag = grid_cell_facetag (fine_grid); + + // we are returning values from this matrix + const rlw_int cell_faces = grid_cell_faces (fine_grid); + + // loop through all faces for this element; face_ndx is the local + // index amongst faces for just this one element. + for (int local_face = 0; + local_face < cell_facetag.size (glob_elem_id); + ++local_face) { + + // if we found a match, then return this; don't look any more + if (cell_facetag[glob_elem_id][local_face] == target_tag) { + + // return the (global) index of the face, not the tag! + return cell_faces[glob_elem_id][local_face]; + } + } + + // in a structured grid we expect to find every face + throw OPM_EXC ("Element %d does not have face #%d", glob_elem_id, target_tag); + } + + /** + * Get absolute elevation (z-coordinate) of a face. This uses the + * elevation at the centroid as representative of the entire face. + * + * @param glob_elem_id Element index in the fine grid. + * @param s Side to locate. + * @return Elevation for the midpoint of this face. + */ + double find_zcoord (int glob_elem_id, const Side3D& s) { + // find the desired face for this element + const int face_ndx = find_face (glob_elem_id, s); + + // get the z-coordinate for it + const int z_ndx = face_ndx * Dim3D::COUNT + Dim3D::Z.val; + const double z = fine_grid.face_centroids[z_ndx]; + return z; + } + + /** + * Height of a particular element. + * + * @param glob_elem_id Element index in the fine grid. + * @return Difference between center of top and bottom face. + */ + double find_height (int glob_elem_id) { + // get the z-coordinate for each the top and bottom face for this element + const double up_z = find_zcoord (glob_elem_id, UP); + const double down_z = find_zcoord (glob_elem_id, DOWN); + + // the side that is down should have the z coordinate with highest magnitude + const double height = down_z - up_z; + return height; + } + + void create_heights () { + // allocate memory to hold the heights + ts.dz = new double [fine_grid.number_of_cells]; + ts.h = new double [fine_grid.number_of_cells]; + ts.z0 = new double [ts.number_of_cells]; + ts.h_tot = new double [ts.number_of_cells]; + + // view that lets us treat it as a matrix + const rlw_int blk_id (ts.number_of_cells, ts.col_cellpos, ts.col_cells); + const rlw_double dz (ts.number_of_cells, ts.col_cellpos, ts.dz); + const rlw_double h (ts.number_of_cells, ts.col_cellpos, ts.h); + + // find all measures per column + for (int col = 0; col < blk_id.cols (); ++col) { + // reference height for this column (if there is any elements) + if (blk_id.size (col)) { + const int top_ndx = blk_id[col][0]; + ts.z0[col] = find_zcoord (top_ndx, UP); + } + + // reset height for each column + double accum = 0.; + + // height of each element in the column element + double* const dz_col = dz[col]; + double* const h_col = h[col]; + for (int col_elem = 0; col_elem < blk_id.size (col); ++col_elem) { + h_col[col_elem] = accum; + accum += dz_col[col_elem] = find_height (blk_id[col][col_elem]); + } + + // store total accumulated height at the end for each column + ts.h_tot[col] = accum; + } + } +}; + +TopSurf* +TopSurf::create (const UnstructuredGrid& fine_grid) { + unique_ptr ts (new TopSurf); + + // outsource the entire construction to a builder object + TopSurfBuilder (fine_grid, *(ts.get ())); + compute_geometry (ts.get ()); + + // client owns pointer to constructed grid from this point + return ts.release (); +} + +TopSurf::TopSurf () + : col_cells (0) + , col_cellpos (0) + , fine_col (0) + , dz (0) + , z0 (0) + , h_tot (0) { + // zero initialize all members that come from UnstructuredGrid + // since that struct is a C struct, it doesn't have a ctor + dimensions = 0; + number_of_cells = 0; + number_of_faces = 0; + number_of_nodes = 0; + face_nodes = 0; + face_nodepos = 0; + face_cells = 0; + cell_faces = 0; + cell_facepos = 0; + node_coordinates = 0; + face_centroids = 0; + face_areas = 0; + face_normals = 0; + cell_centroids = 0; + cell_volumes = 0; + global_cell = 0; + cartdims[0] = 0; + cartdims[1] = 0; + cartdims[2] = 0; + cell_facetag = 0; +} + +TopSurf::~TopSurf () { + // deallocate memory that may have been created. if the dtor is + // called from throwing an exception, the members should be zero + // initialized, so it's OK to send them to delete. + delete [] face_nodes; + delete [] face_nodepos; + delete [] face_cells; + delete [] cell_faces; + delete [] cell_facepos; + delete [] node_coordinates; + delete [] face_centroids; + delete [] face_areas; + delete [] face_normals; + delete [] cell_volumes; + delete [] global_cell; + delete [] cell_facetag; + // these are the extra members that are TopSurf specific + delete [] col_cells; + delete [] col_cellpos; + delete [] fine_col; + delete [] dz; + delete [] h; + delete [] z0; + delete [] h_tot; +} diff --git a/opm/grid/verteq/topsurf.hpp b/opm/grid/verteq/topsurf.hpp new file mode 100644 index 0000000000..70f0cc11bd --- /dev/null +++ b/opm/grid/verteq/topsurf.hpp @@ -0,0 +1,177 @@ +#ifndef OPM_VERTEQ_TOPSURF_HPP_INCLUDED +#define OPM_VERTEQ_TOPSURF_HPP_INCLUDED + +// Copyright (C) 2013 Uni Research AS +// This file is licensed under the GNU General Public License v3.0 + +#ifndef OPM_GRID_HEADER_INCLUDED +#include +#endif + +namespace Opm { + +/** + * Two-dimensional top surface of a full, three-dimensional grid. + * + * This grid is set up such that each cell is an upscaling of all the cells + * in each column in the full grid. It also contains a mean to map results + * from this grid back to the full grid. + * + * The full grid is also referred to as the fine grid, and this grid as the + * coarse grid, or upscaled, grid. + * + * Note: Do NOT call destroy_grid () when done with this structure; it will + * only clean up half of it. Wrap it is a smart pointer that calls the + * destructor. + */ +struct TopSurf : public UnstructuredGrid { + virtual ~TopSurf (); + + /** + * Indices of the columns' underlaying cells in the full grid. + * + * Consecutive indices from the _fine_ grid, not this one, for each of the + * columns, i.e. cells in _this_ grid, as one flat list. + * + * The values are sorted in z-order, starting from the top and moving + * downwards to the bottom. This is useful because you can keep a running + * counter for the depth, filling items as you go. (For this to be really + * useful, the original grid should be reordered so that cells in the + * z-direction are closer). + * + * Use this field together with the col_cellpos to iterate through a column + * in the fine grid. + * + * @example + * @code{.cpp} + * TopSurf* ts = ...; + * rlw_int col_cells (ts->number_of_cells, ts->col_cellpos, ts->col_cells); + * for (int col = 0; col < col_cells.cols(); ++col) { + * for (int block = 0; block < col_cells.size (col); ++block) { + * ... col_cells[col][block] ... + * } + * } + * @endcode + * + * @see TopSurf::column, TopSurf::col_cellpos + */ + int* col_cells; + + /** + * Number of cells in the columns preceeding each one. + * + * For each column c, the number col_cellpos[c] is the number of cells in + * the _full_ grid that belongs to the columns 0..(c-1). + * + * This arrangement means that col_cellpos[c] is the index into col_cells + * of the first fine cell, whereas col_cellpos[c+1] is the index into + * col_cells of the last fine cell. + * + * @see TopSurf::column, TopSurf::col_cellpos + */ + int* col_cellpos; + + /** + * Maximum vertical resolution, in blocks. + * + * This holds the largest number of blocks there is in any column in the + * fine grid, i.e. max_vert_res >= col_cellpos[i+1] - col_cellpos[i], for + * any i in [0,number_of_cells-1]. Use this measure to allocate sufficient + * space for temporary storage that holds column data. + */ + int max_vert_res; + + /** + * Mapping from underlaying fine grid into top surface grid. + * + * For each element e in the fine grid, cell_col[e] is the index of the + * column/cell in the top surface. The number of cells in the fine grid + * can be found in col_cellpos[number_of_cells+1]. + * + * Note: The indices in this array is NOT cell indices of this grid, + * but rather of the underlaying grid. Instead, it are the values that + * are stored in this array which are element identities. + * + * @see TopSurf::col_cellpos, TopSurf::col_cells + */ + int* fine_col; + + /** + * Height in each fine grid block, setup in columns. + * + * For each column, there is a consecutive list of height for each fine + * grid block in that column. This array has the same format as the + * col_cells run-length matrix, and the heights given here correspond to + * the indices in that matrix. + * + * The height of a block is defined as the z-coordinate difference + * between the centroid of the top face and the centroid of the bottom + * face. + * + * @see TopSurf::col_cells + */ + double* dz; + + /** + * Reference height for each column. + * + * This is a flat array with number_of_elements items. + * + * The reference height of a column is defined as the z-coordinate of + * the centroid of the top face of the upper block in the column. From + * these values and (a subset of) the values in face_centroids it is + * possible to recreate the 2.5D surface of the top. + */ + double* z0; + + /** + * Height from top of column down to each fine grid block. + * + * This is the accumulated sum of all dz preceeding this block, in each + * column. The first entry is thus always zero. + */ + double* h; + + /** + * Accumulated height of all blocks in each column. + * + * This is a flat array with number_of_elements items. + * + * Sum of all the heights of the blocks in each column. Since a gap between + * blocks is a violation of the vertical equilibrium assumption, this value + * should be the z-coordinate of the bottom face of the last block, up to + * numerical rounding errors. + * + * @see TopSurf::dz, TopSurf::z0 + */ + double* h_tot; + + /** + * Create an upscaled grid based on a full, three-dimensional grid. + * + * @param fine Grid that should be upscaled. + * + * This must be a three-dimensional, Cartesian grid which has an active + * cluster which is without holes and which is convex. The UnstructuredGrid + * structure is used because it is the lingua franca of grids in the + * simulator, not because this method will handle every possible grid. + * + * This pointer is NOT adopted. The caller must still dispose of the grid. + * + * @return Upscaled, fine grid. + * + * The caller have the responsibility of disposing this grid; no other + * references will initially exist. + */ + static TopSurf* create (const UnstructuredGrid& fine); + +private: + /** + * @brief You are not meant to construct these yourself; use create (). + */ + TopSurf (); +}; + +} // namespace Opm + +#endif // OPM_VERTEQ_TOPSURF_HPP_INCLUDED diff --git a/opm/grid/verteq/utility/exc.cpp b/opm/grid/verteq/utility/exc.cpp new file mode 100644 index 0000000000..5a8fe36b53 --- /dev/null +++ b/opm/grid/verteq/utility/exc.cpp @@ -0,0 +1,95 @@ +/** + * Copyright (C) 2013 Uni Research AS. + * This code is licensed under The GNU General Public License v3.0 or later. + */ +#include +#include + +#include +#include +#include +#include +#include // va_list, va_start, va_end +#include // vsnprintf +#include // std::string +#ifdef __MSVC__ +# include // alloca +#endif + +namespace Opm { +namespace Exc { + +// this symbol must be exported so that other DSO can get it +struct SYMBOL_IS_EXPORTED message; + +typedef boost::error_info message_t; + +Base& +Base::operator () (char const* fmt, ...) { + // number of bytes necessary to render this string + va_list params; + va_start (params, fmt); + int count = vsnprintf (0, 0, fmt, params) + 1; + va_end (params); + + // allocate these bytes; this is necessary since we aren't allowed + // to touch the internal buffer of a std::string. we make up for it + // by allocating on the stack instead of on the free store! + size_t bytes = count * sizeof (char); + char* buf = static_cast (alloca (bytes)); + + // print the string to those bytes. note that it is important to + // restart the argument list, it cannot be reused from above! + va_start (params, fmt); + vsnprintf (buf, bytes, fmt, params); + va_end (params); + + // put a copy of the buffer inside the exception structure + // see + *this << message_t (std::string (buf)); + + // allow the exception object to be chained + return *this; +} + +// this message is used if nothing else is specified. it is in the +// anonymous namespace so that no-one else is able to use it. +namespace { +char const* UNSPECIFIED = ""; +} + +char const* +Base::what () const throw () { + // retrieve the stored reason, or a generic message if there is none + std::string const* str = boost::get_error_info (*this); + return str ? str->c_str () : UNSPECIFIED; +} + +std::string +diag_what (std::exception const& ex) { + // header + std::stringstream buf; + buf << "Error"; + // std::exception info; this contains the reason + char const* what = ex.what(); + if (what != UNSPECIFIED) { + buf << ": " << what; + } + // boost::exception info; this contains the location + boost::exception const* bex = dynamic_cast (&ex); + if (bex) { + char const* const* file = boost::get_error_info (*bex); + int const* line = boost::get_error_info (*bex); + char const* const* func = boost::get_error_info (*bex); + if (func) { + buf << ", in function " << *func; + } + if (file && line) { + buf << ", at " << *file << ":" << *line; + } + } + return buf.str(); +} + +} /* namespace Exc */ +} /* namespace Opm */ diff --git a/opm/grid/verteq/utility/exc.hpp b/opm/grid/verteq/utility/exc.hpp new file mode 100644 index 0000000000..bd01e6e756 --- /dev/null +++ b/opm/grid/verteq/utility/exc.hpp @@ -0,0 +1,106 @@ +/** + * Copyright (C) 2013 Uni Research AS. + * This code is licensed under The GNU General Public License v3.0 or later. + */ +#include // std::exception +#ifndef UUID_274DA366004E11DCB1DDFE2E56D89593 +#include // boost::exception +#endif +#ifndef UUID_8D22C4CA9CC811DCAA9133D256D89593 +#include // boost::throw_{function,file,line} +#endif + +#ifndef OPM_VERTEQ_VISIBILITY_INCLUDED +#include +#endif +#if defined (opmverteq_EXPORTS) +# define OPM_VERTEQ_PUBLIC SYMBOL_IS_EXPORTED +#else +# define OPM_VERTEQ_PUBLIC SYMBOL_IS_IMPORTED +#endif +#define OPM_VERTEQ_PRIVATE SYMBOL_IS_LOCALDEF + +// enable printf-style attribute for MSVC +#ifdef _MSC_VER +# if _MSC_VER >= 1400 +# define _USE_ATTRIBUTES_FOR_SAL 1 +# include +# endif +#endif + +namespace Opm { +namespace Exc { +/** + * Base class for exceptions thrown from the OPM framework. + * + * Throw exceptions like this: + * + * throw OPM_EXC ("Magic value = %d", 42); + * + * Catch exceptions like this: + * + * catch (std::exception& ex) { + * std::cerr << OPM_WHAT (ex) << std::endl; + * exit (-1); + * } + */ +struct OPM_VERTEQ_PUBLIC Base + : public virtual std::exception + , public virtual boost::exception { + + /** + * @brief Add a message to the string + * @param fmt printf-style format string. The rest of the parameters + * are formatted according to this string. + * @return Same exception object so it can be chained. (In particular, + * more boost::error_info objects may be added). + */ + virtual Base& operator () ( +#ifdef _MSC_VER +# if _MSC_VER >= 1400 + __format_string +# endif +#endif + char const* fmt, ...) + // there is an implicit this parameter at the start, so the format + // string is the second parameter, and the variable argument list + // starts with the third. +#ifdef __GNUC__ + __attribute__ ((format (printf, 2, 3))) +#endif + ; + + /** + * @brief Message created at the throw-site about the error. + * @return String that can be printed in the log about the exception. + */ + virtual char const* what () const throw (); +}; + +/** + * @brief Retrieve information about the code that failed + * @param ex Exception that was thrown + * @return Text containing error information and location + */ +std::string OPM_VERTEQ_PUBLIC diag_what (std::exception const& ex); + +/** +* Create a new exception object, possibly filled with location +* information if a debug build was done. +*/ +#ifdef DEBUG +// const_cast is necessary because the operator<< returns a +// const, and we may want it mutable to add more information + #define OPM_EXC const_cast (Opm::Exc::Base ()\ + << ::boost::throw_function (BOOST_CURRENT_FUNCTION)\ + << ::boost::throw_file (__FILE__)\ + << ::boost::throw_line (static_cast (__LINE__))\ +) +#else + #define OPM_EXC Opm::Exc::Base () +#endif + +#define OPM_WHAT(ex) Opm::Exc::diag_what (ex) + +} /* namespace Opm::Exc */ +} /* namespace Opm */ diff --git a/opm/grid/verteq/utility/runlen.cpp b/opm/grid/verteq/utility/runlen.cpp new file mode 100644 index 0000000000..da0d13497d --- /dev/null +++ b/opm/grid/verteq/utility/runlen.cpp @@ -0,0 +1,17 @@ +#include "runlen.hpp" +#include +using namespace Opm; + +rlw_int Opm::grid_cell_facetag (const UnstructuredGrid& g) { + // tag for faces in cells + return rlw_int (g.number_of_cells, + g.cell_facepos, + g.cell_facetag); +} + +rlw_int Opm::grid_cell_faces (const UnstructuredGrid& g) { + // id for faces in cells + return rlw_int (g.number_of_cells, + g.cell_facepos, + g.cell_faces); +} diff --git a/opm/grid/verteq/utility/runlen.hpp b/opm/grid/verteq/utility/runlen.hpp new file mode 100644 index 0000000000..1e9c30c42b --- /dev/null +++ b/opm/grid/verteq/utility/runlen.hpp @@ -0,0 +1,199 @@ +#ifndef OPM_VERTEQ_RUNLEN_HPP_INCLUDED +#define OPM_VERTEQ_RUNLEN_HPP_INCLUDED + +// Copyright (C) 2013 Uni Research AS +// This file is licensed under the GNU General Public License v3.0 + +// forward declaration +struct UnstructuredGrid; + +namespace Opm { + +/** + * Regards a set of (member) variables as a run-length encoded matrix. + * + * Each column can have a variable number of rows. Although the *values* + * of the matrix can be changed, its sparsity cannot, i.e. one cannot + * remove or add new elements to a column. + * + * Use this class to access and iterate over a run-length encoded matrix + * in the format that is used by UnstructuredGrid without having to worry + * about getting the indexing right. + * + * @tparam T Datatype for the extra data that should be stored for + * each element, e.g. double. + * + * @example + * @code{.cpp} + * RunLenView faces_in_cell ( + * g.number_of_cells, + * g.cell_facepos, + * g.cell_faces + * ); + * + * int num_local_faces = faces_in_cell.size (cellno); + * int first_local_face = faces_in_cell [cellno] [0]; + * @endcode + * + * Notice if you want to loop through every item and know where you are + * (because you intend to use this as an index in another matrix), you + * can do: + * + * @example + * @code{.cpp} + * RunLenView cell_faces ( + * g.number_of_cells, + * g.cell_facepos, + * g.cell_faces + * ); + * + * for (int cell = 0; cell < cell_faces.cols(); ++cell) { + * for (int local_face = 0; local_face < cell_faces.size (cell); ++local_face) { + * int face_id = cell_faces[cell][local_face]; + * } + * } + * @endcode + * + * @see Opm::RunLenData + */ +template +class RunLenView { +protected: + /** + * Size information. pos has num_of_cols+1 items, pos[i] contains + * the starting index of the data values for column i. (Since there + * is one more element than there are columns, the last one is the + * total number of elements). The number 0 is explicitly stored in + * the first column to avoid special processing. + */ + int num_of_cols; + int* pos; + + /** + * Data for each of the individual elements, stored consecutively + * for each column located together, followed by the next column. + */ + T* data; + +public: + /** + * Construct a view into the run-length encoded matrix. The view is + * only defined as long as the underlaying structures exist! + * + * @param number Number of cells (elements, faces, nodes) + * @param pos_ptr Table of starting indices. + * If columns are "foo" and rows are "bar", then this + * is the member called "foo_barpos". + * @param values Actual data storage. + * If columns are "foo" and rows are "bar", then this + * is the member called "foo_bars". + */ + RunLenView (int num_cols, int* pos_ptr, T* values) + // store them locally for later use + : num_of_cols (num_cols) + , pos (pos_ptr) + , data (values) { + } + + /** + * Create another view of the same data. + * + * @param rhs View to a run-length-encoded matrix + */ + RunLenView (const RunLenView& rhs) + // copy all fields verbatim + : num_of_cols (rhs.num_of_cols) + , pos (rhs.pos) + , data (rhs.data) { + } + + /** + * Access a column directly. + * + * @param col Index of the column to get + * @return Pointer to the start of the column + */ + T* operator [] (int col) const { + return &data [pos [col]]; + } + + /** + * Number of columns that are stored in the entire matrix. + * + * @return Number of columns. + */ + int cols () const { + return num_of_cols; + } + + /** + * Number of elements that are stored in one particular column. + * + * @param col Index of the column. + * @return Number of elements. + */ + int size (int col) const { + return pos [col + 1] - pos [col]; + } + + /** + * Quick accessor to get the last element. When we store accumulated + * data in the array, this will quickly give us the total. + * + * Note that this is NOT the end iterator for the column. + * + * @param col Index of the column + * @return Value of the last element. If there is no elements in + * this column, then the return value is undefined. + */ + T& last (int col) const { + return data [pos [col + 1] - 1]; + } +}; + +/** + * Allocate a new vector of data for each element, accessible as + * a zig-zag matrix. + * + * Use this kind of matrix when you want to enhance the grid structure + * with some information per element, using the existing format. + * + * @see Opm::RunLenView + */ +template +struct RunLenData : public RunLenView { + /** + * Allocate a matrix based on sizes specified elsewhere. This is + * useful if you want to supply with your own data. + * + * @param number Number of entities (rows). + * @param pos_ptr Number of elements to allocate in front of each + * entity; starting index in the array for this row. + * + * @see Opm::RunLenView::RunLenView + */ + RunLenData (int number, int* pos_ptr) + // allocate a new vector for the data, containing the needed + // number of elements. note that there is only one new + // operation is the parameter list, so there is no leakage if + // an out-of-memory exception is thrown. + : RunLenView (number, pos_ptr, new T [pos_ptr [number]]) { + } + + ~RunLenData () { + // this member is initialized with data allocated in our ctor + delete [] RunLenView ::data; + } +}; + +// shorthands for most used types +typedef const RunLenView rlw_int; +typedef const RunLenView rlw_double; + +// access common run-length encoded matrices in a grid structure +rlw_int grid_cell_facetag (const UnstructuredGrid& g); +rlw_int grid_cell_faces (const UnstructuredGrid& g); + +} /* namespace Opm */ + +#endif /* OPM_VERTEQ_RUNLEN_HPP_INCLUDED */ diff --git a/opm/grid/verteq/utility/visibility.h b/opm/grid/verteq/utility/visibility.h new file mode 100644 index 0000000000..149a5f23c7 --- /dev/null +++ b/opm/grid/verteq/utility/visibility.h @@ -0,0 +1,42 @@ +#ifndef OPM_VERTEQ_VISIBILITY_INCLUDED +#define OPM_VERTEQ_VISIBILITY_INCLUDED + +/** + * Macros to encapsulate symbol visibility on various platforms. + * You need to define a separate macro for your package; symbols + * that are exported in one dynamic shared object, may of course + * be imported by another one! + * + * Usage: + * + * #include + * + * #if defined (foo_EXPORTS) + * # define FOO_PUBLIC SYMBOL_IS_EXPORTED + * #else + * # define FOO_PUBLIC SYMBOL_IS_IMPORTED + * #endif + * #define FOO_PRIVATE SYMBOL_IS_LOCALDEF + * + * struct FOO_PUBLIC Bar { + * }; + * + * int FOO_PRIVATE mumble (); + */ +#if defined (_WIN32) +# define SYMBOL_IS_EXPORTED __declspec (dllexport) +# define SYMBOL_IS_IMPORTED __declspec (dllimport) +# define SYMBOL_IS_LOCALDEF +#else +# if __GNUC__ >= 4 +# define SYMBOL_IS_EXPORTED __attribute__ ((visibility ("default"))) +# define SYMBOL_IS_IMPORTED __attribute__ ((visibility ("default"))) +# define SYMBOL_IS_LOCALDEF __attribute__ ((visibility ("hidden"))) +# else +# define SYMBOL_IS_EXPORTED +# define SYMBOL_IS_IMPORTED +# define SYMBOL_IS_LOCALDEF +# endif +#endif + +#endif /* OPM_VERTEQ_VISIBILITY_INCLUDED */ diff --git a/tests/test_cpgrid.cpp b/tests/test_cpgrid.cpp index e11178c08f..277c383e54 100644 --- a/tests/test_cpgrid.cpp +++ b/tests/test_cpgrid.cpp @@ -5,12 +5,17 @@ #include #include +#include #include #include #include +#include +#include + +#include #define DISABLE_DEPRECATED_METHOD_CHECK 1 using Dune::referenceElement; //grid check assume usage of Dune::Geometry @@ -24,6 +29,7 @@ using Dune::referenceElement; //grid check assume usage of Dune::Geometry #include + template void testGridIteration( const GridView& gridView, const int nElem ) { @@ -85,20 +91,20 @@ void testGridIteration( const GridView& gridView, const int nElem ) } + template void testGrid(Grid& grid, const std::string& name, const size_t nElem, const size_t nVertices) { - typedef typename Grid::LeafGridView GridView; - /* + //testVerteq( grid ); - try { - gridcheck( grid ); - } - catch ( const Dune::Exception& e) - { - std::cerr << "Warning: " << e.what() << std::endl; - } -*/ + typedef typename Grid::LeafGridView GridView; + // try { + // gridcheck( grid ); + // } + // catch ( const Dune::Exception& e) + // { + // std::cerr << "Warning: " << e.what() << std::endl; + // } std::cout << name << std::endl; testGridIteration( grid.leafGridView(), nElem ); @@ -132,16 +138,94 @@ void testGrid(Grid& grid, const std::string& name, const size_t nElem, const siz } + int main(int argc, char** argv ) { // initialize MPI Dune::MPIHelper::instance( argc, argv ); + std::stringstream dgfFile; + // create unit cube with 4 cells in each direction + dgfFile << "DGF" << std::endl; + dgfFile << "Interval" << std::endl; + dgfFile << "0 0 0" << std::endl; + dgfFile << "4 4 4" << std::endl; + dgfFile << "4 4 4" << std::endl; + dgfFile << "#" << std::endl; - // test CpGrid - typedef Dune::CpGrid Grid; + // test PolyhedralGrid + { + typedef Dune::PolyhedralGrid< 3, 3 > Grid; +#if HAVE_ECL_INPUT + const char *deckString = + "RUNSPEC\n" + "METRIC\n" + "DIMENS\n" + "4 4 4 /\n" + "GRID\n" + "DXV\n" + "4*1 /\n" + "DYV\n" + "4*1 /\n" + "DZ\n" + "16*1 /\n" + "TOPS\n" + "16*100.0 /\n"; + + Opm::Parser parser; + const auto deck = parser.parseString(deckString); + std::vector porv; + + std::cout << "Check ecl grid" << std::endl; + Grid grid(deck, porv); + gridcheck( grid ); + //testGrid( grid, "polyhedralgrid", 8, 27 ); + /* + Opm::TopSurf* ts; + ts = Opm::TopSurf::create (grid); + std::cout << ts->dimensions << std::endl; + std::cout << ts->number_of_cells <<" " << ts->number_of_faces << " " << ts->number_of_nodes << " " << std::endl; + //for (int i = 0; i < ts->number_of_nodes*ts->dimensions; ++i) + // std::cout << ts->node_coordinates[i] << std::endl; + typedef Dune::PolyhedralGrid< 2, 2 > Grid2D; + + std::cout << "tsDune for " << std::endl; + Grid2D tsDune (*ts); + std::cout << "tsDune after " << std::endl; + testGrid ( tsDune, "ts", 27 ); + */ +#endif + { + std::cout <<"Check 3d grid" << std::endl; + Dune::GridPtr< Grid > gridPtr( dgfFile ); + //testGrid( *gridPtr, "polyhedralgrid-dgf", std::pow(3, int(Grid::dimension)) ); + gridcheck( *gridPtr ); + } + + { + std::cout <<"Check 2d grid" << std::endl; + std::stringstream dgfFile; + // create unit cube with 8 cells in each direction + dgfFile << "DGF" << std::endl; + dgfFile << "Interval" << std::endl; + dgfFile << "0 0" << std::endl; + dgfFile << "1 1" << std::endl; + dgfFile << "2 2" << std::endl; + dgfFile << "#" << std::endl; + typedef Dune::PolyhedralGrid< 2, 2 > Grid; + Dune::GridPtr< Grid > gridPtr( dgfFile ); + + std::cout << "Grididm = " << int(Grid::dimension) << std::endl; + size_t nVx = 9; //std::pow(int(3), int(Grid::dimension)); + //testGrid( *gridPtr, "polyhedralgrid-dgf", 4, nVx); + gridcheck( *gridPtr ); + } + } + // test CpGrid + { + typedef Dune::CpGrid Grid; #if HAVE_ECL_INPUT - const char *deckString = + const char *deckString = "RUNSPEC\n" "METRIC\n" "DIMENS\n" @@ -156,29 +240,20 @@ int main(int argc, char** argv ) "TOPS\n" "8*100.0 /\n"; - Opm::Parser parser; - const auto deck = parser.parseString(deckString); - std::vector porv; + Opm::Parser parser; + const auto deck = parser.parseString(deckString); + std::vector porv; - Grid grid; - const int* actnum = deck.hasKeyword("ACTNUM") ? deck.getKeyword("ACTNUM").getIntData().data() : nullptr; - Opm::EclipseGrid ecl_grid(deck , actnum); + Grid grid; + const int* actnum = deck.hasKeyword("ACTNUM") ? deck.getKeyword("ACTNUM").getIntData().data() : nullptr; + Opm::EclipseGrid ecl_grid(deck , actnum); - grid.processEclipseFormat(&ecl_grid, false, false, false, porv); - testGrid( grid, "CpGrid_ecl", 8, 27 ); + grid.processEclipseFormat(&ecl_grid, false, false, false, porv); + testGrid( grid, "CpGrid_ecl", 8, 27 ); #endif - std::stringstream dgfFile; - // create unit cube with 8 cells in each direction - dgfFile << "DGF" << std::endl; - dgfFile << "Interval" << std::endl; - dgfFile << "0 0 0" << std::endl; - dgfFile << "4 4 4" << std::endl; - dgfFile << "4 4 4" << std::endl; - dgfFile << "#" << std::endl; - - Dune::GridPtr< Grid > gridPtr( dgfFile ); - testGrid( *gridPtr, "CpGrid_dgf", 64, 125 ); - + Dune::GridPtr< Grid > gridPtr( dgfFile ); + testGrid( *gridPtr, "CpGrid_dgf", 64, 125 ); + } return 0; }