From a64d066ad3788b668df4580266ad056e791c7416 Mon Sep 17 00:00:00 2001 From: norce Date: Fri, 8 Mar 2019 10:37:10 +0100 Subject: [PATCH 01/22] Incuded opm-verteq 2d grid into opm-grid. --- CMakeLists_files.cmake | 11 +- opm/grid/cpgrid/dgfparser.hh | 2 +- opm/grid/polyhedralgrid/grid.hh | 18 +- opm/grid/verteq/nav.cpp | 92 +++ opm/grid/verteq/nav.hpp | 425 ++++++++++++++ opm/grid/verteq/topsurf.cpp | 826 +++++++++++++++++++++++++++ opm/grid/verteq/topsurf.hpp | 177 ++++++ opm/grid/verteq/utility/exc.cpp | 95 +++ opm/grid/verteq/utility/exc.hpp | 106 ++++ opm/grid/verteq/utility/runlen.cpp | 17 + opm/grid/verteq/utility/runlen.hpp | 199 +++++++ opm/grid/verteq/utility/visibility.h | 42 ++ tests/grid_test.cc | 31 +- 13 files changed, 2031 insertions(+), 10 deletions(-) create mode 100644 opm/grid/verteq/nav.cpp create mode 100644 opm/grid/verteq/nav.hpp create mode 100644 opm/grid/verteq/topsurf.cpp create mode 100644 opm/grid/verteq/topsurf.hpp create mode 100644 opm/grid/verteq/utility/exc.cpp create mode 100644 opm/grid/verteq/utility/exc.hpp create mode 100644 opm/grid/verteq/utility/runlen.cpp create mode 100644 opm/grid/verteq/utility/runlen.hpp create mode 100644 opm/grid/verteq/utility/visibility.h diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 24ee6395b..ef448ca69 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -57,6 +57,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) @@ -96,7 +100,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_geom2d.cpp tests/test_gridutilities.cpp tests/test_minpvprocessor.cpp -# tests/grid_test.cc + tests/grid_test.cc tests/p2pcommunicator_test.cc tests/test_repairzcorn.cpp tests/test_sparsetable.cpp @@ -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/cpgrid/dgfparser.hh b/opm/grid/cpgrid/dgfparser.hh index c8d5f12b9..649cee84d 100644 --- a/opm/grid/cpgrid/dgfparser.hh +++ b/opm/grid/cpgrid/dgfparser.hh @@ -21,7 +21,7 @@ #define DUNE_CPGRID_DGFPARSER_HH #include -#include +#include namespace Dune { diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index f23d63e08..6b067765d 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1113,9 +1113,12 @@ namespace Dune const int numCells = size( 0 ); cellVertices_.resize( numCells ); + std::cout << numCells << std::endl; + // sort vertices such that they comply with the dune cube reference element if( grid_.cell_facetag ) { + std::cout << "facetag" << std::endl; typedef std::array KeyType; std::map< const KeyType, const int > vertexFaceTags; const int vertexFacePattern [8][3] = { @@ -1140,6 +1143,8 @@ namespace Dune vertexFaceTags.insert( std::make_pair( key, i ) ); } + std::cout << "cells before" << std::endl; + for (int c = 0; c < numCells; ++c) { typedef std::map vertexmap_t; @@ -1167,15 +1172,20 @@ namespace Dune } } } + std::cout << "face node done" << std::endl; + std::cout << cell_pts.size() << std::endl; typedef std::map< int, std::set > vertexlist_t; vertexlist_t vertexList; - for( int faceTag = 0; faceTag<6; ++faceTag ) + for( int faceTag = 0; faceTag cell_pts; diff --git a/opm/grid/verteq/nav.cpp b/opm/grid/verteq/nav.cpp new file mode 100644 index 000000000..8efc315cc --- /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 000000000..61583127d --- /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 000000000..3a371e152 --- /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 000000000..0be3b3037 --- /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 000000000..5a8fe36b5 --- /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 000000000..bd01e6e75 --- /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 000000000..da0d13497 --- /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 000000000..1e9c30c42 --- /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 000000000..149a5f23c --- /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/grid_test.cc b/tests/grid_test.cc index 78ae8f2df..b7fc086d9 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -4,13 +4,15 @@ #include #include -#include -#include +#include +#include #include #include #include -#include +#include + +#include #define DISABLE_DEPRECATED_METHOD_CHECK 1 #if DUNE_VERSION_NEWER(DUNE_GRID,2,5) @@ -116,6 +118,7 @@ void testGrid(Grid& grid, const std::string& name) std::cerr << "Warning: " << e.what() << std::endl; } #endif + std::cout << name << std::endl; testGridIteration( grid.leafGridView() ); @@ -126,11 +129,11 @@ void testGrid(Grid& grid, const std::string& name) std::cout << "VertexMapper.size(): " << mapper.size() << "\n"; if (mapper.size() != 27) { std::cout << "Wrong size of vertex mapper. Expected 27!\n"; - std::abort(); + //std::abort(); } // VTKWriter does not work with geometry type none at the moment - if( grid.geomTypes( 0 )[ 0 ].isCube() ) + if( true || grid.geomTypes( 0 )[ 0 ].isCube() ) { std::cout << "create vtkWriter\n"; typedef Dune::VTKWriter VtkWriter; @@ -165,8 +168,7 @@ int main(int argc, char** argv ) #if HAVE_ECL_INPUT Opm::Parser parser; - Opm::ParseContext parseContext; - const auto deck = parser.parseString(deckString , parseContext); + const auto deck = parser.parseString(deckString); std::vector porv; #endif @@ -176,7 +178,21 @@ int main(int argc, char** argv ) #if HAVE_ECL_INPUT Grid grid(deck, porv); testGrid( grid, "polyhedralgrid" ); + 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 før " << std::endl; + Grid2D tsDune (*ts); + std::cout << "tsDune etter " << std::endl; + testGrid ( tsDune, "ts"); + #endif + + //Dune::GridPtr< Grid > gridPtr( dgfFile ); //testGrid( *gridPtr, "polyhedralgrid-dgf" ); } @@ -194,6 +210,7 @@ int main(int argc, char** argv ) Opm::UgGridHelpers::createEclipseGrid( grid , ecl_grid ); testGrid( grid, "cpgrid2" ); + #endif Dune::GridPtr< Grid > gridPtr( dgfFile ); //testGrid( *gridPtr, "cpgrid-dgf" ); From bc7bd9d2b21f0ed8616ca39defb67e804b9b6535 Mon Sep 17 00:00:00 2001 From: norce Date: Fri, 8 Mar 2019 10:37:10 +0100 Subject: [PATCH 02/22] Incuded opm-verteq 2d grid into opm-grid. --- CMakeLists_files.cmake | 11 +- opm/grid/cpgrid/dgfparser.hh | 2 +- opm/grid/polyhedralgrid/grid.hh | 18 +- opm/grid/verteq/nav.cpp | 92 +++ opm/grid/verteq/nav.hpp | 425 ++++++++++++++ opm/grid/verteq/topsurf.cpp | 826 +++++++++++++++++++++++++++ opm/grid/verteq/topsurf.hpp | 177 ++++++ opm/grid/verteq/utility/exc.cpp | 95 +++ opm/grid/verteq/utility/exc.hpp | 106 ++++ opm/grid/verteq/utility/runlen.cpp | 17 + opm/grid/verteq/utility/runlen.hpp | 199 +++++++ opm/grid/verteq/utility/visibility.h | 42 ++ tests/grid_test.cc | 31 +- 13 files changed, 2031 insertions(+), 10 deletions(-) create mode 100644 opm/grid/verteq/nav.cpp create mode 100644 opm/grid/verteq/nav.hpp create mode 100644 opm/grid/verteq/topsurf.cpp create mode 100644 opm/grid/verteq/topsurf.hpp create mode 100644 opm/grid/verteq/utility/exc.cpp create mode 100644 opm/grid/verteq/utility/exc.hpp create mode 100644 opm/grid/verteq/utility/runlen.cpp create mode 100644 opm/grid/verteq/utility/runlen.hpp create mode 100644 opm/grid/verteq/utility/visibility.h diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 24ee6395b..ef448ca69 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -57,6 +57,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) @@ -96,7 +100,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_geom2d.cpp tests/test_gridutilities.cpp tests/test_minpvprocessor.cpp -# tests/grid_test.cc + tests/grid_test.cc tests/p2pcommunicator_test.cc tests/test_repairzcorn.cpp tests/test_sparsetable.cpp @@ -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/cpgrid/dgfparser.hh b/opm/grid/cpgrid/dgfparser.hh index c8d5f12b9..649cee84d 100644 --- a/opm/grid/cpgrid/dgfparser.hh +++ b/opm/grid/cpgrid/dgfparser.hh @@ -21,7 +21,7 @@ #define DUNE_CPGRID_DGFPARSER_HH #include -#include +#include namespace Dune { diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index f23d63e08..6b067765d 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1113,9 +1113,12 @@ namespace Dune const int numCells = size( 0 ); cellVertices_.resize( numCells ); + std::cout << numCells << std::endl; + // sort vertices such that they comply with the dune cube reference element if( grid_.cell_facetag ) { + std::cout << "facetag" << std::endl; typedef std::array KeyType; std::map< const KeyType, const int > vertexFaceTags; const int vertexFacePattern [8][3] = { @@ -1140,6 +1143,8 @@ namespace Dune vertexFaceTags.insert( std::make_pair( key, i ) ); } + std::cout << "cells before" << std::endl; + for (int c = 0; c < numCells; ++c) { typedef std::map vertexmap_t; @@ -1167,15 +1172,20 @@ namespace Dune } } } + std::cout << "face node done" << std::endl; + std::cout << cell_pts.size() << std::endl; typedef std::map< int, std::set > vertexlist_t; vertexlist_t vertexList; - for( int faceTag = 0; faceTag<6; ++faceTag ) + for( int faceTag = 0; faceTag cell_pts; diff --git a/opm/grid/verteq/nav.cpp b/opm/grid/verteq/nav.cpp new file mode 100644 index 000000000..8efc315cc --- /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 000000000..61583127d --- /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 000000000..3a371e152 --- /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 000000000..0be3b3037 --- /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 000000000..5a8fe36b5 --- /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 000000000..bd01e6e75 --- /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 000000000..da0d13497 --- /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 000000000..1e9c30c42 --- /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 000000000..149a5f23c --- /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/grid_test.cc b/tests/grid_test.cc index 78ae8f2df..b7fc086d9 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -4,13 +4,15 @@ #include #include -#include -#include +#include +#include #include #include #include -#include +#include + +#include #define DISABLE_DEPRECATED_METHOD_CHECK 1 #if DUNE_VERSION_NEWER(DUNE_GRID,2,5) @@ -116,6 +118,7 @@ void testGrid(Grid& grid, const std::string& name) std::cerr << "Warning: " << e.what() << std::endl; } #endif + std::cout << name << std::endl; testGridIteration( grid.leafGridView() ); @@ -126,11 +129,11 @@ void testGrid(Grid& grid, const std::string& name) std::cout << "VertexMapper.size(): " << mapper.size() << "\n"; if (mapper.size() != 27) { std::cout << "Wrong size of vertex mapper. Expected 27!\n"; - std::abort(); + //std::abort(); } // VTKWriter does not work with geometry type none at the moment - if( grid.geomTypes( 0 )[ 0 ].isCube() ) + if( true || grid.geomTypes( 0 )[ 0 ].isCube() ) { std::cout << "create vtkWriter\n"; typedef Dune::VTKWriter VtkWriter; @@ -165,8 +168,7 @@ int main(int argc, char** argv ) #if HAVE_ECL_INPUT Opm::Parser parser; - Opm::ParseContext parseContext; - const auto deck = parser.parseString(deckString , parseContext); + const auto deck = parser.parseString(deckString); std::vector porv; #endif @@ -176,7 +178,21 @@ int main(int argc, char** argv ) #if HAVE_ECL_INPUT Grid grid(deck, porv); testGrid( grid, "polyhedralgrid" ); + 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 før " << std::endl; + Grid2D tsDune (*ts); + std::cout << "tsDune etter " << std::endl; + testGrid ( tsDune, "ts"); + #endif + + //Dune::GridPtr< Grid > gridPtr( dgfFile ); //testGrid( *gridPtr, "polyhedralgrid-dgf" ); } @@ -194,6 +210,7 @@ int main(int argc, char** argv ) Opm::UgGridHelpers::createEclipseGrid( grid , ecl_grid ); testGrid( grid, "cpgrid2" ); + #endif Dune::GridPtr< Grid > gridPtr( dgfFile ); //testGrid( *gridPtr, "cpgrid-dgf" ); From 7bed2aff9bd130f1eee58d940783719897923f79 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 15 Mar 2019 14:02:32 +0100 Subject: [PATCH 03/22] [feature][PolyhedralGrid] merge from gridfactory branch. This also adds a proper DGF implementation for Polyhedron and Polygon. --- opm/grid/polyhedralgrid/dgfparser.hh | 313 +++++++++++++++++-- opm/grid/polyhedralgrid/entity.hh | 2 +- opm/grid/polyhedralgrid/grid.hh | 364 +++++++++++++++++----- opm/grid/polyhedralgrid/gridfactory.hh | 264 ++++++++++++++++ opm/grid/polyhedralgrid/idset.hh | 22 +- opm/grid/polyhedralgrid/indexset.hh | 14 + opm/grid/polyhedralgrid/polyhedralmesh.hh | 244 +++++++++++++++ tests/grid_test.cc | 14 +- 8 files changed, 1112 insertions(+), 125 deletions(-) create mode 100644 opm/grid/polyhedralgrid/gridfactory.hh create mode 100644 opm/grid/polyhedralgrid/polyhedralmesh.hh diff --git a/opm/grid/polyhedralgrid/dgfparser.hh b/opm/grid/polyhedralgrid/dgfparser.hh index c426cad89..9dd9ff764 100644 --- a/opm/grid/polyhedralgrid/dgfparser.hh +++ b/opm/grid/polyhedralgrid/dgfparser.hh @@ -3,19 +3,134 @@ #ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH #define DUNE_POLYHEDRALGRID_DGFPARSER_HH +#include +#include + #include +#include +#include #include -#include +#if DUNE_VERSION_NEWER(DUNE_GRID,2,5) +#include +#endif + +#include + +#if HAVE_OPM_PARSER +#include +#include +#include +#endif namespace Dune { -#warning TODO: non-trivial DGFGridFactory + namespace dgf + { + +#if ! DUNE_VERSION_NEWER(DUNE_GRID,2,5) + namespace PolyhedralGrid + { + + // PolygonBlock + // ------------ + + struct PolygonBlock + : public BasicBlock + { + PolygonBlock ( std::istream &in, int numVtx, int vtxOfs ) + : BasicBlock( in, "Polygon" ), vtxBegin_( vtxOfs ), vtxEnd_( vtxOfs + numVtx ) + {} + + int get ( std::vector< std::vector< int > > &polygons ) + { + reset(); + std::vector< int > polygon; + while( getnextline() ) + { + polygon.clear(); + for( int vtxIdx; getnextentry( vtxIdx ); ) + { + if( (vtxBegin_ > vtxIdx) || (vtxIdx >= vtxEnd_) ) + DUNE_THROW( DGFException, "Error in " << *this << ": Invalid vertex index (" << vtxIdx << " not int [" << vtxBegin_ << ", " << vtxEnd_ << "[)" ); + polygon.push_back( vtxIdx - vtxBegin_ ); + } + + polygons.push_back( polygon ); + } + return polygons.size(); + } + + private: + int vtxBegin_, vtxEnd_; + }; + + + + // PolyhedronBlock + // --------------- + + struct PolyhedronBlock + : public BasicBlock + { + explicit PolyhedronBlock ( std::istream &in, int numPolys ) + : BasicBlock( in, "Polyhedron" ), numPolys_( numPolys ) + {} + + int get ( std::vector< std::vector< int > > &polyhedra ) + { + reset(); + std::vector< int > polyhedron; + int minPolyId = 1; + while( getnextline() ) + { + polyhedron.clear(); + for( int polyIdx; getnextentry( polyIdx ); ) + { + if( (polyIdx < 0) || (polyIdx > numPolys_) ) + DUNE_THROW( DGFException, "Error in " << *this << ": Invalid polygon index (" << polyIdx << " not int [0, " << numPolys_ << "])" ); + + minPolyId = std::min( minPolyId, polyIdx ); + polyhedron.push_back( polyIdx ); + } + + polyhedra.push_back( polyhedron ); + } + + // substract minimal number to have 0 starting numbering + if( minPolyId > 0 ) + { + const size_t polySize = polyhedra.size(); + for( size_t i=0; i struct DGFGridFactory< PolyhedralGrid< dim, dimworld > > @@ -27,80 +142,224 @@ namespace Dune typedef typename Grid::template Codim<0>::Entity Element; typedef typename Grid::template Codim::Entity Vertex; - explicit DGFGridFactory ( std::istream &input, - MPICommunicator comm = MPIHelper::getCommunicator() ) - : grid_( nullptr ) + explicit DGFGridFactory ( std::istream &input, MPICommunicator comm = MPIHelper::getCommunicator() ) + : grid_() { + input.clear(); + input.seekg( 0 ); + if( !input ) + DUNE_THROW( DGFException, "Error resetting input stream" ); + generate( input ); } - explicit DGFGridFactory ( const std::string &filename, - MPICommunicator comm = MPIHelper::getCommunicator() ) - : grid_( nullptr ) + explicit DGFGridFactory ( const std::string &filename, MPICommunicator comm = MPIHelper::getCommunicator() ) + : grid_() { - } + std::ifstream input( filename ); + if( !input ) + DUNE_THROW( DGFException, "Macrofile '" << filename << "' not found" ); - Grid *grid () const - { - return grid_; +#if HAVE_OPM_PARSER + if( !DuneGridFormatParser::isDuneGridFormat( input ) ) + { + Opm::Parser parser; + Opm::ParseContext parseContext; + const auto deck = parser.parseFile(filename, parseContext); + std::vector porv; + + grid_.reset( new Grid( deck, porv ) ); + return ; + } + else +#endif + { + generate( input ); + } } + Grid *grid () const { return grid_.operator->(); } + template< class Intersection > bool wasInserted ( const Intersection &intersection ) const { - return false; + return false; } template< class Intersection > int boundaryId ( const Intersection &intersection ) const { - return false; + return false; } - bool haveBoundaryParameters () const - { - return false; - } + bool haveBoundaryParameters () const { return false; } template< int codim > int numParameters () const { - return 0; + //return (codim == dimension ? numVtxParams_ : 0);; + return 0; } template< class Intersection > const typename DGFBoundaryParameter::type & boundaryParameter ( const Intersection &intersection ) const { - return DGFBoundaryParameter::defaultValue(); + return DGFBoundaryParameter::defaultValue();; } template< class Entity > std::vector< double > ¶meter ( const Entity &entity ) { - static std::vector dummy; - return dummy; + static std::vector< double > dummy; + return dummy; } private: - Grid *grid_; + int readVertices ( std::istream &input, std::vector< std::vector< double > > &vertices ) + { + int dimWorld = 3; + dgf::VertexBlock vtxBlock( input, dimWorld ); + if( !vtxBlock.isactive() ) + DUNE_THROW( DGFException, "Vertex block not found" ); + + vtxBlock.get( vertices, vtxParams_, numVtxParams_ ); + return vtxBlock.offset(); + } + + std::vector< std::vector< int > > readPolygons ( std::istream &input, int numVtx, int vtxOfs ) + { + dgf::PolygonBlock polygonBlock( input, numVtx, vtxOfs ); + if( !polygonBlock.isactive() ) + DUNE_THROW( DGFException, "Polygon block not found" ); + + std::vector< std::vector< int > > polygons; + polygonBlock.get( polygons ); + return polygons; + } + + std::vector< std::vector< int > > readPolyhedra ( std::istream &input, int numPolygons ) + { + dgf::PolyhedronBlock polyhedronBlock( input, numPolygons ); + std::vector< std::vector< int > > polyhedra; + if( polyhedronBlock.isactive() ) + { + polyhedronBlock.get( polyhedra ); + } + return polyhedra; + } + + template< class Iterator > + void copy ( Iterator begin, Iterator end, double *dest ) + { + for( ; begin != end; ++begin ) + dest = std::copy( begin->begin(), begin->end(), dest ); + } + + template< class Iterator > + void copy ( Iterator begin, Iterator end, int *dest, int *offset ) + { + int size = 0; + for( ; begin != end; ++begin ) + { + *(offset++) = size; + size += begin->size(); + dest = std::copy( begin->begin(), begin->end(), dest ); + } + *offset = size; + } + + void generate ( std::istream &input ) + { + if( !DuneGridFormatParser::isDuneGridFormat( input ) ) + { + DUNE_THROW( DGFException, "Not in DGF format" ); + } + + typedef std::vector< std::vector< double > > CoordinateVectorType; + CoordinateVectorType nodes; + const int vtxOfs = readVertices( input, nodes ); + + typedef std::vector< std::vector< int > > IndexVectorType; + IndexVectorType faces = readPolygons ( input, nodes.size(), vtxOfs ); + IndexVectorType cells = readPolyhedra( input, faces.size() ); + + if( cells.empty() ) + { + DUNE_THROW( DGFException, "Polyhedron block not found" ); + } + + typedef GridFactory< Grid > GridFactoryType; + typedef typename GridFactoryType :: Coordinate Coordinate ; + + GridFactoryType gridFactory; + + const int nNodes = nodes.size(); + Coordinate node( 0 ); + for( int i=0; i numbers; + + const int nFaces = faces.size(); + for(int i = 0; i < nFaces; ++ i ) + { + // copy values into appropriate data type + std::vector& face = faces[ i ]; + numbers.resize( face.size() ); + std::copy( face.begin(), face.end(), numbers.begin() ); + gridFactory.insertElement( type, numbers ); + } + + // free memory + //faces.swap( IndexVectorType() ); + + // insert cells with type none/dim + type.makeNone( Grid::dimension ); + + const int nCells = cells.size(); + for(int i = 0; i < nCells; ++ i ) + { + // copy values into appropriate data type + std::vector& cell = cells[ i ]; + numbers.resize( cell.size() ); + std::copy( cell.begin(), cell.end(), numbers.begin() ); + gridFactory.insertElement( type, numbers ); + } + //cells.swap( IndexVectorType() ); + + grid_ = gridFactory.createGrid(); + } + + std::unique_ptr< Grid > grid_; + int numVtxParams_; + std::vector< std::vector< double > > vtxParams_; }; // DGFGridInfo for PolyhedralGrid - // ---------------------- + // ------------------------------ template< int dim, int dimworld > struct DGFGridInfo< PolyhedralGrid< dim, dimworld > > { static int refineStepsForHalf () { - return 0; + return 0; } static double refineWeight () { - return 0; + return 0; } }; diff --git a/opm/grid/polyhedralgrid/entity.hh b/opm/grid/polyhedralgrid/entity.hh index ceff39df8..f372dd01c 100644 --- a/opm/grid/polyhedralgrid/entity.hh +++ b/opm/grid/polyhedralgrid/entity.hh @@ -97,7 +97,7 @@ namespace Dune */ GeometryType type () const { - return data()->geomTypes(codimension)[0]; + return data()->cellGeometryType(seed_); } /** \brief obtain the level of this entity */ diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 6b067765d..9bd87e0aa 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -3,7 +3,6 @@ #ifndef DUNE_POLYHEDRALGRID_GRID_HH #define DUNE_POLYHEDRALGRID_GRID_HH -#include #include #include @@ -25,11 +24,13 @@ #include #include #include +#include // Re-enable warnings. #include #include + #include #include #include @@ -38,6 +39,8 @@ namespace Dune { + + // PolyhedralGridFamily // ------------ @@ -66,19 +69,11 @@ namespace Dune typedef PolyhedralGridIntersectionIterator< const Grid > LeafIntersectionIteratorImpl; typedef PolyhedralGridIntersectionIterator< const Grid > LevelIntersectionIteratorImpl; -#if DUNE_VERSION_NEWER(DUNE_GRID,2,3) typedef Dune::Intersection< const Grid, LeafIntersectionImpl > LeafIntersection; typedef Dune::Intersection< const Grid, LevelIntersectionImpl > LevelIntersection; typedef Dune::IntersectionIterator< const Grid, LeafIntersectionIteratorImpl, LeafIntersectionImpl > LeafIntersectionIterator; typedef Dune::IntersectionIterator< const Grid, LevelIntersectionIteratorImpl, LevelIntersectionImpl > LevelIntersectionIterator; -#else - typedef Dune::Intersection< const Grid, PolyhedralGridIntersection > LeafIntersection; - typedef Dune::Intersection< const Grid, PolyhedralGridIntersection > LevelIntersection; - - typedef Dune::IntersectionIterator< const Grid, PolyhedralGridIntersectionIterator, PolyhedralGridIntersection > LeafIntersectionIterator; - typedef Dune::IntersectionIterator< const Grid, PolyhedralGridIntersectionIterator, PolyhedralGridIntersection > LevelIntersectionIterator; -#endif typedef PolyhedralGridIterator< 0, const Grid, All_Partition > HierarchicIteratorImpl; typedef Dune::EntityIterator< 0, const Grid, HierarchicIteratorImpl > HierarchicIterator; @@ -134,13 +129,13 @@ namespace Dune typedef Dune::GridView< PolyhedralGridViewTraits< dim, dimworld, ctype, pitype > > LevelGridView; }; - typedef typename Partition::LevelGridView LevelGridView; - typedef typename Partition::LeafGridView LeafGridView; - + typedef typename Partition< All_Partition >::LeafGridView LeafGridView; + typedef typename Partition< All_Partition >::LevelGridView LevelGridView; }; }; + // PolyhedralGrid // -------------- @@ -159,7 +154,7 @@ namespace Dune < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > > /** \endcond */ { - typedef PolyhedralGrid< dim, dimworld > Grid; + typedef PolyhedralGrid< dim, dimworld, coord_t > Grid; typedef GridDefaultImplementation < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > > Base; @@ -173,7 +168,31 @@ namespace Dune destroy_grid( grdPtr ); } }; + public: + typedef PolyhedralMesh< dim, dimworld, coord_t, int > PolyhedralMeshType; + + typedef std::unique_ptr< UnstructuredGridType, UnstructuredGridDeleter > UnstructuredGridPtr; + + 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 ); + if( !grid ) + DUNE_THROW( GridError, "Unable to allocate grid" ); + return UnstructuredGridPtr( grid ); + } + + static void + computeGeometry ( UnstructuredGridPtr& ug ) + { + // get C pointer to UnstructuredGrid + UnstructuredGrid* ugPtr = ug.operator ->(); + + // compute geometric quantities like cell volume and face normals + compute_geometry( ugPtr ); + } + /** \cond */ typedef PolyhedralGridFamily< dim, dimworld, coord_t > GridFamily; /** \endcond */ @@ -301,7 +320,7 @@ namespace Dune * \param[in] poreVolumes vector with pore volumes (default = empty) */ explicit PolyhedralGrid ( const Opm::Deck& deck, - const std::vector& poreVolumes = std::vector ()) + const std::vector& poreVolumes = std::vector ()) : gridPtr_( createGrid( deck, poreVolumes ) ), grid_( *gridPtr_ ), comm_( *this ), @@ -313,6 +332,24 @@ namespace Dune } #endif + /** \brief constructor + * + * \note The grid will take ownership of the supplied grid pointer. + * + * \param[in] ug pointer to UnstructuredGrid + */ + explicit PolyhedralGrid ( UnstructuredGridPtr &&gridPtr ) + : gridPtr_( std::move( gridPtr ) ), + grid_( *gridPtr_ ), + polyhedralMesh_( grid_ ), + comm_( *this ), + leafIndexSet_( *this ), + globalIdSet_( *this ), + localIdSet_( *this ) + { + init(); + } + /** \brief constructor * * The references to ug are stored in the grid. @@ -375,23 +412,7 @@ namespace Dune */ int size ( int codim ) const { - if( codim == 0 ) - { - return grid_.number_of_cells; - } - else if ( codim == 1 ) - { - return grid_.number_of_faces; - } - else if ( codim == dim ) - { - return grid_.number_of_nodes; - } - else - { - std::cerr << "Warning: codimension " << codim << " not available in PolyhedralGrid" << std::endl; - return 0; - } + return polyhedralMesh_.size( codim ); } /** \brief obtain number of entites on a level @@ -778,6 +799,11 @@ namespace Dune return grid_.global_cell; } + const int* globalCellPtr() const + { + return grid_.global_cell; + } + void getIJK(const int c, std::array& ijk) const { int gc = globalCell()[c]; @@ -787,7 +813,6 @@ namespace Dune } protected: - #if HAVE_ECL_INPUT UnstructuredGridType* createGrid( const Opm::Deck& deck, const std::vector< double >& poreVolumes ) const { @@ -816,21 +841,33 @@ namespace Dune g.actnum = actnum.data(); g.mapaxes = mapaxes.data(); + if (!poreVolumes.empty() && (eclipseGrid->getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) + { + Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]); + const std::vector& minpvv = eclipseGrid->getMinpvVector(); + // Currently the pinchProcessor is not used and only opmfil is supported + // The polyhedralgrid only only supports the opmfil option + //bool opmfil = eclipseGrid->getMinpvMode() == Opm::MinpvMode::OpmFIL; + bool opmfil = true; + const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2]; + std::vector thickness(cartGridSize); + for (size_t i = 0; i < cartGridSize; ++i) { + thickness[i] = eclipseGrid->getCellThicknes(i); + } + const double z_tolerance = eclipseGrid->isPinchActive() ? eclipseGrid->getPinchThresholdThickness() : 0.0; + mp.process(thickness, z_tolerance, poreVolumes, minpvv, actnum, opmfil, zcorn.data()); + } + + /* if (!poreVolumes.empty() && (eclipseGrid->getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) { Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]); - const std::vector& minpvv = eclipseGrid->getMinpvVector(); + const double minpv_value = eclipseGrid->getMinpvValue(); // Currently the pinchProcessor is not used and only opmfil is supported - // The polyhedralgrid only only supports the opmfil option //bool opmfil = eclipseGrid->getMinpvMode() == Opm::MinpvMode::OpmFIL; bool opmfil = true; - const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2]; - std::vector thickness(cartGridSize); - for (size_t i = 0; i < cartGridSize; ++i) { - thickness[i] = eclipseGrid->getCellThicknes(i); - } - const double z_tolerance = eclipseGrid->isPinchActive() ? eclipseGrid->getPinchThresholdThickness() : 0.0; - mp.process(thickness, z_tolerance, poreVolumes, minpvv, actnum, opmfil, zcorn.data()); + mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data()); } + */ const double z_tolerance = eclipseGrid->isPinchActive() ? eclipseGrid->getPinchThresholdThickness() : 0.0; @@ -842,6 +879,7 @@ namespace Dune } #endif + public: using Base::getRealImplementation; @@ -879,19 +917,14 @@ namespace Dune switch (codim) { case 0: - { - const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ]; - return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); - } case 1: { - const int faceVertex = grid_.face_nodes[grid_.face_nodepos[seed.index()] + i]; - return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex ); + const int vxIndex = polyhedralMesh_.subEntity( seed.index(), i, codim ); + return polyhedralMesh_.coordinate( vxIndex ); } case dim: { - const int coordIndex = GlobalCoordinate :: dimension * seed.index(); - return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); + return polyhedralMesh_.coordinate( seed.index() ); } } return GlobalCoordinate( 0 ); @@ -901,15 +934,33 @@ namespace Dune int subEntities( const EntitySeed& seed, const int codim ) const { const int index = seed.index(); - switch (codim) + if( seed.codimension == 0 ) { - case 0: - return 1; - case 1: - return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; - case dim: - return cellVertices_[ index ].size(); + switch (codim) + { + case 0: + return 1; + case 1: + return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; + case 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 ]; + } } + else if ( seed.codimension == dim ) + { + return 1 ; + } + return 0; } @@ -947,6 +998,26 @@ namespace Dune return EntitySeed(); } + template + typename Codim::EntitySeed + subEntitySeed( const typename Codim<1>::EntitySeed& faceSeed, const int i ) const + { + assert( i>= 0 && i::EntitySeed EntitySeed; + if ( codim == 1 ) + { + return EntitySeed( faceSeed.index() ); + } + else if ( codim == dim ) + { + return EntitySeed( grid_.face_nodes[ grid_.face_nodepos[ faceSeed.index() ] + i ] ); + } + else + { + DUNE_THROW(NotImplemented,"codimension not available"); + } + } + const std::vector< GeometryType > &geomTypes ( const unsigned int codim ) const { static std::vector< GeometryType > emptyDummy; @@ -958,6 +1029,18 @@ namespace Dune return emptyDummy; } + template < class Seed > + GeometryType cellGeometryType( const Seed& seed ) const + { + if( Seed::codimension == 0 && ! cellGeomTypes_.empty() ) + { + assert( seed.index() < cellGeomTypes_.size() ); + return cellGeomTypes_[ seed.index() ]; + } + else + return geomTypes( Seed::codimension )[ 0 ]; + } + int indexInInside( const typename Codim<0>::EntitySeed& seed, const int i ) const { return ( grid_.cell_facetag ) ? cartesianIndexInInside( seed, i ) : i; @@ -1044,6 +1127,15 @@ namespace Dune template GlobalCoordinate centroids( const EntitySeed& seed ) const { + if( ! seed.isValid() ) + return GlobalCoordinate( 0 ); + + const int index = seed.index(); + const int codim = EntitySeed::codimension; + assert( index >= 0 && index < size( codim ) ); + return polyhedralMesh_.center( index, codim ); + + /* const int index = GlobalCoordinate :: dimension * seed.index(); const int codim = EntitySeed::codimension; assert( index >= 0 && index < size( codim ) * GlobalCoordinate :: dimension ); @@ -1065,6 +1157,7 @@ namespace Dune DUNE_THROW(InvalidStateException,"codimension not implemented"); return GlobalCoordinate( 0 ); } + */ } GlobalCoordinate copyToGlobalCoordinate( const double* coords ) const @@ -1080,24 +1173,16 @@ namespace Dune template double volumes( const EntitySeed& seed ) const { - const int index = seed.index(); const int codim = EntitySeed::codimension; - if( codim == 0 ) - { - return grid_.cell_volumes[ index ]; - } - else if ( codim == 1 ) - { - return grid_.face_areas[ index ]; - } - else if ( codim == dim ) + if( codim == dim || ! seed.isValid() ) { return 1.0; } else { - DUNE_THROW(InvalidStateException,"codimension not implemented"); - return 0.0; + const int index = seed.index(); + assert( seed.isValid() ); + return polyhedralMesh_.volume( index, codim ); } } @@ -1112,13 +1197,9 @@ namespace Dune // setup list of cell vertices const int numCells = size( 0 ); cellVertices_.resize( numCells ); - - std::cout << numCells << std::endl; - // sort vertices such that they comply with the dune cube reference element if( grid_.cell_facetag ) { - std::cout << "facetag" << std::endl; typedef std::array KeyType; std::map< const KeyType, const int > vertexFaceTags; const int vertexFacePattern [8][3] = { @@ -1143,8 +1224,6 @@ namespace Dune vertexFaceTags.insert( std::make_pair( key, i ) ); } - std::cout << "cells before" << std::endl; - for (int c = 0; c < numCells; ++c) { typedef std::map vertexmap_t; @@ -1172,20 +1251,15 @@ namespace Dune } } } - std::cout << "face node done" << std::endl; - std::cout << cell_pts.size() << std::endl; typedef std::map< int, std::set > vertexlist_t; vertexlist_t vertexList; for( int faceTag = 0; faceTag::max(); - std::cout << "no facetag" << std::endl; for (int c = 0; c < numCells; ++c) { std::set cell_pts; @@ -1244,7 +1315,116 @@ namespace Dune cellVertices_[ c ].resize( cell_pts.size() ); std::copy(cell_pts.begin(), cell_pts.end(), cellVertices_[ c ].begin() ); + maxVx = std::max( maxVx, int( cell_pts.size() ) ); + minVx = std::min( minVx, int( cell_pts.size() ) ); + } + + if( minVx == maxVx && maxVx == 4 ) + { + for (int c = 0; c < numCells; ++c) + { + assert( cellVertices_[ c ].size() == 4 ); + GlobalCoordinate center( 0 ); + GlobalCoordinate p[ 4 ]; + for( int i=0; i<4; ++i ) + { + const int vertex = cellVertices_[ c ][ i ]; + + for( int d=0; d matrix( 0 ); + matrix [0][0] = p[1][0] - p[0][0] ; + matrix [0][1] = p[1][1] - p[0][1] ; + matrix [0][2] = p[1][2] - p[0][2] ; + + matrix [1][0] = p[2][0] - p[0][0] ; + matrix [1][1] = p[2][1] - p[0][1] ; + matrix [1][2] = p[2][2] - p[0][2] ; + + matrix [2][0] = p[3][0] - p[0][0] ; + matrix [2][1] = p[3][1] - p[0][1] ; + matrix [2][2] = p[3][2] - p[0][2] ; + + grid_.cell_volumes[ c ] = std::abs( matrix.determinant() )/ 6.0; + } } + + // check face normals + { + typedef Dune::FieldVector< double, dim > Coordinate; + const int faces = grid_.number_of_faces; + for( int face = 0 ; face < faces; ++face ) + { + const int a = grid_.face_cells[ 2*face ]; + const int b = grid_.face_cells[ 2*face + 1 ]; + + assert( a >=0 || b >=0 ); + + if( grid_.face_areas[ face ] < 0 ) + std::abort(); + + Coordinate centerDiff( 0 ); + if( b >= 0 ) + { + for( int d=0; d= 0 ) + { + for( int d=0; d gridPtr_; + UnstructuredGridPtr gridPtr_; const UnstructuredGridType& grid_; + PolyhedralMeshType polyhedralMesh_; + CollectiveCommunication comm_; std::array< int, 3 > cartDims_; std::vector< std::vector< GeometryType > > geomTypes_; @@ -1285,6 +1476,9 @@ namespace Dune std::vector< GlobalCoordinate > unitOuterNormals_; + // geometry type of each cell if existing (if not then all are polyhedral) + std::vector< GeometryType > cellGeomTypes_; + mutable LeafIndexSet leafIndexSet_; mutable GlobalIdSet globalIdSet_; mutable LocalIdSet localIdSet_; diff --git a/opm/grid/polyhedralgrid/gridfactory.hh b/opm/grid/polyhedralgrid/gridfactory.hh new file mode 100644 index 000000000..c02a1cb6f --- /dev/null +++ b/opm/grid/polyhedralgrid/gridfactory.hh @@ -0,0 +1,264 @@ +// -*- mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=2 sw=2 sts=2: +#ifndef DUNE_POLYHEDRALGRID_GRIDFACTORY_HH +#define DUNE_POLYHEDRALGRID_GRIDFACTORY_HH + +#include +#include + +#include +#include + +#include +#include + +namespace Dune +{ + + + // GridFactory for PolyhedralGrid + // --------------------------------- + + template< int dim, int dimworld > + class GridFactory< PolyhedralGrid< dim, dimworld > > + : public GridFactoryInterface< PolyhedralGrid< dim, dimworld > > + { + public: + typedef PolyhedralGrid< dim, dimworld > Grid; + + const static int dimension = Grid::dimension; + const static int dimensionworld = Grid::dimensionworld; + typedef typename Grid::ctype ctype; + + typedef MPIHelper::MPICommunicator MPICommunicatorType; + typedef typename Grid::template Codim<0>::Entity Element; + typedef typename Grid::template Codim::Entity Vertex; + + typedef Dune::FieldVector CoordinateType; + typedef CoordinateType Coordinate; + +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) + typedef ToUniquePtr UniquePtrType; +#else // #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) + typedef Grid* UniquePtrType; +#endif // #else // #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) + + + /** \brief Default constructor */ + explicit GridFactory ( const MPICommunicatorType &communicator = MPIHelper::getCommunicator() ) + : nodes_(), + faces_(), + cells_() + {} + + virtual void insertVertex(const CoordinateType& pos) + { + nodes_.push_back( pos ); + } + + /** \brief Insert an element into the coarse grid + \param type The GeometryType of the new element + \param items The items are usually the vertex numbers of the inserted + element. If the geometry type is none the these can be face numbers. + + \note If the GeometryType is none then faces need to be inserted separately + using this method and passing a GeometryType with dimension - 1 + (with respect to the Grid's dimension). + */ + virtual void insertElement(const GeometryType& type, + const std::vector& items) + { + if( type.isNone() ) + { + // copy into vector of integers + std::vector< int > numbers( items.size() ); + std::copy( items.begin(), items.end(), numbers.begin() ); + + if( type.dim() == dimension-1 ) + { + faces_.push_back( numbers ); + } + else if( type.dim() == dimension ) + { + // note vertices holds the face + // numbers in this case + cells_.push_back( numbers ); + } + else + { + DUNE_THROW(Dune::NotImplemented,"insertElement not implemented for type " << type ); + } + } + else // use ReferenceElement to insert faces + { + + + } + } + + virtual void insertElement(const GeometryType& type, + const std::vector& vertices, + const shared_ptr,FieldVector > >& elementParametrization) + { + std::cerr << "Warning: elementParametrization is being ignored in insertElement!" << std::endl; + insertElement( type, vertices ); + } + + void insertBoundarySegment(const std::vector& vertices) + { + DUNE_THROW(NotImplemented,"yet"); + } + + UniquePtrType createGrid() + { + std::vector< CoordinateType >& nodes = nodes_; + std::vector< std::vector< int > >& faces = faces_; + std::vector< std::vector< int > >& cells = cells_; + + if( cells.empty() ) + { + DUNE_THROW( GridError, "No cells found for PolyhedralGrid" ); + } + + const auto sumSize = [] ( std::size_t s, const std::vector< int > &v ) { return s + v.size(); }; + const std::size_t numFaceNodes = std::accumulate( faces.begin(), faces.end(), std::size_t( 0 ), sumSize ); + const std::size_t numCellFaces = std::accumulate( cells.begin(), cells.end(), std::size_t( 0 ), sumSize ); + + typename Grid::UnstructuredGridPtr ug = + Grid::allocateGrid( cells.size(), faces.size(), numFaceNodes, numCellFaces, nodes.size() ); + + // copy faces + { +#ifndef NDEBUG + std::map< std::vector< int >, std::vector< int > > faceMap; +#endif + + const int nFaces = faces.size(); + // set all face_cells values to -2 as default + std::fill( ug->face_cells, ug->face_cells + 2*nFaces, -1 ); + + int facepos = 0; + std::vector< int > faceVertices; + faceVertices.reserve( 30 ); + for( int face = 0; face < nFaces; ++face ) + { + //std::cout << "face " << face << ": "; + faceVertices.clear(); + ug->face_nodepos[ face ] = facepos; + const int nVertices = faces[ face ].size(); + for( int vx = 0; vx < nVertices; ++vx, ++facepos ) + { + //std::cout << " " << faces[ face ][ vx ]; + ug->face_nodes[ facepos ] = faces[ face ][ vx ]; + faceVertices.push_back( faces[ face ][ vx ] ); + } + //std::cout << std::endl; + +#ifndef NDEBUG + // sort vertices + std::sort( faceVertices.begin(), faceVertices.end() ); + // make sure each face only exists once + faceMap[ faceVertices ].push_back( face ); + assert( faceMap[ faceVertices ].size() == 1 ); +#endif + } + ug->face_nodepos[ nFaces ] = facepos ; + } + + // copy cells + { + const int nCells = cells.size(); + int cellpos = 0; + for( int cell = 0; cell < nCells; ++cell ) + { + //std::cout << "Cell " << cell << ": "; + ug->cell_facepos[ cell ] = cellpos; + const int nFaces = cells[ cell ].size(); + for( int f = 0; f < nFaces; ++f, ++cellpos ) + { + const int face = cells[ cell ][ f ]; + // std::cout << " " << face ; + ug->cell_faces[ cellpos ] = face; + + // TODO find cells for each face + if( ug->face_cells[ 2*face ] == -1 ) + { + ug->face_cells[ 2*face ] = cell; + } + else // if ( ug->face_cells[ 2*face+1 ] == -1 ) + { + //assert( ug->face_cells[ 2*face+1 ] == -1 ); + ug->face_cells[ 2*face+1 ] = cell; + } + } + //std::cout << std::endl; + } + ug->cell_facepos[ nCells ] = cellpos ; + } + + // copy node coordinates + { + const int nNodes = nodes.size(); + int nodepos = 0; + for( int vx = 0 ; vx < nNodes; ++vx ) + { + for( int d=0; dnode_coordinates[ nodepos ] = nodes[ vx ][ d ]; + } + } + + /* + for( int i=0; iface_cells[ 2*i ] << " " << + ug->face_cells[ 2*i+1] << std::endl; + } + */ + + // free cell face tag since it's not a cartesian grid + if( ug->cell_facetag ) + { + std::free( ug->cell_facetag ); + ug->cell_facetag = nullptr ; + for( int i=0; i<3; ++i ) ug->cartdims[ i ] = 0; + } + + // compute geometric quntities like cell volume and face normals + Grid::computeGeometry( ug ); + + // check normal direction + { + const int faces = ug->number_of_faces; + for( int face = 0 ; face < faces; ++face ) + { + const int a = ug->face_cells[ 2*face ]; + const int b = ug->face_cells[ 2*face + 1 ]; + Coordinate centerDiff( 0 ); + Coordinate normal( 0 ); + for( int d=0; dcell_centroids[ b*dim + d ] - ug->cell_centroids[ a*dim + d ]; + normal[ d ] = ug->face_normals[ face*dim + d ]; + } + + // if diff and normal point in different direction, flip faces + if( centerDiff * normal > 0 ) + { + ug->face_cells[ 2*face ] = b; + ug->face_cells[ 2*face + 1 ] = a; + } + } + } + + return new Grid( std::move( ug ) ); + } + + protected: + std::vector< CoordinateType > nodes_; + std::vector< std::vector< int > > faces_; + std::vector< std::vector< int > > cells_; + }; + +} // namespace Dune + +#endif // #ifndef DUNE_POLYHEDRALGRID_DGFPARSER_HH diff --git a/opm/grid/polyhedralgrid/idset.hh b/opm/grid/polyhedralgrid/idset.hh index 57b1672fe..47928e949 100644 --- a/opm/grid/polyhedralgrid/idset.hh +++ b/opm/grid/polyhedralgrid/idset.hh @@ -24,7 +24,8 @@ namespace Dune typedef IdSet< Grid, This, IdType > Base; PolyhedralGridIdSet (const Grid& grid) - : grid_(grid) + : grid_( grid ), + globalCellPtr_( grid_.globalCellPtr() ) {} //! id meethod for entity and specific codim @@ -32,12 +33,22 @@ namespace Dune IdType id ( const typename Traits::template Codim< codim >::Entity &entity ) const { const int index = entity.seed().index(); - if (codim == 0) - return grid_.globalCell()[ index ]; + // in case + if (codim == 0 && globalCellPtr_ ) + return globalCellPtr_[ index ]; else return index; } +#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 @@ -59,9 +70,9 @@ namespace Dune if( codim == 0 ) return id( entity ); else if ( codim == 1 ) - return id( Grid::getRealImplementation( entity ).template subEntity< 1 > ( i ) ); + return id( entity.template subEntity< 1 >( i ) ); else if ( codim == dim ) - return id( Grid::getRealImplementation( entity ).template subEntity< dim > ( i ) ); + return id( entity.template subEntity< dim >( i ) ); else { DUNE_THROW(NotImplemented,"codimension not available"); @@ -71,6 +82,7 @@ namespace Dune protected: const Grid& grid_; + const int* globalCellPtr_; }; } // namespace Dune diff --git a/opm/grid/polyhedralgrid/indexset.hh b/opm/grid/polyhedralgrid/indexset.hh index 34d745a29..e542f8a49 100644 --- a/opm/grid/polyhedralgrid/indexset.hh +++ b/opm/grid/polyhedralgrid/indexset.hh @@ -52,6 +52,20 @@ namespace Dune return grid().getRealImplementation(entity).index(); } +#if ! DUNE_VERSION_NEWER(DUNE_GRID,2,4) + template< int cd > + IndexType index ( const typename Traits::template Codim< cd >::EntityPointer &entityPointer ) const + { + return index( *entityPointer ); + } +#endif + + template< int cd > + IndexType subIndex ( const typename Traits::template Codim< cd >::Entity &entity, int i, unsigned int codim ) const + { + return subIndex( entity, i, codim ); + } + template< class Entity > IndexType subIndex ( const Entity &entity, int i, unsigned int codim ) const { diff --git a/opm/grid/polyhedralgrid/polyhedralmesh.hh b/opm/grid/polyhedralgrid/polyhedralmesh.hh new file mode 100644 index 000000000..8c41d337b --- /dev/null +++ b/opm/grid/polyhedralgrid/polyhedralmesh.hh @@ -0,0 +1,244 @@ +#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 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/tests/grid_test.cc b/tests/grid_test.cc index b7fc086d9..1b0e025d0 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -185,18 +185,17 @@ int main(int argc, char** argv ) //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 før " << std::endl; + std::cout << "tsDune for " << std::endl; Grid2D tsDune (*ts); - std::cout << "tsDune etter " << std::endl; + std::cout << "tsDune after " << std::endl; testGrid ( tsDune, "ts"); #endif - - - //Dune::GridPtr< Grid > gridPtr( dgfFile ); - //testGrid( *gridPtr, "polyhedralgrid-dgf" ); + Dune::GridPtr< Grid > gridPtr( dgfFile ); + testGrid( *gridPtr, "polyhedralgrid-dgf" ); } + /* // test CpGrid { typedef Dune::CpGrid Grid; @@ -212,8 +211,9 @@ int main(int argc, char** argv ) testGrid( grid, "cpgrid2" ); #endif - Dune::GridPtr< Grid > gridPtr( dgfFile ); + //Dune::GridPtr< Grid > gridPtr( dgfFile ); //testGrid( *gridPtr, "cpgrid-dgf" ); } + */ return 0; } From 9834617973f29c625081dc1397dbf6d543cc5b1a Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 15 Mar 2019 15:57:36 +0100 Subject: [PATCH 04/22] Implemented DGF for normal grids. --- opm/grid/polyhedralgrid/dgfparser.hh | 114 +++++++++++++++++++++++++-- opm/grid/polyhedralgrid/grid.hh | 4 +- tests/grid_test.cc | 6 +- 3 files changed, 114 insertions(+), 10 deletions(-) diff --git a/opm/grid/polyhedralgrid/dgfparser.hh b/opm/grid/polyhedralgrid/dgfparser.hh index 9dd9ff764..9fb8e63c6 100644 --- a/opm/grid/polyhedralgrid/dgfparser.hh +++ b/opm/grid/polyhedralgrid/dgfparser.hh @@ -9,6 +9,8 @@ #include #include +#include + #include #include @@ -217,7 +219,7 @@ namespace Dune private: int readVertices ( std::istream &input, std::vector< std::vector< double > > &vertices ) { - int dimWorld = 3; + int dimWorld = Grid::dimensionworld ; dgf::VertexBlock vtxBlock( input, dimWorld ); if( !vtxBlock.isactive() ) DUNE_THROW( DGFException, "Vertex block not found" ); @@ -270,22 +272,120 @@ namespace Dune void generate ( std::istream &input ) { - if( !DuneGridFormatParser::isDuneGridFormat( input ) ) + DuneGridFormatParser dgf( 0, 1 ); + + dgf.element = DuneGridFormatParser::General; + dgf.dimgrid = dim; + dgf.dimw = dimworld; + + if( !dgf.isDuneGridFormat( input ) ) { DUNE_THROW( DGFException, "Not in DGF format" ); } + if( !dgf.readDuneGrid( input, dim, dimworld ) ) + DUNE_THROW( InvalidStateException, "DGF file not recognized on second call." ); + typedef std::vector< std::vector< double > > CoordinateVectorType; CoordinateVectorType nodes; - const int vtxOfs = readVertices( input, nodes ); typedef std::vector< std::vector< int > > IndexVectorType; - IndexVectorType faces = readPolygons ( input, nodes.size(), vtxOfs ); - IndexVectorType cells = readPolyhedra( input, faces.size() ); + IndexVectorType faces; + IndexVectorType cells; - if( cells.empty() ) + if( dgf.nofelements == 0 ) { - DUNE_THROW( DGFException, "Polyhedron block not found" ); + const int vtxOfs = readVertices( input, nodes ); + + faces = readPolygons ( input, nodes.size(), vtxOfs ); + cells = readPolyhedra( input, faces.size() ); + + if( cells.empty() ) + { + DUNE_THROW( DGFException, "Polyhedron block not found" ); + } + } + else // convert dgf input to polyhedral + { + nodes.resize( dgf.nofvtx ); + // copy vertices + std::copy( dgf.vtx.begin(), dgf.vtx.end(), nodes.begin() ); + + for( const auto& node : nodes ) + { + for( size_t i=0; i face_t; + std::map< face_t, int > tmpFaces; + + const int nFaces = (nVx == dim+1) ? dim+1 : 2*dim; + + Dune::GeometryType type( (nVx == dim+1) ? + Impl :: SimplexTopology< dim > :: type :: id : + Impl :: CubeTopology < dim > :: type :: id, + dim ); + + const auto& refElem = Dune::referenceElement< typename Grid::ctype, dim >( type ); + + cells.resize( dgf.nofelements ); + + face_t face; + int faceNo = 0; + for( int n = 0; n < dgf.nofelements; ++n ) + { + const auto& elem = dgf.elements[ n ]; + auto& cell = cells[ n ]; + assert( elem.size() == nVx ); + cell.resize( nFaces ); + for(int f=0; fsecond; + + cell[ f ] = myFaceNo; + } + + /* + for( const auto& c : cell ) + { + std::cout << c << " "; + } + std::cout << std::endl; + */ + } + + faces.clear(); + faces.resize( tmpFaces.size() ); + for( const auto& face : tmpFaces ) + { + faces[ face.second ] = face.first ; + /* + for( const auto& f : face.first ) + { + std::cout << f << " "; + } + std::cout << std::endl; + */ + } } typedef GridFactory< Grid > GridFactoryType; diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 9bd87e0aa..b90f0b57d 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1410,8 +1410,8 @@ namespace Dune normal[ d ] = grid_.face_normals[ face*dim + d ]; } - if( normal.two_norm() < 1e-10 ) - std::abort(); + //if( normal.two_norm() < 1e-10 ) + // std::abort(); if( centerDiff.two_norm() < 1e-10 ) std::abort(); diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 1b0e025d0..5204cf6d3 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -163,7 +163,9 @@ int main(int argc, char** argv ) dgfFile << "Interval" << std::endl; dgfFile << "0 0 0" << std::endl; dgfFile << "1 1 1" << std::endl; - dgfFile << "8 8 8" << std::endl; + dgfFile << "2 2 2" << std::endl; + dgfFile << "#" << std::endl; + dgfFile << "Simplex" << std::endl; dgfFile << "#" << std::endl; #if HAVE_ECL_INPUT @@ -175,6 +177,7 @@ int main(int argc, char** argv ) // test PolyhedralGrid { typedef Dune::PolyhedralGrid< 3, 3 > Grid; + /* #if HAVE_ECL_INPUT Grid grid(deck, porv); testGrid( grid, "polyhedralgrid" ); @@ -191,6 +194,7 @@ int main(int argc, char** argv ) testGrid ( tsDune, "ts"); #endif + */ Dune::GridPtr< Grid > gridPtr( dgfFile ); testGrid( *gridPtr, "polyhedralgrid-dgf" ); } From 73b299e97eed93be2df8e3c574772cb7c2c04fbe Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Wed, 24 Apr 2019 16:43:18 +0200 Subject: [PATCH 05/22] [bugfix][DGFFactory] release grid pointer once accessed from GridPtr to avoid false destruction. --- opm/grid/polyhedralgrid/dgfparser.hh | 281 +++++++++++----------- opm/grid/polyhedralgrid/grid.hh | 135 ++++++++++- opm/grid/polyhedralgrid/polyhedralmesh.hh | 9 + tests/grid_test.cc | 6 +- 4 files changed, 277 insertions(+), 154 deletions(-) diff --git a/opm/grid/polyhedralgrid/dgfparser.hh b/opm/grid/polyhedralgrid/dgfparser.hh index 9fb8e63c6..5660f4272 100644 --- a/opm/grid/polyhedralgrid/dgfparser.hh +++ b/opm/grid/polyhedralgrid/dgfparser.hh @@ -5,6 +5,8 @@ #include #include +#include +#include #include #include @@ -20,7 +22,7 @@ #include -#if HAVE_OPM_PARSER +#if HAVE_ECL_INPUT #include #include #include @@ -145,7 +147,8 @@ namespace Dune typedef typename Grid::template Codim::Entity Vertex; explicit DGFGridFactory ( std::istream &input, MPICommunicator comm = MPIHelper::getCommunicator() ) - : grid_() + : gridPtr_(), + grid_( nullptr ) { input.clear(); input.seekg( 0 ); @@ -155,21 +158,21 @@ namespace Dune } explicit DGFGridFactory ( const std::string &filename, MPICommunicator comm = MPIHelper::getCommunicator() ) - : grid_() + : gridPtr_(), + grid_( nullptr ) { std::ifstream input( filename ); if( !input ) DUNE_THROW( DGFException, "Macrofile '" << filename << "' not found" ); -#if HAVE_OPM_PARSER +#if HAVE_ECL_INPUT if( !DuneGridFormatParser::isDuneGridFormat( input ) ) { Opm::Parser parser; - Opm::ParseContext parseContext; - const auto deck = parser.parseFile(filename, parseContext); + const auto deck = parser.parseString( filename ); std::vector porv; - grid_.reset( new Grid( deck, porv ) ); + gridPtr_.reset( new Grid( deck, porv ) ); return ; } else @@ -179,7 +182,15 @@ namespace Dune } } - Grid *grid () const { return grid_.operator->(); } + Grid *grid () const + { + if( ! grid_ ) + { + // set pointer to grid and avoid grid being deleted + grid_ = gridPtr_.release(); + } + return grid_; + } template< class Intersection > bool wasInserted ( const Intersection &intersection ) const @@ -272,29 +283,33 @@ namespace Dune void generate ( std::istream &input ) { - DuneGridFormatParser dgf( 0, 1 ); - - dgf.element = DuneGridFormatParser::General; - dgf.dimgrid = dim; - dgf.dimw = dimworld; - - if( !dgf.isDuneGridFormat( input ) ) + // check whether an interval block is active, otherwise read polyhedrons + dgf::IntervalBlock intervalBlock( input ); + if( intervalBlock.isactive() ) { - DUNE_THROW( DGFException, "Not in DGF format" ); - } + if( intervalBlock.numIntervals() != 1 ) + DUNE_THROW( DGFException, "Currently, CpGrid can only handle 1 interval block." ); - if( !dgf.readDuneGrid( input, dim, dimworld ) ) - DUNE_THROW( InvalidStateException, "DGF file not recognized on second call." ); + if( intervalBlock.dimw() != dimworld ) + DUNE_THROW( DGFException, "CpGrid cannot handle an interval of dimension " << intervalBlock.dimw() << "." ); + const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 ); - typedef std::vector< std::vector< double > > CoordinateVectorType; - CoordinateVectorType nodes; + std::vector< double > spacing( dimworld ); + for( int i=0; i > IndexVectorType; - IndexVectorType faces; - IndexVectorType cells; - - if( dgf.nofelements == 0 ) + gridPtr_.reset( new Grid( interval.n, spacing ) ); + return ; + } + else // polyhedral input { + typedef std::vector< std::vector< double > > CoordinateVectorType; + CoordinateVectorType nodes; + + typedef std::vector< std::vector< int > > IndexVectorType; + IndexVectorType faces; + IndexVectorType cells; + const int vtxOfs = readVertices( input, nodes ); faces = readPolygons ( input, nodes.size(), vtxOfs ); @@ -302,144 +317,134 @@ namespace Dune if( cells.empty() ) { - DUNE_THROW( DGFException, "Polyhedron block not found" ); + DUNE_THROW( DGFException, "Polyhedron block not found!" ); } - } - else // convert dgf input to polyhedral - { - nodes.resize( dgf.nofvtx ); - // copy vertices - std::copy( dgf.vtx.begin(), dgf.vtx.end(), nodes.begin() ); - for( const auto& node : nodes ) - { - for( size_t i=0; i face_t; - std::map< face_t, int > tmpFaces; + typedef GridFactory< Grid > GridFactoryType; + typedef typename GridFactoryType :: Coordinate Coordinate ; - const int nFaces = (nVx == dim+1) ? dim+1 : 2*dim; + GridFactoryType gridFactory; - Dune::GeometryType type( (nVx == dim+1) ? - Impl :: SimplexTopology< dim > :: type :: id : - Impl :: CubeTopology < dim > :: type :: id, - dim ); + const int nNodes = nodes.size(); + Coordinate node( 0 ); + for( int i=0; i( type ); + gridFactory.insertVertex( node ); + } + //nodes.swap( CoordinateVectorType() ); - cells.resize( dgf.nofelements ); + // insert faces with type none/dim-1 + GeometryType type; + type.makeNone( Grid::dimension - 1 ); + std::vector< unsigned int > numbers; - face_t face; - int faceNo = 0; - for( int n = 0; n < dgf.nofelements; ++n ) + const int nFaces = faces.size(); + for(int i = 0; i < nFaces; ++ i ) { - const auto& elem = dgf.elements[ n ]; - auto& cell = cells[ n ]; - assert( elem.size() == nVx ); - cell.resize( nFaces ); - for(int f=0; fsecond; + // copy values into appropriate data type + std::vector& face = faces[ i ]; + numbers.resize( face.size() ); + std::copy( face.begin(), face.end(), numbers.begin() ); + gridFactory.insertElement( type, numbers ); + } - cell[ f ] = myFaceNo; - } + // free memory + //faces.swap( IndexVectorType() ); - /* - for( const auto& c : cell ) - { - std::cout << c << " "; - } - std::cout << std::endl; - */ + // insert cells with type none/dim + type.makeNone( Grid::dimension ); + + const int nCells = cells.size(); + for(int i = 0; i < nCells; ++ i ) + { + // copy values into appropriate data type + std::vector& cell = cells[ i ]; + numbers.resize( cell.size() ); + std::copy( cell.begin(), cell.end(), numbers.begin() ); + gridFactory.insertElement( type, numbers ); } + //cells.swap( IndexVectorType() ); - faces.clear(); - faces.resize( tmpFaces.size() ); - for( const auto& face : tmpFaces ) + gridPtr_ = gridFactory.createGrid(); + } // end else branch + + // alternative conversion to polyhedral format that does not work yet. +#if 0 + else // convert dgf input to polyhedral { - faces[ face.second ] = face.first ; - /* - for( const auto& f : face.first ) + nodes.resize( dgf.nofvtx ); + // copy vertices + std::copy( dgf.vtx.begin(), dgf.vtx.end(), nodes.begin() ); + + for( const auto& node : nodes ) { - std::cout << f << " "; + for( size_t i=0; i GridFactoryType; - typedef typename GridFactoryType :: Coordinate Coordinate ; - GridFactoryType gridFactory; + const unsigned int nVx = dgf.elements[ 0 ].size(); - const int nNodes = nodes.size(); - Coordinate node( 0 ); - for( int i=0; i face_t; + std::map< face_t, int > tmpFaces; - gridFactory.insertVertex( node ); - } - //nodes.swap( CoordinateVectorType() ); + const int nFaces = (nVx == dim+1) ? dim+1 : 2*dim; - // insert faces with type none/dim-1 - GeometryType type; - type.makeNone( Grid::dimension - 1 ); - std::vector< unsigned int > numbers; + Dune::GeometryType type( (nVx == dim+1) ? + Impl :: SimplexTopology< dim > :: type :: id : + Impl :: CubeTopology < dim > :: type :: id, + dim ); - const int nFaces = faces.size(); - for(int i = 0; i < nFaces; ++ i ) - { - // copy values into appropriate data type - std::vector& face = faces[ i ]; - numbers.resize( face.size() ); - std::copy( face.begin(), face.end(), numbers.begin() ); - gridFactory.insertElement( type, numbers ); - } + const auto& refElem = Dune::referenceElement< typename Grid::ctype, dim >( type ); - // free memory - //faces.swap( IndexVectorType() ); + cells.resize( dgf.nofelements ); - // insert cells with type none/dim - type.makeNone( Grid::dimension ); + face_t face; + int faceNo = 0; + for( int n = 0; n < dgf.nofelements; ++n ) + { + const auto& elem = dgf.elements[ n ]; + auto& cell = cells[ n ]; + assert( elem.size() == nVx ); + cell.resize( nFaces ); + for(int f=0; fsecond; - const int nCells = cells.size(); - for(int i = 0; i < nCells; ++ i ) - { - // copy values into appropriate data type - std::vector& cell = cells[ i ]; - numbers.resize( cell.size() ); - std::copy( cell.begin(), cell.end(), numbers.begin() ); - gridFactory.insertElement( type, numbers ); - } - //cells.swap( IndexVectorType() ); + cell[ f ] = myFaceNo; + } - grid_ = gridFactory.createGrid(); + /* + for( const auto& c : cell ) + { + std::cout << c << " "; + } + std::cout << std::endl; + */ + } +#endif } - std::unique_ptr< Grid > grid_; + mutable std::unique_ptr< Grid > gridPtr_; + mutable Grid* grid_; int numVtxParams_; std::vector< std::vector< double > > vtxParams_; }; diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index b90f0b57d..40da066a4 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -32,6 +32,7 @@ #include #include +#include #include #include #include @@ -159,8 +160,10 @@ namespace Dune typedef GridDefaultImplementation < dim, dimworld, coord_t, PolyhedralGridFamily< dim, dimworld, coord_t > > Base; + public: typedef UnstructuredGrid UnstructuredGridType; + protected: struct UnstructuredGridDeleter { inline void operator () ( UnstructuredGridType* grdPtr ) @@ -332,6 +335,23 @@ namespace Dune } #endif + /** \brief constructor + * + * \param[in] deck Opm Eclipse deck + * \param[in] poreVolumes vector with pore volumes (default = empty) + */ + explicit PolyhedralGrid ( const std::vector< int >& n, + const std::vector< double >& dx ) + : gridPtr_( createGrid( n, dx ) ), + grid_( *gridPtr_ ), + comm_( *this ), + leafIndexSet_( *this ), + globalIdSet_( *this ), + localIdSet_( *this ) + { + init(); + } + /** \brief constructor * * \note The grid will take ownership of the supplied grid pointer. @@ -341,7 +361,7 @@ namespace Dune explicit PolyhedralGrid ( UnstructuredGridPtr &&gridPtr ) : gridPtr_( std::move( gridPtr ) ), grid_( *gridPtr_ ), - polyhedralMesh_( grid_ ), + //polyhedralMesh_( grid_ ), comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), @@ -412,7 +432,24 @@ namespace Dune */ int size ( int codim ) const { - return polyhedralMesh_.size( codim ); + //return polyhedralMesh_.size( codim ); + if( codim == 0 ) + { + return grid_.number_of_cells; + } + else if ( codim == 1 ) + { + return grid_.number_of_faces; + } + else if ( codim == dim ) + { + return grid_.number_of_nodes; + } + else + { + std::cerr << "Warning: codimension " << codim << " not available in PolyhedralGrid" << std::endl; + return 0; + } } /** \brief obtain number of entites on a level @@ -879,6 +916,24 @@ namespace Dune } #endif + UnstructuredGridType* createGrid( const std::vector< int >& n, const std::vector< double >& dx ) const + { + UnstructuredGridType* cgrid = nullptr ; + assert( int(n.size()) == dim ); + if( dim == 2 ) + { + cgrid = create_grid_cart2d( n[ 0 ], n[ 1 ], dx[ 0 ], dx[ 1 ] ); + } + else if ( dim == 3 ) + { + cgrid = create_grid_hexa3d( n[ 0 ], n[ 1 ], n[ 2 ], dx[ 0 ], dx[ 1 ], dx[ 2 ] ); + } + + if (!cgrid) { + OPM_THROW(std::runtime_error, "Failed to construct grid."); + } + return cgrid; + } public: using Base::getRealImplementation; @@ -899,7 +954,7 @@ namespace Dune } case 1: { - return 0;//grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; + return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; } case dim: { @@ -913,6 +968,27 @@ namespace Dune GlobalCoordinate corner( const EntitySeed& seed, const int i ) const { + const int codim = EntitySeed :: codimension; + switch (codim) + { + case 0: + { + const int coordIndex = GlobalCoordinate :: dimension * cellVertices_[ seed.index() ][ i ]; + return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); + } + case 1: + { + const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + i]; + return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex ); + } + case dim: + { + const int coordIndex = GlobalCoordinate :: dimension * seed.index(); + return copyToGlobalCoordinate( grid_.node_coordinates + coordIndex ); + } + } + + /* const int codim = EntitySeed :: codimension; switch (codim) { @@ -927,6 +1003,7 @@ namespace Dune return polyhedralMesh_.coordinate( seed.index() ); } } + */ return GlobalCoordinate( 0 ); } @@ -1130,12 +1207,6 @@ namespace Dune if( ! seed.isValid() ) return GlobalCoordinate( 0 ); - const int index = seed.index(); - const int codim = EntitySeed::codimension; - assert( index >= 0 && index < size( codim ) ); - return polyhedralMesh_.center( index, codim ); - - /* const int index = GlobalCoordinate :: dimension * seed.index(); const int codim = EntitySeed::codimension; assert( index >= 0 && index < size( codim ) * GlobalCoordinate :: dimension ); @@ -1157,7 +1228,6 @@ namespace Dune DUNE_THROW(InvalidStateException,"codimension not implemented"); return GlobalCoordinate( 0 ); } - */ } GlobalCoordinate copyToGlobalCoordinate( const double* coords ) const @@ -1173,19 +1243,36 @@ namespace Dune template double volumes( const EntitySeed& seed ) const { - const int codim = EntitySeed::codimension; + static const int codim = EntitySeed::codimension; if( codim == dim || ! seed.isValid() ) { return 1.0; } else { - const int index = seed.index(); assert( seed.isValid() ); - return polyhedralMesh_.volume( index, codim ); + + if( codim == 0 ) + { + return grid_.cell_volumes[ seed.index() ]; + } + else if ( codim == 1 ) + { + return grid_.face_areas[ seed.index() ]; + } + else + { + DUNE_THROW(InvalidStateException,"codimension not implemented"); + return 0.0; + } + + //const int index = seed.index(); + //assert( seed.isValid() ); + //return polyhedralMesh_.volume( index, codim ); } } + protected: void init() { // copy Cartesian dimensions @@ -1196,6 +1283,7 @@ namespace Dune // setup list of cell vertices const int numCells = size( 0 ); + cellVertices_.resize( numCells ); // sort vertices such that they comply with the dune cube reference element if( grid_.cell_facetag ) @@ -1461,8 +1549,29 @@ namespace Dune unitOuterNormals_[ face ] = normal; } + + //print( std::cout, grid_ ); + } + + void print( std::ostream& out, const UnstructuredGridType& grid ) const + { + const int numCells = grid.number_of_cells; + for( int c=0; c Date: Mon, 6 May 2019 11:28:53 +0200 Subject: [PATCH 06/22] Re-enable verteq test. --- tests/grid_test.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 744f12515..e34f5e1e4 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -177,7 +177,6 @@ int main(int argc, char** argv ) // test PolyhedralGrid { typedef Dune::PolyhedralGrid< 3, 3 > Grid; - /* #if HAVE_ECL_INPUT Grid grid(deck, porv); testGrid( grid, "polyhedralgrid" ); @@ -194,7 +193,6 @@ int main(int argc, char** argv ) testGrid ( tsDune, "ts"); #endif - */ Dune::GridPtr< Grid > gridPtr( dgfFile ); testGrid( *gridPtr, "polyhedralgrid-dgf" ); } From 791667a9480e35e5a4d620a76da59c3444e46de9 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Mon, 6 May 2019 11:31:52 +0200 Subject: [PATCH 07/22] Remove debug output. --- opm/grid/polyhedralgrid/grid.hh | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index e17fced0d..f55e3af4c 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1289,7 +1289,6 @@ namespace Dune // sort vertices such that they comply with the dune cube reference element if( grid_.cell_facetag ) { - std::cout << "facetag" << std::endl; typedef std::array KeyType; std::map< const KeyType, const int > vertexFaceTags; const int vertexFacePattern [8][3] = { @@ -1314,8 +1313,6 @@ namespace Dune vertexFaceTags.insert( std::make_pair( key, i ) ); } - std::cout << "cells before" << std::endl; - for (int c = 0; c < numCells; ++c) { typedef std::map vertexmap_t; @@ -1343,15 +1340,12 @@ namespace Dune } } } - std::cout << "face node done" << std::endl; - std::cout << cell_pts.size() << std::endl; typedef std::map< int, std::set > vertexlist_t; vertexlist_t vertexList; for( int faceTag = 0; faceTag Date: Tue, 7 May 2019 15:20:27 +0200 Subject: [PATCH 08/22] [bugfix][TopSurf] include correct header. --- opm/grid/verteq/topsurf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/grid/verteq/topsurf.hpp b/opm/grid/verteq/topsurf.hpp index 0be3b3037..70f0cc11b 100644 --- a/opm/grid/verteq/topsurf.hpp +++ b/opm/grid/verteq/topsurf.hpp @@ -5,7 +5,7 @@ // This file is licensed under the GNU General Public License v3.0 #ifndef OPM_GRID_HEADER_INCLUDED -#include +#include #endif namespace Opm { From 20968e526ecc0fab8b0b4c61e756d31e2116d892 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Mon, 13 May 2019 14:18:57 +0200 Subject: [PATCH 09/22] [bugfix][GridFactory] make code compile with DUNE 2.6 --- opm/grid/polyhedralgrid/dgfparser.hh | 5 +++++ opm/grid/polyhedralgrid/gridfactory.hh | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/opm/grid/polyhedralgrid/dgfparser.hh b/opm/grid/polyhedralgrid/dgfparser.hh index 5660f4272..99366e816 100644 --- a/opm/grid/polyhedralgrid/dgfparser.hh +++ b/opm/grid/polyhedralgrid/dgfparser.hh @@ -368,7 +368,12 @@ namespace Dune } //cells.swap( IndexVectorType() ); +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7) gridPtr_ = gridFactory.createGrid(); +#else + gridPtr_.reset( gridFactory.createGrid() ); +#endif + } // end else branch // alternative conversion to polyhedral format that does not work yet. diff --git a/opm/grid/polyhedralgrid/gridfactory.hh b/opm/grid/polyhedralgrid/gridfactory.hh index c02a1cb6f..10ee16e19 100644 --- a/opm/grid/polyhedralgrid/gridfactory.hh +++ b/opm/grid/polyhedralgrid/gridfactory.hh @@ -37,7 +37,7 @@ namespace Dune typedef Dune::FieldVector CoordinateType; typedef CoordinateType Coordinate; -#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7) typedef ToUniquePtr UniquePtrType; #else // #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) typedef Grid* UniquePtrType; From d0cbcf0a45ed274f21a310a923adc9122c1dd1e7 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Wed, 29 May 2019 10:35:36 +0200 Subject: [PATCH 10/22] Initial column iterator utility for verteq. --- opm/grid/common/VerteqColumnUtility.hpp | 56 ++++++++++ opm/grid/polyhedralgrid/grid.hh | 34 +++++- .../polyhedralgrid/verteqcolumnutility.hh | 102 ++++++++++++++++++ tests/grid_test.cc | 33 +++++- 4 files changed, 223 insertions(+), 2 deletions(-) create mode 100644 opm/grid/common/VerteqColumnUtility.hpp create mode 100644 opm/grid/polyhedralgrid/verteqcolumnutility.hh diff --git a/opm/grid/common/VerteqColumnUtility.hpp b/opm/grid/common/VerteqColumnUtility.hpp new file mode 100644 index 000000000..c0472ff79 --- /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/grid.hh b/opm/grid/polyhedralgrid/grid.hh index f55e3af4c..9c712d7da 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -38,6 +38,10 @@ #include #include +#include + +#include + namespace Dune { @@ -162,6 +166,7 @@ namespace Dune public: typedef UnstructuredGrid UnstructuredGridType; + typedef Opm::TopSurf TopSurfaceGridType; protected: struct UnstructuredGridDeleter @@ -326,6 +331,7 @@ namespace Dune const std::vector& poreVolumes = std::vector ()) : gridPtr_( createGrid( deck, poreVolumes ) ), grid_( *gridPtr_ ), + topSurfaceGrid_( nullptr ), comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), @@ -344,6 +350,7 @@ namespace Dune const std::vector< double >& dx ) : gridPtr_( createGrid( n, dx ) ), grid_( *gridPtr_ ), + topSurfaceGrid_( nullptr ), comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), @@ -361,6 +368,7 @@ namespace Dune explicit PolyhedralGrid ( UnstructuredGridPtr &&gridPtr ) : gridPtr_( std::move( gridPtr ) ), grid_( *gridPtr_ ), + topSurfaceGrid_( nullptr ), //polyhedralMesh_( grid_ ), comm_( *this ), leafIndexSet_( *this ), @@ -380,6 +388,26 @@ namespace Dune explicit PolyhedralGrid ( const UnstructuredGridType& grid ) : gridPtr_(), grid_( grid ), + topSurfaceGrid_( nullptr ), + comm_( *this ), + leafIndexSet_( *this ), + globalIdSet_( *this ), + localIdSet_( *this ) + { + 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 ), @@ -396,6 +424,8 @@ namespace Dune /** \} */ + const TopSurfaceGridType* topSurfaceGrid() const { return topSurfaceGrid_; } + /** \name Size Methods * \{ */ @@ -1576,11 +1606,12 @@ namespace Dune } - protected: UnstructuredGridPtr gridPtr_; const UnstructuredGridType& grid_; + const TopSurfaceGridType* topSurfaceGrid_; + PolyhedralMeshType polyhedralMesh_; CollectiveCommunication comm_; @@ -1695,5 +1726,6 @@ namespace Dune #include #include #include +#include #endif // #ifndef DUNE_POLYHEDRALGRID_GRID_HH diff --git a/opm/grid/polyhedralgrid/verteqcolumnutility.hh b/opm/grid/polyhedralgrid/verteqcolumnutility.hh new file mode 100644 index 000000000..1b2ddfd4c --- /dev/null +++ b/opm/grid/polyhedralgrid/verteqcolumnutility.hh @@ -0,0 +1,102 @@ +#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_cells[ entity.impl().index() ]; + return IteratorType( topSurfaceGrid_, + topSurfaceGrid_->col_cellpos[ colIdx ] ); + } + else + return end( entity ); + } + + IteratorType end( const Entity& entity ) const + { + if( topSurfaceGrid_ ) + { + const int colIdx = topSurfaceGrid_->col_cells[ entity.impl().index() ]; + return IteratorType( topSurfaceGrid_, + topSurfaceGrid_->col_cellpos[ colIdx+1 ] ); + } + else + return IteratorType( nullptr, -1 ); + } + }; + +} // end namespace Dune +#endif diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 972f69256..f3a421f52 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -14,6 +14,8 @@ #include +#include + #define DISABLE_DEPRECATED_METHOD_CHECK 1 #if DUNE_VERSION_NEWER(DUNE_GRID,2,5) #include @@ -54,6 +56,7 @@ void testGridIteration( const GridView& gridView ) int numElem = 0; ElemIterator elemIt = gridView.template begin<0>(); ElemIterator elemEndIt = gridView.template end<0>(); + /* for (; elemIt != elemEndIt; ++elemIt) { const Geometry& elemGeom = elemIt->geometry(); if (std::abs(elemGeom.volume() - 1.0) > 1e-8) @@ -95,16 +98,43 @@ void testGridIteration( const GridView& gridView ) } } - if (numIs != 6) + if (numIs != 2 * GridView::dimension ) std::cout << "number of intersections is wrong for element " << numElem << "\n"; ++ numElem; } + */ if (numElem != 2*2*2) std::cout << "number of elements is wrong: " << numElem << "\n"; } +template +void testVerteq(const Grid& grid) +{ + typedef typename Grid::LeafGridView GridView; + GridView gv = grid.leafGridView(); + + Dune::VerteqColumnUtility< Grid > verteqUtil ( grid ); + + std::cout << "Start checking vertical equilibirum utility." << std::endl; + + const auto end = gv.template end< 0 > (); + for( auto it = gv.template begin< 0 > (); it != end; ++it ) + { + const auto& entity = *it ; + const auto endCol = verteqUtil.end( entity ); + for( auto col = verteqUtil.begin( entity ); col != endCol; ++col ) + { + const auto& collCell = *col; + std::cout << "Column cell [ " << collCell.index() + << " ]: h = " << collCell.h() + << " dz = " << collCell.dz() << std::endl; + } + } +} + + template void testGrid(Grid& grid, const std::string& name) { @@ -150,6 +180,7 @@ void testGrid(Grid& grid, const std::string& name) vtkWriter.write(name, Dune::VTK::ascii); } + testVerteq( grid ); } int main(int argc, char** argv ) From d7c917a3578d588932863faf27d1757132ed0a94 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Wed, 29 May 2019 12:10:29 +0200 Subject: [PATCH 11/22] [bugfix][VerteqUtil] Fix some issues with the vertical column access. --- opm/grid/polyhedralgrid/geometry.hh | 4 ++-- opm/grid/polyhedralgrid/grid.hh | 1 + opm/grid/polyhedralgrid/verteqcolumnutility.hh | 15 +++++++++------ tests/grid_test.cc | 9 +++++---- 4 files changed, 17 insertions(+), 12 deletions(-) diff --git a/opm/grid/polyhedralgrid/geometry.hh b/opm/grid/polyhedralgrid/geometry.hh index 68252aa6e..a4b7d951f 100644 --- a/opm/grid/polyhedralgrid/geometry.hh +++ b/opm/grid/polyhedralgrid/geometry.hh @@ -77,7 +77,7 @@ namespace Dune GlobalCoordinate global(const LocalCoordinate& local) const { const GeometryType geomType = type(); - if( geomType.isCube() ) + if( geomType.isHexahedron() ) { assert( mydimension == 3 ); assert( coorddimension == 3 ); @@ -118,7 +118,7 @@ namespace Dune } return xyz; } - else if ( geomType.isNone() ) + else if ( geomType.isNone() || geomType.isCube() ) { // if no geometry type return the center of the element return center(); diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 9c712d7da..b72305932 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -413,6 +413,7 @@ namespace Dune globalIdSet_( *this ), localIdSet_( *this ) { + std::cout << "Creating TopSurfaceGrid" << topSurfaceGrid_ << std::endl; init(); } diff --git a/opm/grid/polyhedralgrid/verteqcolumnutility.hh b/opm/grid/polyhedralgrid/verteqcolumnutility.hh index 1b2ddfd4c..f3ea18dda 100644 --- a/opm/grid/polyhedralgrid/verteqcolumnutility.hh +++ b/opm/grid/polyhedralgrid/verteqcolumnutility.hh @@ -77,9 +77,9 @@ namespace Dune { if( topSurfaceGrid_ ) { - const int colIdx = topSurfaceGrid_->col_cells[ entity.impl().index() ]; - return IteratorType( topSurfaceGrid_, - topSurfaceGrid_->col_cellpos[ colIdx ] ); + const int colIdx = topSurfaceGrid_->col_cellpos[ entity.impl().index() ]; + std::cout << "begin " << colIdx << std::endl; + return IteratorType( topSurfaceGrid_, colIdx ); } else return end( entity ); @@ -89,12 +89,15 @@ namespace Dune { if( topSurfaceGrid_ ) { - const int colIdx = topSurfaceGrid_->col_cells[ entity.impl().index() ]; - return IteratorType( topSurfaceGrid_, - topSurfaceGrid_->col_cellpos[ colIdx+1 ] ); + 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 ); + } } }; diff --git a/tests/grid_test.cc b/tests/grid_test.cc index f3a421f52..1a55d592e 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -56,7 +56,6 @@ void testGridIteration( const GridView& gridView ) int numElem = 0; ElemIterator elemIt = gridView.template begin<0>(); ElemIterator elemEndIt = gridView.template end<0>(); - /* for (; elemIt != elemEndIt; ++elemIt) { const Geometry& elemGeom = elemIt->geometry(); if (std::abs(elemGeom.volume() - 1.0) > 1e-8) @@ -103,7 +102,6 @@ void testGridIteration( const GridView& gridView ) ++ numElem; } - */ if (numElem != 2*2*2) std::cout << "number of elements is wrong: " << numElem << "\n"; @@ -123,6 +121,7 @@ void testVerteq(const Grid& grid) for( auto it = gv.template begin< 0 > (); it != end; ++it ) { const auto& entity = *it ; + std::cout << "Start column for entity " << entity.impl().index() << std::endl; const auto endCol = verteqUtil.end( entity ); for( auto col = verteqUtil.begin( entity ); col != endCol; ++col ) { @@ -132,12 +131,15 @@ void testVerteq(const Grid& grid) << " dz = " << collCell.dz() << std::endl; } } + std::cout << "Finished checking vertical equilibirum utility." << std::endl; } template void testGrid(Grid& grid, const std::string& name) { + testVerteq( grid ); + typedef typename Grid::LeafGridView GridView; #if DUNE_VERSION_NEWER(DUNE_GRID,2,5) try { @@ -180,7 +182,6 @@ void testGrid(Grid& grid, const std::string& name) vtkWriter.write(name, Dune::VTK::ascii); } - testVerteq( grid ); } int main(int argc, char** argv ) @@ -223,11 +224,11 @@ int main(int argc, char** argv ) Grid2D tsDune (*ts); std::cout << "tsDune after " << std::endl; testGrid ( tsDune, "ts"); + return 0; #endif Dune::GridPtr< Grid > gridPtr( dgfFile ); testGrid( *gridPtr, "polyhedralgrid-dgf" ); - } /* From f19861360ba8b990bcfc78c5ab8e18acd659962d Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 13 Sep 2019 14:55:50 +0200 Subject: [PATCH 12/22] Squash --- opm/grid/polyhedralgrid/geometry.hh | 229 ++++++++++++++++------------ opm/grid/polyhedralgrid/grid.hh | 5 +- opm/grid/polyhedralgrid/gridview.hh | 1 - tests/grid_test.cc | 57 +++++-- 4 files changed, 178 insertions(+), 114 deletions(-) diff --git a/opm/grid/polyhedralgrid/geometry.hh b/opm/grid/polyhedralgrid/geometry.hh index a4b7d951f..b3f30926a 100644 --- a/opm/grid/polyhedralgrid/geometry.hh +++ b/opm/grid/polyhedralgrid/geometry.hh @@ -9,6 +9,7 @@ #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 5 ) #include +#include #else #include #include @@ -41,6 +42,68 @@ namespace Dune typedef Dune::FieldVector< ctype, coorddimension > GlobalCoordinate; typedef Dune::FieldVector< ctype, mydimension > LocalCoordinate; + typedef typename Grid::Traits::ExtraData ExtraData; + typedef typename Grid::Traits::template Codim::EntitySeed EntitySeed; + + template< class ct > + struct PolyhedralMultiLinearGeometryTraits + : public Dune::MultiLinearGeometryTraits< ct > + { + struct Storage + { + struct Iterator + : public Dune::ForwardIteratorFacade< Iterator, GlobalCoordinate, GlobalCoordinate > + { + const Storage* data_; + int count_; + explicit Iterator( const Storage* ptr, int count ) : data_( ptr ), count_( count ) {} + + GlobalCoordinate dereference() const { return data_->corner( count_ ); } + void increment() { ++count_; } + + bool equals( const Iterator& other ) const { return count_ == other.count_; } + }; + + ExtraData data_; + // host geometry object + EntitySeed seed_; + + Storage( ExtraData data, EntitySeed seed ) + : data_( data ), seed_( seed ) + {} + + Storage( ExtraData data ) + : data_( data ), seed_() + {} + + ExtraData data() const { return data_; } + bool isValid () const { return seed_.isValid(); } + + GlobalCoordinate operator [] (const int i) const { return corner( i ); } + + Iterator begin() const { return Iterator(this, 0); } + Iterator end () const { return Iterator(this, corners()); } + + int corners () const { return data()->corners( seed_ ); } + GlobalCoordinate corner ( const int i ) const { return data()->corner( seed_, i ); } + GlobalCoordinate center () const { return data()->centroids( seed_ ); } + + ctype volume() const { return data()->volumes( seed_ ); } + }; + + template + struct CornerStorage + { + typedef Storage Type; + }; + }; + + typedef Dune::MultiLinearGeometry< ctype, mydimension, coorddimension, PolyhedralMultiLinearGeometryTraits > + MultiLinearGeometryType; + + typedef typename PolyhedralMultiLinearGeometryTraits< ctype > ::template + CornerStorage::Type CornerStorageType; + //! type of jacobian inverse transposed typedef FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed; @@ -54,137 +117,105 @@ namespace Dune typedef Dune::GenericGeometry::MatrixHelper< Dune::GenericGeometry::DuneCoordTraits< ctype > > MatrixHelperType; #endif - typedef typename Grid::Traits::ExtraData ExtraData; - typedef typename Grid::Traits::template Codim::EntitySeed EntitySeed; - explicit PolyhedralGridBasicGeometry ( ExtraData data ) - : data_( data ), - seed_( ) + : storage_( data ) {} PolyhedralGridBasicGeometry ( ExtraData data, const EntitySeed& seed ) - : data_( data ), - seed_( seed ) - {} + : storage_( data, seed ) + { + GeometryType myType = type(); + if( ! myType.isNone() && storage_.isValid() ) + { + if( myType.isLine() && storage_.corners() != 2 ) + { + std::cout << myType << " " << storage_.corners() << std::endl; + std::abort(); + } + + //std::cout << myType << " " << storage_.corners() << std::endl; + geometryImpl_.reset( new MultiLinearGeometryType(myType, storage_) ); + } + } GeometryType type () const { return data()->geomTypes(codimension)[0]; } - bool affine () const { return false; } + bool affine () const { return (geometryImpl_) ? geometryImpl_->affine() : false; } - int corners () const { return data()->corners( seed_ ); } - GlobalCoordinate corner ( const int i ) const { return data()->corner( seed_, i ); } - GlobalCoordinate center () const { return data()->centroids( seed_ ); } + int corners () const { return storage_.corners(); } + GlobalCoordinate corner ( const int i ) const { return storage_.corner( i ); } + GlobalCoordinate center () const { return storage_.center(); } GlobalCoordinate global(const LocalCoordinate& local) const { - const GeometryType geomType = type(); - if( geomType.isHexahedron() ) - { - assert( mydimension == 3 ); - assert( coorddimension == 3 ); - - // uvw = { (1-u, 1-v, 1-w), (u, v, w) } - LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local }; - uvw[0] -= local; - - //const ReferenceElement< ctype , mydimension > & refElement = - // ReferenceElements< ctype, mydimension >::general( type() ); - - // Access pattern for uvw matching ordering of corners. - const int pat[8][3] = { { 0, 0, 0 }, - { 1, 0, 0 }, - { 0, 1, 0 }, - { 1, 1, 0 }, - { 0, 0, 1 }, - { 1, 0, 1 }, - { 0, 1, 1 }, - { 1, 1, 1 } }; - - const int nCorners = corners(); - //refElement.size( mydimension ); - - GlobalCoordinate xyz(0.0); - for (int i = 0; i < nCorners ; ++i) - { - GlobalCoordinate cornerContrib = corner(i); - //LocalCoordinate refCorner = refElement.position(i,mydimension); - double factor = 1.0; - for (int j = 0; j < mydimension; ++j) - { - //factor *= uvw[ refCorner[ j ] ][ j ]; - factor *= uvw[ pat[ i ][ j ] ][ j ]; - } - cornerContrib *= factor; - xyz += cornerContrib; - } - return xyz; - } - else if ( geomType.isNone() || geomType.isCube() ) + if( geometryImpl_ ) { - // if no geometry type return the center of the element - return center(); - } - else - { - DUNE_THROW(NotImplemented,"global for geometry type " << geomType << " is not supported!"); - return GlobalCoordinate( 0 ); + return geometryImpl_->global( local ); } + + return center(); } /// Mapping from the cell to the reference domain. /// May be slow. - LocalCoordinate local(const GlobalCoordinate& y) const + LocalCoordinate local(const GlobalCoordinate& global) const { - const GeometryType geomType = type(); - if( geomType.isCube() ) + if( geometryImpl_ ) { - // This code is modified from dune/grid/genericgeometry/mapping.hh - // \todo: Implement direct computation. - const ctype epsilon = 1e-12; - const ReferenceElement< ctype , mydimension > & refElement = - ReferenceElements< ctype, mydimension >::general(type()); - - LocalCoordinate x = refElement.position(0,0); - LocalCoordinate dx; - do { - // DF^n dx^n = F^n, x^{n+1} -= dx^n - JacobianTransposed JT = jacobianTransposed(x); - GlobalCoordinate z = global(x); - z -= y; - MatrixHelperType::template xTRightInvA(JT, z, dx ); - x -= dx; - } while (dx.two_norm2() > epsilon*epsilon); - return x; + return geometryImpl_->local( global ); } - else if ( geomType.isNone() ) + + // if no geometry type return a vector filled with 1 + return LocalCoordinate( 1 ); + } + + ctype integrationElement ( const LocalCoordinate &local ) const + { + if( geometryImpl_ ) { - // if no geometry type return a vector filled with 1 - return LocalCoordinate( 1 ); + return geometryImpl_->integrationElement( local ); } - else + return volume(); + } + + ctype volume () const + { + if( geometryImpl_ ) { - DUNE_THROW(NotImplemented,"local for geometry type " << geomType << " is not supported!"); - return LocalCoordinate( 0 ); + return geometryImpl_->volume(); } + return storage_.volume(); } - ctype integrationElement ( const LocalCoordinate & ) const { return volume(); } - ctype volume () const { return data()->volumes( seed_ ); } - #if DUNE_VERSION_NEWER(DUNE_GRID,2,4) - JacobianTransposed jacobianTransposed ( const LocalCoordinate & ) const + JacobianTransposed jacobianTransposed ( const LocalCoordinate & local ) const { + if( geometryImpl_ ) + { + return geometryImpl_->jacobianTransposed( local ); + } + DUNE_THROW(NotImplemented,"jacobianTransposed not implemented"); return JacobianTransposed( 0 ); } - JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate & ) const + JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate & local ) const { + if( geometryImpl_ ) + { + return geometryImpl_->jacobianInverseTransposed( local ); + } + DUNE_THROW(NotImplemented,"jacobianInverseTransposed not implemented"); return JacobianInverseTransposed( 0 ); } #else const JacobianTransposed& jacobianTransposed ( const LocalCoordinate &local ) const { + if( geometryImpl_ ) + { + return geometryImpl_->jacobianTransposed( local ); + } + DUNE_THROW(NotImplemented,"jacobianTransposed not implemented"); static const JacobianTransposed jac( 0 ); return jac; @@ -192,18 +223,22 @@ namespace Dune const JacobianInverseTransposed& jacobianInverseTransposed ( const LocalCoordinate &local ) const { + if( geometryImpl_ ) + { + return geometryImpl_->jacobianInverseTransposed( local ); + } + DUNE_THROW(NotImplemented,"jacobianInverseTransposed not implemented"); static const JacobianInverseTransposed jac( 0 ); return jac; } #endif - ExtraData data() const { return data_; } + ExtraData data() const { return storage_.data(); } protected: - ExtraData data_; - // host geometry object - EntitySeed seed_; + CornerStorageType storage_; + std::shared_ptr< MultiLinearGeometryType > geometryImpl_; }; diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index b72305932..282a8c7ab 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -985,7 +985,8 @@ namespace Dune } case 1: { - return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; + //return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; + return grid_.face_nodepos[ index+1 ] - grid_.face_nodepos[ index ]; } case dim: { @@ -1009,7 +1010,7 @@ namespace Dune } case 1: { - const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + i]; + const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + i ]; return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex ); } case dim: diff --git a/opm/grid/polyhedralgrid/gridview.hh b/opm/grid/polyhedralgrid/gridview.hh index 6057c57a6..7e926b31f 100644 --- a/opm/grid/polyhedralgrid/gridview.hh +++ b/opm/grid/polyhedralgrid/gridview.hh @@ -57,7 +57,6 @@ namespace Dune PolyhedralGridView ( const Grid &grid, const int level = 0 ) : grid_( &grid ) { - static_cast(&level); } const Grid &grid () const diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 1a55d592e..20a632745 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -44,7 +45,7 @@ const char *deckString = "4*100.0 /\n"; template -void testGridIteration( const GridView& gridView ) +void testGridIteration( const GridView& gridView, const int nElem ) { typedef typename GridView::template Codim<0>::Iterator ElemIterator; typedef typename GridView::IntersectionIterator IsIt; @@ -103,10 +104,11 @@ void testGridIteration( const GridView& gridView ) ++ numElem; } - if (numElem != 2*2*2) - std::cout << "number of elements is wrong: " << numElem << "\n"; + if (numElem != nElem ) + std::cout << "number of elements is wrong: " << numElem << ", expected " << nElem << std::endl; } +/* template void testVerteq(const Grid& grid) { @@ -133,12 +135,13 @@ void testVerteq(const Grid& grid) } std::cout << "Finished checking vertical equilibirum utility." << std::endl; } +*/ template -void testGrid(Grid& grid, const std::string& name) +void testGrid(Grid& grid, const std::string& name, const size_t nElem, const size_t nVertices) { - testVerteq( grid ); + //testVerteq( grid ); typedef typename Grid::LeafGridView GridView; #if DUNE_VERSION_NEWER(DUNE_GRID,2,5) @@ -152,15 +155,15 @@ void testGrid(Grid& grid, const std::string& name) #endif std::cout << name << std::endl; - testGridIteration( grid.leafGridView() ); + testGridIteration( grid.leafGridView(), nElem ); std::cout << "create vertex mapper\n"; Dune::MultipleCodimMultipleGeomTypeMapper mapper(grid.leafGridView()); std::cout << "VertexMapper.size(): " << mapper.size() << "\n"; - if (mapper.size() != 27) { - std::cout << "Wrong size of vertex mapper. Expected 27!\n"; + if (mapper.size() != nVertices ) { + std::cout << "Wrong size of vertex mapper. Expected " << nVertices << "!" << std::endl; //std::abort(); } @@ -201,17 +204,20 @@ int main(int argc, char** argv ) //dgfFile << "#" << std::endl; #if HAVE_ECL_INPUT + /* Opm::Parser parser; const auto deck = parser.parseString(deckString); std::vector porv; + */ #endif // test PolyhedralGrid { typedef Dune::PolyhedralGrid< 3, 3 > Grid; #if HAVE_ECL_INPUT + /* Grid grid(deck, porv); - testGrid( grid, "polyhedralgrid" ); + testGrid( grid, "polyhedralgrid", 8, 27 ); Opm::TopSurf* ts; ts = Opm::TopSurf::create (grid); std::cout << ts->dimensions << std::endl; @@ -223,12 +229,35 @@ int main(int argc, char** argv ) std::cout << "tsDune for " << std::endl; Grid2D tsDune (*ts); std::cout << "tsDune after " << std::endl; - testGrid ( tsDune, "ts"); - return 0; - + testGrid ( tsDune, "ts", 27 ); + */ #endif - Dune::GridPtr< Grid > gridPtr( dgfFile ); - testGrid( *gridPtr, "polyhedralgrid-dgf" ); + /* + { + std::cout <<"Check 3d grid" << std::endl; + Dune::GridPtr< Grid > gridPtr( dgfFile ); + testGrid( *gridPtr, "polyhedralgrid-dgf", std::pow(3, Grid::dimension) ); + } + */ + + { + 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 << "2 2" << 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 ); + } } /* From ccf4cd196045b4b825b4eb7cb43fad48112f9a6f Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 13 Sep 2019 16:39:59 +0200 Subject: [PATCH 13/22] added grid_test. --- CMakeLists_files.cmake | 3 ++- opm/grid/polyhedralgrid/grid.hh | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index ef448ca69..e99b52b05 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -100,7 +100,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_geom2d.cpp tests/test_gridutilities.cpp tests/test_minpvprocessor.cpp - tests/grid_test.cc + #tests/grid_test.cc tests/p2pcommunicator_test.cc tests/test_repairzcorn.cpp tests/test_sparsetable.cpp @@ -132,6 +132,7 @@ list (APPEND TEST_DATA_FILES list (APPEND EXAMPLE_SOURCE_FILES examples/finitevolume/finitevolume.cc examples/mirror_grid.cpp + tests/grid_test.cc ) # programs listed here will not only be compiled, but also marked for diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 282a8c7ab..654c8cde5 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1050,6 +1050,7 @@ namespace Dune case 0: return 1; case 1: + //std::cout << "Subent codim 1 " << grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ] << std::endl; return grid_.cell_facepos[ index+1 ] - grid_.cell_facepos[ index ]; case dim: return cellVertices_[ index ].size(); From a62b05caf427f16721fffd7a229a8ed06a7ced65 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Wed, 2 Oct 2019 11:22:13 +0200 Subject: [PATCH 14/22] Still not working. --- opm/grid/polyhedralgrid/grid.hh | 6 ++++++ opm/grid/polyhedralgrid/indexset.hh | 3 +++ 2 files changed, 9 insertions(+) diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 654c8cde5..1c48c8c47 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1322,6 +1322,7 @@ namespace Dune // sort vertices such that they comply with the dune cube reference element if( grid_.cell_facetag ) { + std::cout << "Sort vertices " << std::endl; typedef std::array KeyType; std::map< const KeyType, const int > vertexFaceTags; const int vertexFacePattern [8][3] = { @@ -1411,6 +1412,11 @@ namespace Dune cellVertices_[ c ][ (*vx).second ] = (*it).first ; } } + + for( int c=0; c ( i ) ); else if ( codim == dimension ) + { + std::cout << "Return vertex number " << std::endl; return index( grid().getRealImplementation( entity ).template subEntity< dimension > ( i ) ); + } else { DUNE_THROW(NotImplemented,"codimension not available"); From 74a5f6589483368d502603e61e9d468e11ce3ab0 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Sat, 26 Oct 2019 22:01:57 +0200 Subject: [PATCH 15/22] Make compile with master. --- opm/grid/polyhedralgrid/grid.hh | 14 +++++--------- tests/grid_test.cc | 2 +- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 1c48c8c47..8b0411125 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -890,19 +890,15 @@ namespace Dune const auto eclipseGrid = std::make_shared(deck, rawactnum); struct grdecl g; - std::vector actnum; - std::vector coord; - std::vector zcorn; - std::vector mapaxes; g.dims[0] = eclipseGrid->getNX(); g.dims[1] = eclipseGrid->getNY(); g.dims[2] = eclipseGrid->getNZ(); - eclipseGrid->exportMAPAXES( mapaxes ); - eclipseGrid->exportCOORD( coord ); - eclipseGrid->exportZCORN( zcorn ); - eclipseGrid->exportACTNUM( actnum ); + std::vector mapaxes = eclipseGrid->getMAPAXES( ); + std::vector coord = eclipseGrid->getCOORD( ); + std::vector zcorn = eclipseGrid->getZCORN( ); + std::vector actnum = eclipseGrid->getACTNUM( ); g.coord = coord.data(); g.zcorn = zcorn.data(); @@ -920,7 +916,7 @@ namespace Dune const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2]; std::vector thickness(cartGridSize); for (size_t i = 0; i < cartGridSize; ++i) { - thickness[i] = eclipseGrid->getCellThicknes(i); + thickness[i] = eclipseGrid->getCellThickness(i); } const double z_tolerance = eclipseGrid->isPinchActive() ? eclipseGrid->getPinchThresholdThickness() : 0.0; mp.process(thickness, z_tolerance, poreVolumes, minpvv, actnum, opmfil, zcorn.data()); diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 20a632745..657a72fb9 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -198,7 +198,7 @@ int main(int argc, char** argv ) dgfFile << "Interval" << std::endl; dgfFile << "0 0 0" << std::endl; dgfFile << "2 2 2" << std::endl; - dgfFile << "2 2 2" << std::endl; + dgfFile << "1 1 2" << std::endl; dgfFile << "#" << std::endl; //dgfFile << "Simplex" << std::endl; //dgfFile << "#" << std::endl; From 43b075bd4401b31110d858f193d9550dc2fd1ed0 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Sat, 26 Oct 2019 22:36:29 +0200 Subject: [PATCH 16/22] Bugfix in IdSet, make ids globally unique by adding the size of previous codim sets. --- opm/grid/polyhedralgrid/idset.hh | 22 +++++++++++++++++----- opm/grid/polyhedralgrid/indexset.hh | 1 - tests/grid_test.cc | 8 ++++---- 3 files changed, 21 insertions(+), 10 deletions(-) diff --git a/opm/grid/polyhedralgrid/idset.hh b/opm/grid/polyhedralgrid/idset.hh index 47928e949..9bbf2a2b9 100644 --- a/opm/grid/polyhedralgrid/idset.hh +++ b/opm/grid/polyhedralgrid/idset.hh @@ -13,12 +13,13 @@ namespace Dune template< int dim, int dimworld, typename coord_t > class PolyhedralGridIdSet - : public IdSet< PolyhedralGrid< dim, dimworld, coord_t >, PolyhedralGridIdSet< dim, dimworld, coord_t >, /*IdType=*/int > + : public IdSet< PolyhedralGrid< dim, dimworld, coord_t >, PolyhedralGridIdSet< dim, dimworld, coord_t >, size_t /*IdType=int*/ > { public: typedef PolyhedralGrid< dim, dimworld, coord_t > Grid; typedef typename std::remove_const< Grid >::type::Traits Traits; - typedef typename Traits::Index IdType; + //typedef typename Traits::Index IdType; + typedef size_t IdType; typedef PolyhedralGridIdSet< dim, dimworld, coord_t > This; typedef IdSet< Grid, This, IdType > Base; @@ -26,7 +27,13 @@ namespace Dune PolyhedralGridIdSet (const Grid& grid) : grid_( grid ), globalCellPtr_( grid_.globalCellPtr() ) - {} + { + codimOffset_[ 0 ] = 0; + for( int i=1; i<=dim; ++i ) + { + codimOffset_[ i ] = codimOffset_[ i-1 ] + grid.size( i-1 ); + } + } //! id meethod for entity and specific codim template< int codim > @@ -35,9 +42,11 @@ namespace Dune const int index = entity.seed().index(); // in case if (codim == 0 && globalCellPtr_ ) - return globalCellPtr_[ index ]; + return size_t( globalCellPtr_[ index ] ); else - return index; + { + return codimOffset_[ codim ] + index; + } } #if ! DUNE_VERSION_NEWER(DUNE_GRID,2,4) @@ -72,7 +81,9 @@ namespace Dune else if ( codim == 1 ) return id( entity.template subEntity< 1 >( i ) ); else if ( codim == dim ) + { return id( entity.template subEntity< dim >( i ) ); + } else { DUNE_THROW(NotImplemented,"codimension not available"); @@ -83,6 +94,7 @@ namespace Dune protected: const Grid& grid_; const int* globalCellPtr_; + size_t codimOffset_[ dim+1 ]; }; } // namespace Dune diff --git a/opm/grid/polyhedralgrid/indexset.hh b/opm/grid/polyhedralgrid/indexset.hh index 65e45043f..7240c995f 100644 --- a/opm/grid/polyhedralgrid/indexset.hh +++ b/opm/grid/polyhedralgrid/indexset.hh @@ -75,7 +75,6 @@ namespace Dune return index( grid().getRealImplementation( entity ).template subEntity< 1 > ( i ) ); else if ( codim == dimension ) { - std::cout << "Return vertex number " << std::endl; return index( grid().getRealImplementation( entity ).template subEntity< dimension > ( i ) ); } else diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 657a72fb9..a8686605e 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -198,7 +198,7 @@ int main(int argc, char** argv ) dgfFile << "Interval" << std::endl; dgfFile << "0 0 0" << std::endl; dgfFile << "2 2 2" << std::endl; - dgfFile << "1 1 2" << std::endl; + dgfFile << "1 1 1" << std::endl; dgfFile << "#" << std::endl; //dgfFile << "Simplex" << std::endl; //dgfFile << "#" << std::endl; @@ -247,15 +247,15 @@ int main(int argc, char** argv ) dgfFile << "DGF" << std::endl; dgfFile << "Interval" << std::endl; dgfFile << "0 0" << std::endl; - dgfFile << "2 2" << std::endl; - dgfFile << "2 2" << std::endl; + dgfFile << "1 1" << std::endl; + dgfFile << "1 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); + //testGrid( *gridPtr, "polyhedralgrid-dgf", 4, nVx); gridcheck( *gridPtr ); } } From 4dcfb59241b3ccdd8e87e1f99be73804c0f56fac Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 7 Nov 2019 09:43:52 +0100 Subject: [PATCH 17/22] Make grid_test work with DUNEs grid checks. --- opm/grid/UnstructuredGrid.c | 52 +++++++++++++++++++++++++++++++++ opm/grid/UnstructuredGrid.h | 2 ++ opm/grid/cart_grid.c | 2 +- opm/grid/polyhedralgrid/grid.hh | 1 + tests/grid_test.cc | 2 +- 5 files changed, 57 insertions(+), 2 deletions(-) diff --git a/opm/grid/UnstructuredGrid.c b/opm/grid/UnstructuredGrid.c index 271376795..3b258f6ca 100644 --- a/opm/grid/UnstructuredGrid.c +++ b/opm/grid/UnstructuredGrid.c @@ -166,6 +166,58 @@ 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(")\n"); + } + + 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 c2f8a543a..48cdf5ede 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/cart_grid.c b/opm/grid/cart_grid.c index 8226d1840..10a0d5ed7 100644 --- a/opm/grid/cart_grid.c +++ b/opm/grid/cart_grid.c @@ -609,8 +609,8 @@ fill_cart_topology_2d(struct UnstructuredGrid *G) /* Faces with y-normal */ for (j=0; j Grid; Dune::GridPtr< Grid > gridPtr( dgfFile ); From 0e8d8ec66bc5deb4c97ab4483aa65e2b560933ee Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 7 Nov 2019 14:41:49 +0100 Subject: [PATCH 18/22] Bugfix hasBoundaryIntersections and boundarySegmentIndex. --- opm/grid/UnstructuredGrid.c | 28 ++++++++- opm/grid/polyhedralgrid/entity.hh | 12 +--- opm/grid/polyhedralgrid/entityseed.hh | 4 +- opm/grid/polyhedralgrid/grid.hh | 82 ++++++++++++++++++++++--- opm/grid/polyhedralgrid/intersection.hh | 3 +- tests/grid_test.cc | 2 +- 6 files changed, 105 insertions(+), 26 deletions(-) diff --git a/opm/grid/UnstructuredGrid.c b/opm/grid/UnstructuredGrid.c index 3b258f6ca..e67bbf4ba 100644 --- a/opm/grid/UnstructuredGrid.c +++ b/opm/grid/UnstructuredGrid.c @@ -189,7 +189,11 @@ void print_grid(const struct UnstructuredGrid *grd ) { printf(" %d ",grd->face_nodes[ vx ]); } - printf(")\n"); + printf(")"); + int bnd = grd->face_cells[ 2*f ] < 0 || grd->face_cells[ 2*f+1 ] < 0; + if( bnd ) + printf(" boundary"); + printf("\n"); } for( int c=0; cnumber_of_cells; ++c ) @@ -217,6 +221,28 @@ void print_grid(const struct UnstructuredGrid *grd ) } printf(")\n"); } + + printf("Neighbors \n"); + for( int c=0; cnumber_of_cells; ++c ) + { + printf("Cell %d \n ", 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 ) + { + int face = grd->cell_faces[ f ]; + int bnd = grd->face_cells[ 2*face ] < 0 || grd->face_cells[ 2*face+1 ] < 0; + printf(" Face %d is boundary %d",f-grd->cell_facepos[ c ],bnd); + + /* + 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 diff --git a/opm/grid/polyhedralgrid/entity.hh b/opm/grid/polyhedralgrid/entity.hh index f372dd01c..fcfc02235 100644 --- a/opm/grid/polyhedralgrid/entity.hh +++ b/opm/grid/polyhedralgrid/entity.hh @@ -296,17 +296,7 @@ namespace Dune bool hasBoundaryIntersections () const { - const int faces = subEntities( 1 ); - for( int i=0; itemplate subEntity< 1 >( i ).seed().isValid() ) -#else - if( !this->template subEntity< 1 >( i )->seed().isValid() ) -#endif - return true; - } - return false; + return data()->hasBoundaryIntersections( this->seed() ); } #if ! DUNE_VERSION_NEWER(DUNE_GRID,2,4) diff --git a/opm/grid/polyhedralgrid/entityseed.hh b/opm/grid/polyhedralgrid/entityseed.hh index 1a2ad5398..3466aeef5 100644 --- a/opm/grid/polyhedralgrid/entityseed.hh +++ b/opm/grid/polyhedralgrid/entityseed.hh @@ -41,7 +41,9 @@ namespace Dune int index () const { return index_ ; } - bool isValid() const { return index_ != defaultIndex; } + // check that index is valid, which means >= 0 + // boundary faces can be arbitrary number < 0 + bool isValid() const { return index_ > defaultIndex; } bool equals(const PolyhedralGridEntitySeed& other) const { return index_ == other.index_; } diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 8ac6e0e80..eb1320a47 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -335,7 +335,8 @@ namespace Dune comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), - localIdSet_( *this ) + localIdSet_( *this ), + nBndSegments_( 0 ) { init(); } @@ -354,7 +355,8 @@ namespace Dune comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), - localIdSet_( *this ) + localIdSet_( *this ), + nBndSegments_( 0 ) { init(); } @@ -373,7 +375,8 @@ namespace Dune comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), - localIdSet_( *this ) + localIdSet_( *this ), + nBndSegments_( 0 ) { init(); } @@ -392,7 +395,8 @@ namespace Dune comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), - localIdSet_( *this ) + localIdSet_( *this ), + nBndSegments_( 0 ) { init(); } @@ -411,7 +415,8 @@ namespace Dune comm_( *this ), leafIndexSet_( *this ), globalIdSet_( *this ), - localIdSet_( *this ) + localIdSet_( *this ), + nBndSegments_( 0 ) { std::cout << "Creating TopSurfaceGrid" << topSurfaceGrid_ << std::endl; init(); @@ -513,7 +518,7 @@ namespace Dune */ size_t numBoundarySegments () const { - return 0; + return nBndSegments_; } /** \} */ @@ -1125,6 +1130,42 @@ namespace Dune } } + bool hasBoundaryIntersections(const typename Codim<0>::EntitySeed& seed ) const + { + const int faces = subEntities( seed, 1 ); + for( int f=0; ftemplate subEntitySeed<1>( seed, f ); + if( isBoundaryFace( faceSeed ) ) + return true; + } + return false; + } + + bool isBoundaryFace(const int face ) const + { + assert( face >= 0 && face < grid_.number_of_faces ); + const int facePos = 2 * face; + return ((grid_.face_cells[ facePos ] < 0) || (grid_.face_cells[ facePos+1 ] < 0)); + } + + bool isBoundaryFace(const typename Codim<1>::EntitySeed& faceSeed ) const + { + assert( faceSeed.valid() ); + return isBoundaryFace( faceSeed.index() ); + } + + int boundarySegmentIndex(const typename Codim<0>::EntitySeed& seed, const int face ) const + { + const auto faceSeed = this->template subEntitySeed<1>( seed, face ); + assert( faceSeed.isValid() ); + const int facePos = 2 * faceSeed.index(); + const int idx = std::min( grid_.face_cells[ facePos ], grid_.face_cells[ facePos+1 ]); + // check that this is actually the boundary + assert( idx < 0 ); + return -(idx+1); // +1 to include 0 boundary segment index + } + const std::vector< GeometryType > &geomTypes ( const unsigned int codim ) const { static std::vector< GeometryType > emptyDummy; @@ -1305,6 +1346,8 @@ namespace Dune protected: void init() { + std::cout << "PolyhedralGrid init" << std::endl; + // copy Cartesian dimensions for( int i=0; i<3; ++i ) { @@ -1581,14 +1624,31 @@ namespace Dune } } // end else of ( grid_.cell_facetag ) + nBndSegments_ = 0; unitOuterNormals_.resize( grid_.number_of_faces ); for( int face = 0; face < grid_.number_of_faces; ++face ) { - const int normalIdx = face * GlobalCoordinate :: dimension ; - GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx ); - normal /= normal.two_norm(); + const int normalIdx = face * GlobalCoordinate :: dimension ; + GlobalCoordinate normal = copyToGlobalCoordinate( grid_.face_normals + normalIdx ); + normal /= normal.two_norm(); + unitOuterNormals_[ face ] = normal; - unitOuterNormals_[ face ] = normal; + if( isBoundaryFace( face ) ) + { + // increase number if boundary segments + ++nBndSegments_; + const int facePos = 2 * face ; + // store negative number to indicate boundary + // the abstract value is the segment index + if( grid_.face_cells[ facePos ] < 0 ) + { + grid_.face_cells[ facePos ] = -nBndSegments_; + } + else if ( grid_.face_cells[ facePos+1 ] < 0 ) + { + grid_.face_cells[ facePos+1 ] = -nBndSegments_; + } + } } //print( std::cout, grid_ ); @@ -1634,6 +1694,8 @@ namespace Dune mutable GlobalIdSet globalIdSet_; mutable LocalIdSet localIdSet_; + size_t nBndSegments_; + private: // no copying PolyhedralGrid ( const PolyhedralGrid& ); diff --git a/opm/grid/polyhedralgrid/intersection.hh b/opm/grid/polyhedralgrid/intersection.hh index 5099a12d2..4d913487f 100644 --- a/opm/grid/polyhedralgrid/intersection.hh +++ b/opm/grid/polyhedralgrid/intersection.hh @@ -102,8 +102,7 @@ namespace Dune size_t boundarySegmentIndex () const { - // not implementable - return -1; + return data()->boundarySegmentIndex( seed_, intersectionIdx_); } LocalGeometry geometryInInside () const diff --git a/tests/grid_test.cc b/tests/grid_test.cc index e98f9936e..e6f27f1a5 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -248,7 +248,7 @@ int main(int argc, char** argv ) dgfFile << "Interval" << std::endl; dgfFile << "0 0" << std::endl; dgfFile << "1 1" << std::endl; - dgfFile << "4 4" << std::endl; + dgfFile << "2 2" << std::endl; dgfFile << "#" << std::endl; typedef Dune::PolyhedralGrid< 2, 2 > Grid; Dune::GridPtr< Grid > gridPtr( dgfFile ); From 8c2152df116985a2c1c70b44aa723d59f17cc012 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 7 Nov 2019 16:11:04 +0100 Subject: [PATCH 19/22] Fix 3d version. --- opm/grid/UnstructuredGrid.c | 85 +++++++++++++++++++++++---------- opm/grid/polyhedralgrid/grid.hh | 9 +++- tests/grid_test.cc | 8 ++-- 3 files changed, 71 insertions(+), 31 deletions(-) diff --git a/opm/grid/UnstructuredGrid.c b/opm/grid/UnstructuredGrid.c index e67bbf4ba..dced04f76 100644 --- a/opm/grid/UnstructuredGrid.c +++ b/opm/grid/UnstructuredGrid.c @@ -181,6 +181,58 @@ void print_grid(const struct UnstructuredGrid *grd ) printf(")\n"); } + static const int map[ 6 ][ 4 ] = { {0, 1, 3, 2}, // face 0 + {0, 1, 3, 2}, // face 1 + {0, 3, 1, 2}, // face 2 + {0, 3, 1, 2}, // face 3 + {0, 1, 3, 2}, // face 4 + {0, 1, 3, 2}}; // face 5 + + if( dim == 3 ) + { + for( int face=0; facenumber_of_faces; ++face ) + { + int c = grd->face_cells[ 2*face ]; + if( c < 0 ) + c = grd->face_cells[ 2*face+1 ]; + if( c < 0 ) exit(1); + + //for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) + const int faces = grd->cell_facepos[ c+1 ] - grd->cell_facepos[ c ]; + if( faces == 6 ) // Cartesian hexahedral case + { + printf("Swap vertices for face %d \n", face); + int i = -1; + for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) + { + if( grd->cell_faces[ f ] == face ) + { + i = grd->cell_facepos[ c+1 ] - f - 1; + break; + } + } + printf("Face is %d\n", i); + int vxPos = grd->face_nodepos[ face ]; + const int vxs[ 4 ] = { grd->face_nodes[ vxPos ], + grd->face_nodes[ vxPos+1 ], + grd->face_nodes[ vxPos+2 ], + grd->face_nodes[ vxPos+3 ] }; + for( int vx=0; vx<4; ++vx ) + { + grd->face_nodes[ vxPos+vx ] = vxs[ map[ i ][ vx ]]; + } + + /* + for( int vx=grd->face_nodepos[ f ]; vxface_nodepos[ f+1 ]; ++vx ) + { + printf(" %d ",grd->face_nodes[ vx ]); + } + printf(")\n"); + */ + } + } + } + for( int f=0; fnumber_of_faces; ++f ) { //for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) @@ -196,42 +248,23 @@ void print_grid(const struct UnstructuredGrid *grd ) printf("\n"); } - 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 ) + if( dim == 2 ) { - 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 ) + for( int c=0; cnumber_of_cells; ++c ) { - 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"); - */ + 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 ] ); } - printf(")\n"); } - printf("Neighbors \n"); for( int c=0; cnumber_of_cells; ++c ) { - printf("Cell %d \n ", 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 ) { - int face = grd->cell_faces[ f ]; - int bnd = grd->face_cells[ 2*face ] < 0 || grd->face_cells[ 2*face+1 ] < 0; - printf(" Face %d is boundary %d",f-grd->cell_facepos[ c ],bnd); + printf(" %d ",grd->cell_faces[ f ]); /* for( int vx=grd->face_nodepos[ f ]; vxface_nodepos[ f+1 ]; ++vx ) diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index eb1320a47..81aecc2a0 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -955,13 +955,13 @@ namespace Dune if( dim == 2 ) { cgrid = create_grid_cart2d( n[ 0 ], n[ 1 ], dx[ 0 ], dx[ 1 ] ); - print_grid( cgrid ); } else if ( dim == 3 ) { cgrid = create_grid_hexa3d( n[ 0 ], n[ 1 ], n[ 2 ], dx[ 0 ], dx[ 1 ], dx[ 2 ] ); } + print_grid( cgrid ); if (!cgrid) { OPM_THROW(std::runtime_error, "Failed to construct grid."); } @@ -1012,6 +1012,13 @@ namespace Dune } case 1: { + static const int map[ 6 ][ 4 ] = { {0, 1, 3, 2}, // face 0 + {0, 1, 3, 2}, // face 1 + {0, 3, 1, 2}, // face 2 + {0, 3, 1, 2}, // face 3 + {0, 1, 3, 2}, // face 4 + {0, 1, 3, 2}}; // face 5 + const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + i ]; return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex ); } diff --git a/tests/grid_test.cc b/tests/grid_test.cc index e6f27f1a5..216e40fe4 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -197,8 +197,8 @@ int main(int argc, char** argv ) dgfFile << "DGF" << std::endl; dgfFile << "Interval" << std::endl; dgfFile << "0 0 0" << std::endl; - dgfFile << "2 2 2" << std::endl; dgfFile << "1 1 1" << std::endl; + dgfFile << "4 4 4" << std::endl; dgfFile << "#" << std::endl; //dgfFile << "Simplex" << std::endl; //dgfFile << "#" << std::endl; @@ -232,13 +232,13 @@ int main(int argc, char** argv ) testGrid ( tsDune, "ts", 27 ); */ #endif - /* { std::cout <<"Check 3d grid" << std::endl; Dune::GridPtr< Grid > gridPtr( dgfFile ); - testGrid( *gridPtr, "polyhedralgrid-dgf", std::pow(3, Grid::dimension) ); + //testGrid( *gridPtr, "polyhedralgrid-dgf", std::pow(3, int(Grid::dimension)) ); + gridcheck( *gridPtr ); } - */ + return 0; { std::cout <<"Check 2d grid" << std::endl; From 8c614207234c30dffb9af0ce991694aa279994c2 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 7 Nov 2019 17:01:28 +0100 Subject: [PATCH 20/22] Fixed 2d/3d DUNE grid mapping issue. --- opm/grid/UnstructuredGrid.c | 54 ++------------------------------- opm/grid/cart_grid.c | 2 +- opm/grid/polyhedralgrid/grid.hh | 27 ++++++++++------- tests/grid_test.cc | 1 - 4 files changed, 19 insertions(+), 65 deletions(-) diff --git a/opm/grid/UnstructuredGrid.c b/opm/grid/UnstructuredGrid.c index dced04f76..6ede525c8 100644 --- a/opm/grid/UnstructuredGrid.c +++ b/opm/grid/UnstructuredGrid.c @@ -181,58 +181,6 @@ void print_grid(const struct UnstructuredGrid *grd ) printf(")\n"); } - static const int map[ 6 ][ 4 ] = { {0, 1, 3, 2}, // face 0 - {0, 1, 3, 2}, // face 1 - {0, 3, 1, 2}, // face 2 - {0, 3, 1, 2}, // face 3 - {0, 1, 3, 2}, // face 4 - {0, 1, 3, 2}}; // face 5 - - if( dim == 3 ) - { - for( int face=0; facenumber_of_faces; ++face ) - { - int c = grd->face_cells[ 2*face ]; - if( c < 0 ) - c = grd->face_cells[ 2*face+1 ]; - if( c < 0 ) exit(1); - - //for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) - const int faces = grd->cell_facepos[ c+1 ] - grd->cell_facepos[ c ]; - if( faces == 6 ) // Cartesian hexahedral case - { - printf("Swap vertices for face %d \n", face); - int i = -1; - for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) - { - if( grd->cell_faces[ f ] == face ) - { - i = grd->cell_facepos[ c+1 ] - f - 1; - break; - } - } - printf("Face is %d\n", i); - int vxPos = grd->face_nodepos[ face ]; - const int vxs[ 4 ] = { grd->face_nodes[ vxPos ], - grd->face_nodes[ vxPos+1 ], - grd->face_nodes[ vxPos+2 ], - grd->face_nodes[ vxPos+3 ] }; - for( int vx=0; vx<4; ++vx ) - { - grd->face_nodes[ vxPos+vx ] = vxs[ map[ i ][ vx ]]; - } - - /* - for( int vx=grd->face_nodepos[ f ]; vxface_nodepos[ f+1 ]; ++vx ) - { - printf(" %d ",grd->face_nodes[ vx ]); - } - printf(")\n"); - */ - } - } - } - for( int f=0; fnumber_of_faces; ++f ) { //for( int f=grd->cell_facepos[ c ]; fcell_facepos[ c+1 ]; ++f ) @@ -248,6 +196,7 @@ void print_grid(const struct UnstructuredGrid *grd ) printf("\n"); } + /* if( dim == 2 ) { for( int c=0; cnumber_of_cells; ++c ) @@ -257,6 +206,7 @@ void print_grid(const struct UnstructuredGrid *grd ) SWAP( grd->cell_facetag[ f+1 ], grd->cell_facetag[ f+2 ] ); } } + */ for( int c=0; cnumber_of_cells; ++c ) { diff --git a/opm/grid/cart_grid.c b/opm/grid/cart_grid.c index 10a0d5ed7..8226d1840 100644 --- a/opm/grid/cart_grid.c +++ b/opm/grid/cart_grid.c @@ -609,8 +609,8 @@ fill_cart_topology_2d(struct UnstructuredGrid *G) /* Faces with y-normal */ for (j=0; j GlobalCoordinate - corner( const EntitySeed& seed, const int i ) const + corner( const EntitySeed& seed, int i, const bool swap = false ) const { const int codim = EntitySeed :: codimension; switch (codim) @@ -1012,12 +1012,10 @@ namespace Dune } case 1: { - static const int map[ 6 ][ 4 ] = { {0, 1, 3, 2}, // face 0 - {0, 1, 3, 2}, // face 1 - {0, 3, 1, 2}, // face 2 - {0, 3, 1, 2}, // face 3 - {0, 1, 3, 2}, // face 4 - {0, 1, 3, 2}}; // face 5 + // for faces we need to swap vertices since in UnstructuredGrid + // those are ordered counter clockwise + if( EntitySeed :: dimension == 3 && i > 1 ) + i = 5 - i; const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + i ]; return copyToGlobalCoordinate( grid_.node_coordinates + GlobalCoordinate :: dimension * faceVertex ); @@ -1353,7 +1351,7 @@ namespace Dune protected: void init() { - std::cout << "PolyhedralGrid init" << std::endl; + //std::cout << "PolyhedralGrid init" << std::endl; // copy Cartesian dimensions for( int i=0; i<3; ++i ) @@ -1369,7 +1367,6 @@ namespace Dune // sort vertices such that they comply with the dune cube reference element if( grid_.cell_facetag ) { - std::cout << "Sort vertices " << std::endl; typedef std::array KeyType; std::map< const KeyType, const int > vertexFaceTags; const int vertexFacePattern [8][3] = { @@ -1396,6 +1393,14 @@ namespace Dune for (int c = 0; c < numCells; ++c) { + if( dim == 2 ) + { + // for 2d Cartesian grids the face ordering is wrong + int f = grid_.cell_facepos[ c ]; + std::swap( grid_.cell_faces[ f+1 ], grid_.cell_faces[ f+2 ] ); + std::swap( grid_.cell_facetag[ f+1 ], grid_.cell_facetag[ f+2 ] ); + } + typedef std::map vertexmap_t; typedef typename vertexmap_t :: iterator iterator; diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 216e40fe4..42f074e02 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -238,7 +238,6 @@ int main(int argc, char** argv ) //testGrid( *gridPtr, "polyhedralgrid-dgf", std::pow(3, int(Grid::dimension)) ); gridcheck( *gridPtr ); } - return 0; { std::cout <<"Check 2d grid" << std::endl; From af2ef5ab0e441c2eb8dab80273f0a87a82293e22 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Thu, 7 Nov 2019 17:04:25 +0100 Subject: [PATCH 21/22] Make compile in debug mode. --- opm/grid/polyhedralgrid/grid.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 5a10d6158..2bec3b1a9 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1156,7 +1156,7 @@ namespace Dune bool isBoundaryFace(const typename Codim<1>::EntitySeed& faceSeed ) const { - assert( faceSeed.valid() ); + assert( faceSeed.isValid() ); return isBoundaryFace( faceSeed.index() ); } From b4deed60e9d94f1b6b9ece99993f81f7f2421b45 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 8 Nov 2019 14:37:15 +0100 Subject: [PATCH 22/22] Enable corner point test for Polyhedral and remove warnings. --- opm/grid/polyhedralgrid/grid.hh | 34 ++++++++++++--------------------- tests/grid_test.cc | 22 ++++++++++----------- 2 files changed, 23 insertions(+), 33 deletions(-) diff --git a/opm/grid/polyhedralgrid/grid.hh b/opm/grid/polyhedralgrid/grid.hh index 2bec3b1a9..c737a09f8 100644 --- a/opm/grid/polyhedralgrid/grid.hh +++ b/opm/grid/polyhedralgrid/grid.hh @@ -1000,7 +1000,7 @@ namespace Dune template GlobalCoordinate - corner( const EntitySeed& seed, int i, const bool swap = false ) const + corner( const EntitySeed& seed, const int i ) const { const int codim = EntitySeed :: codimension; switch (codim) @@ -1012,12 +1012,10 @@ namespace Dune } case 1: { - // for faces we need to swap vertices since in UnstructuredGrid - // those are ordered counter clockwise - if( EntitySeed :: dimension == 3 && i > 1 ) - i = 5 - i; - - const int faceVertex = grid_.face_nodes[ grid_.face_nodepos[seed.index() ] + i ]; + // for faces we need to swap vertices in 3d since in UnstructuredGrid + // those are ordered counter clockwise, for 2d this does not matter + const int crner = (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: @@ -1509,8 +1507,8 @@ namespace Dune { assert( cellVertices_[ c ].size() == 4 ); GlobalCoordinate center( 0 ); - GlobalCoordinate p[ 4 ]; - for( int i=0; i<4; ++i ) + GlobalCoordinate p[ dim+1 ]; + for( int i=0; i matrix( 0 ); - matrix [0][0] = p[1][0] - p[0][0] ; - matrix [0][1] = p[1][1] - p[0][1] ; - matrix [0][2] = p[1][2] - p[0][2] ; - - matrix [1][0] = p[2][0] - p[0][0] ; - matrix [1][1] = p[2][1] - p[0][1] ; - matrix [1][2] = p[2][2] - p[0][2] ; - - matrix [2][0] = p[3][0] - p[0][0] ; - matrix [2][1] = p[3][1] - p[0][1] ; - matrix [2][2] = p[3][2] - p[0][2] ; + Dune::GeometryType simplex; + simplex.makeSimplex( dim ); - grid_.cell_volumes[ c ] = std::abs( matrix.determinant() )/ 6.0; + typedef Dune::AffineGeometry< ctype, dim, dimworld> AffineGeometryType; + AffineGeometryType geometry( simplex, p ); + grid_.cell_volumes[ c ] = geometry.volume(); } } diff --git a/tests/grid_test.cc b/tests/grid_test.cc index 42f074e02..0d1561220 100644 --- a/tests/grid_test.cc +++ b/tests/grid_test.cc @@ -33,16 +33,16 @@ const char *deckString = "RUNSPEC\n" "METRIC\n" "DIMENS\n" - "2 2 2 /\n" + "4 4 4 /\n" "GRID\n" "DXV\n" - "2*1 /\n" + "4*1 /\n" "DYV\n" - "2*1 /\n" + "4*1 /\n" "DZ\n" - "8*1 /\n" + "16*1 /\n" "TOPS\n" - "4*100.0 /\n"; + "16*100.0 /\n"; template void testGridIteration( const GridView& gridView, const int nElem ) @@ -204,20 +204,20 @@ int main(int argc, char** argv ) //dgfFile << "#" << std::endl; #if HAVE_ECL_INPUT - /* Opm::Parser parser; const auto deck = parser.parseString(deckString); std::vector porv; - */ #endif // test PolyhedralGrid { typedef Dune::PolyhedralGrid< 3, 3 > Grid; #if HAVE_ECL_INPUT - /* + std::cout << "Check ecl grid" << std::endl; Grid grid(deck, porv); - testGrid( grid, "polyhedralgrid", 8, 27 ); + gridcheck( grid ); + //testGrid( grid, "polyhedralgrid", 8, 27 ); + /* Opm::TopSurf* ts; ts = Opm::TopSurf::create (grid); std::cout << ts->dimensions << std::endl; @@ -259,7 +259,7 @@ int main(int argc, char** argv ) } } - /* +#if 0 // test CpGrid { typedef Dune::CpGrid Grid; @@ -278,6 +278,6 @@ int main(int argc, char** argv ) //Dune::GridPtr< Grid > gridPtr( dgfFile ); //testGrid( *gridPtr, "cpgrid-dgf" ); } - */ +#endif return 0; }