diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index bc0f23a9..c829fd5a 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -34,17 +34,6 @@ list (APPEND MAIN_SOURCE_FILES opm/core/flowdiagnostics/FlowDiagnostics.cpp opm/core/flowdiagnostics/TofDiscGalReorder.cpp opm/core/flowdiagnostics/TofReorder.cpp - opm/core/grid/GridHelpers.cpp - opm/core/grid/GridManager.cpp - opm/core/grid/GridUtilities.cpp - opm/core/grid/cart_grid.c - opm/core/grid/cornerpoint_grid.c - opm/core/grid/cpgpreprocess/facetopology.c - opm/core/grid/cpgpreprocess/geometry.c - opm/core/grid/cpgpreprocess/preprocess.c - opm/core/grid/cpgpreprocess/uniquepoints.c - opm/core/grid/grid.c - opm/core/grid/grid_equal.cpp opm/core/linalg/LinearSolverFactory.cpp opm/core/linalg/LinearSolverInterface.cpp opm/core/linalg/LinearSolverIstl.cpp @@ -112,7 +101,6 @@ list (APPEND MAIN_SOURCE_FILES opm/core/utility/Event.cpp opm/core/utility/MonotCubicInterpolator.cpp opm/core/utility/NullStream.cpp - opm/core/utility/StopWatch.cpp opm/core/utility/VelocityInterpolation.cpp opm/core/utility/WachspressCoord.cpp opm/core/utility/compressedToCartesian.cpp @@ -150,7 +138,6 @@ list (APPEND TEST_SOURCE_FILES tests/test_nonuniformtablelinear.cpp tests/test_parallelistlinformation.cpp tests/test_sparsevector.cpp - tests/test_sparsetable.cpp tests/test_velocityinterpolation.cpp tests/test_quadratures.cpp tests/test_uniformtablelinear.cpp @@ -258,21 +245,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/flowdiagnostics/FlowDiagnostics.hpp opm/core/flowdiagnostics/TofDiscGalReorder.hpp opm/core/flowdiagnostics/TofReorder.hpp - opm/core/grid.h - opm/core/grid/CellQuadrature.hpp - opm/core/grid/ColumnExtract.hpp - opm/core/grid/FaceQuadrature.hpp - opm/core/grid/GridHelpers.hpp - opm/core/grid/GridManager.hpp - opm/core/grid/GridUtilities.hpp - opm/core/grid/MinpvProcessor.hpp - opm/core/grid/PinchProcessor.hpp - opm/core/grid/cart_grid.h - opm/core/grid/cornerpoint_grid.h - opm/core/grid/cpgpreprocess/facetopology.h - opm/core/grid/cpgpreprocess/geometry.h - opm/core/grid/cpgpreprocess/preprocess.h - opm/core/grid/cpgpreprocess/uniquepoints.h opm/core/linalg/LinearSolverFactory.hpp opm/core/linalg/LinearSolverInterface.hpp opm/core/linalg/LinearSolverIstl.hpp @@ -383,9 +355,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/utility/NullStream.hpp opm/core/utility/RegionMapping.hpp opm/core/utility/RootFinders.hpp - opm/core/utility/SparseTable.hpp opm/core/utility/SparseVector.hpp - opm/core/utility/StopWatch.hpp opm/core/utility/UniformTableLinear.hpp opm/core/utility/VelocityInterpolation.hpp opm/core/utility/WachspressCoord.hpp diff --git a/dune.module b/dune.module index 71f8d96f..c711df2c 100644 --- a/dune.module +++ b/dune.module @@ -10,4 +10,4 @@ Label: 2017.04-pre Maintainer: atgeirr@sintef.no MaintainerName: Atgeirr F. Rasmussen Url: http://opm-project.org -Depends: dune-common (>= 2.2) dune-istl (>= 2.2) opm-common opm-parser opm-output opm-material +Depends: dune-common (>= 2.2) dune-istl (>= 2.2) opm-common opm-parser opm-grid opm-output opm-material diff --git a/opm/core/grid.h b/opm/core/grid.h deleted file mode 100644 index c2f8a543..00000000 --- a/opm/core/grid.h +++ /dev/null @@ -1,332 +0,0 @@ -/* - Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_GRID_HEADER_INCLUDED -#define OPM_GRID_HEADER_INCLUDED - -#include -#include - -/** - * \file - * - * Main OPM-Core grid data structure along with helper functions for - * construction, destruction and reading a grid representation from disk. - */ - -#ifdef __cplusplus -extern "C" { -#endif - -/* ----- synopsis of grid.h ---- - -struct UnstructuredGrid -{ - int dimensions; - int number_of_cells; - int number_of_faces; - int number_of_nodes; - int *face_nodes; - int *face_nodepos; - int *face_cells; - int *cell_faces; - int *cell_facepos; - double *node_coordinates; - double *face_centroids; - double *face_areas; - double *face_normals; - double *cell_centroids; - double *cell_volumes; - int *global_cell; - int cartdims[3]; - int *cell_facetag; -}; - -void destroy_grid(struct UnstructuredGrid *g); - -struct UnstructuredGrid * -create_grid_empty(void); - -struct UnstructuredGrid * -allocate_grid(size_t ndims , - size_t ncells , - size_t nfaces , - size_t nfacenodes, - size_t ncellfaces, - size_t nnodes ); - -struct UnstructuredGrid * -read_grid(const char *fname); - - ---- end of synopsis of grid.h ---- -*/ - -/** - Data structure for an unstructured grid, unstructured meaning that - any cell may have an arbitrary number of adjacent cells. The struct - contains both topological and geometrical data. - - The grid consists of a set of cells, which are assumed to partion - the grid domain. A face is defined as the nonempty intersection of - (the closure of) two grid cells (the grid is a cell complex). - - The topology information is limited to some adjacency relations - between cells, faces and nodes only. The data structure does not - contain any information pertaining to edges (except in 2d, where - edges are the same as faces). - - The geometry information is limited to centroids, areas/volumes and - normals. - */ -struct UnstructuredGrid -{ - /** - The topological and geometrical dimensionality of the - grid. Note that we do not support grids that are embedded in - higher-dimensional spaces, such as 2d grids embedded in 3d. - This number must be 2 or 3. - */ - int dimensions; - - /** The number of cells in the grid. */ - int number_of_cells; - /** The number of faces in the grid. */ - int number_of_faces; - /** The number of nodes in the grid. */ - int number_of_nodes; - - /** - Contains for each face, the indices of its adjacent nodes. - The size of the array is equal to the sum over all faces of - each face's number of adjacent nodes, which also is equal to - face_nodepos[number_of_faces]. - */ - int *face_nodes; - /** - For a face f, face_nodepos[f] contains the starting index - for f's nodes in the face_nodes array. - The size of the array is equal to (number_of_faces + 1). - */ - int *face_nodepos; - /** - For a face f, face_cells[2*f] and face_cells[2*f + 1] contain - the cell indices of the cells adjacent to f. The number -1 - indicates the outer boundary. - The order is significant, as it gives the orientation: if - face_cells[2*f] == a and face_cells[2*f + 1] == b, f is - oriented from a to b. The inverse of this mapping is stored in - cell_faces and cell_facepos. - The size of the array is equal to (2*number_of_faces). - */ - int *face_cells; - - /** - Contains for each cell, the indices of its adjacent faces. - The size of the array is equal to the sum over all cells of - each cell's number of adjacent faces, which also is equal to - cell_facepos[number_of_cells]. - */ - int *cell_faces; - /** - For a cell c, cell_facepos[c] contains the starting index - for c's faces in the cell_faces array. - The size of the array is equal to (number_of_cells + 1). - */ - int *cell_facepos; - - /** - Node coordinates, stored consecutively for each node. That is, - for a node i, node_coordinates[dimensions*i + d] contains the - d'th coordinate of node i. - The size of the array is equal to (dimensions*number_of_nodes). - */ - double *node_coordinates; - - /** - Exact or approximate face centroids, stored consecutively for each face. That is, - for a face f, face_centroids[dimensions*f + d] contains the - d'th coordinate of f's centroid. - The size of the array is equal to (dimensions*number_of_faces). - */ - double *face_centroids; - /** - Exact or approximate face areas. - The size of the array is equal to number_of_faces. - */ - double *face_areas; - /** - Exact or approximate face normals, stored consecutively for - each face. That is, for a face f, face_normals[dimensions*f + d] - contains the d'th coordinate of f's normal. - The size of the array is equal to (dimensions*number_of_faces). - - IMPORTANT: the normals are not normalized to have unit length! - They are assumed to always have length equal to the - corresponding face's area. - */ - double *face_normals; - - /** - Exact or approximate cell centroids, stored consecutively for each cell. That is, - for a cell c, cell_centroids[dimensions*c + d] contains the - d'th coordinate of c's centroid. - The size of the array is equal to (dimensions*number_of_cells). - */ - double *cell_centroids; - /** - Exact or approximate cell volumes. - The size of the array is equal to number_of_cells. - */ - double *cell_volumes; - - - /** - If non-null, this array contains the logical cartesian indices - (in a lexicographic ordering) of each cell. - In that case, the array size is equal to number_of_cells. - - This field is intended for grids that have a (possibly - degenerate) logical cartesian structure, for example - cornerpoint grids. - - If null, this indicates that the element indices coincide - with the logical cartesian indices, _or_ that the grid has - no inherent Cartesian structure. Due to this ambiguity, this - field should not be used to check if the grid is Cartesian. - */ - int *global_cell; - - /** - Contains the size of the logical cartesian structure (if any) of the grid. - - This field is intended for grids that have a (possibly - degenerate) logical cartesian structure, for example - cornerpoint grids. - - If the grid is unstructured (non-Cartesian), then at least one - of the items in the (sub-)array cartdims[0..dimensions-1] - _could_ have the value 0 to signal this. - */ - int cartdims[3]; - /** - If non-null, this array contains a number for cell-face - adjacency indicating the face's position with respect to the - cell, in a logical cartesian sense. The tags are in [0, ..., 5] - meaning [I-, I+, J-, J+, K-, K+], where I, J, K are the logical - cartesian principal directions. - The structure of this array is identical to cell_faces, and - cell_facepos indices into cell_facetag as well. - - If non-null, the array size is equal to - cell_facepos[number_of_cells]. - - This field is intended for grids that have a (possibly - degenerate) logical cartesian structure, for example - cornerpoint grids. - */ - int *cell_facetag; - - - /* - This vector is retained to be able to construct an - EclipseGrid representation of the Grid. If the grid - processor actually modifies the elements of the zcorn - vector from the input the modified version is stored here; - otherwise we just use the default. - */ - double * zcorn; -}; - -/** - Destroy and deallocate an UnstructuredGrid and all its data. - - This function assumes that all arrays of the UnstructuredGrid (if - non-null) have been individually allocated by malloc(). They will - be deallocated with free(). - */ -void destroy_grid(struct UnstructuredGrid *g); - -/** - Allocate and initialise an empty UnstructuredGrid. - - This is the moral equivalent of a C++ default constructor. - - \return Fully formed UnstructuredGrid with all fields zero or - NULL as appropriate. NULL in case of - allocation failure. - */ -struct UnstructuredGrid * -create_grid_empty(void); - -/** - Allocate and initialise an UnstructuredGrid where pointers are set - to location with correct size. - - \param[in] ndims Number of physical dimensions. - \param[in] ncells Number of cells. - \param[in] nfaces Number of faces. - \param[in] nfacenodes Size of mapping from faces to nodes. - \param[in] ncellfaces Size of mapping from cells to faces. - (i.e., the number of `half-faces') - \param[in] nnodes Number of Nodes. - - \return Fully formed UnstructuredGrid with all fields except - global_cell allocated, but not filled with meaningful - values. NULL in case of allocation failure. - */ -struct UnstructuredGrid * -allocate_grid(size_t ndims , - size_t ncells , - size_t nfaces , - size_t nfacenodes, - size_t ncellfaces, - size_t nnodes ); - - -/** - Will allocate storage internally in the grid object to hold a copy - of the zcorn data supplied in the second argument. -*/ -void -attach_zcorn_copy(struct UnstructuredGrid* G , - const double * zcorn); - - -/** - * Import a grid from a character representation stored in file. - * - * @param[in] fname File name. - * @return Fully formed UnstructuredGrid with all fields allocated and filled. - * Returns @c NULL in case of allocation failure. - */ -struct UnstructuredGrid * -read_grid(const char *fname); - - - - -bool -grid_equal(const struct UnstructuredGrid * grid1 , const struct UnstructuredGrid * grid2); - -#ifdef __cplusplus -} -#endif - -#endif /* OPM_GRID_HEADER_INCLUDED */ diff --git a/opm/core/grid/CellQuadrature.hpp b/opm/core/grid/CellQuadrature.hpp deleted file mode 100644 index ab0eb043..00000000 --- a/opm/core/grid/CellQuadrature.hpp +++ /dev/null @@ -1,260 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_CELLQUADRATURE_HEADER_INCLUDED -#define OPM_CELLQUADRATURE_HEADER_INCLUDED - -#include -#include -#include -#include - -namespace Opm -{ - - - namespace { - /// Calculates the determinant of a 3 x 3 matrix, represented as - /// three three-dimensional arrays. - inline double determinantOf(const double* a0, - const double* a1, - const double* a2) - { - return - a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - - a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + - a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); - } - - /// Computes the volume of a tetrahedron consisting of 4 vertices - /// with 3-dimensional coordinates - inline double tetVolume(const double* p0, - const double* p1, - const double* p2, - const double* p3) - { - double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; - double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; - double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] }; - return std::fabs(determinantOf(a, b, c) / 6.0); - } - - /// Calculates the area of a triangle consisting of 3 vertices - /// with 2-dimensional coordinates - inline double triangleArea2d(const double* p0, - const double* p1, - const double* p2) - { - double a[2] = { p1[0] - p0[0], p1[1] - p0[1] }; - double b[2] = { p2[0] - p0[0], p2[1] - p0[1] }; - double a_cross_b = a[0]*b[1] - a[1]*b[0]; - return 0.5*std::fabs(a_cross_b); - } - } // anonymous namespace - - - /// A class providing numerical quadrature for cells. - /// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). - /// Note that this class does multiply weights by cell volume, - /// so weights always sum to cell volume. - /// - /// Degree 1 method: - /// Midpoint (centroid) method. - /// n = 1, w_0 = cell volume, x_0 = cell centroid - /// - /// Degree 2 method for 2d (but see the note): - /// Based on subdivision of the cell into triangles, - /// with the centroid as a common vertex, and the triangle - /// edge midpoint rule. - /// Triangle i consists of the centroid C, nodes N_i and N_{i+1}. - /// Its area is A_i. - /// n = 2 * nn (nn = num nodes in face) - /// For i = 0..(nn-1): - /// w_i = 1/3 A_i. - /// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i - /// x_i = (N_i + N_{i+1})/2 - /// x_{nn+i} = (C + N_i)/2 - /// All N and A indices are interpreted cyclic, modulus nn. - /// Note: for simplicity of implementation, we currently use - /// n = 3 * nn - /// For i = 0..(nn-1): - /// w_{3*i + {0,1,2}} = 1/3 A_i - /// x_{3*i} = (N_i + N_{i+1})/2 - /// x_{3*i + {1,2}} = (C + N_{i,i+1})/2 - /// This is simpler, because we can implement it easily - /// based on iteration over faces without requiring any - /// particular (cyclic) ordering. - /// - /// Degree 2 method for 3d: - /// Based on subdivision of each cell face into triangles - /// with the face centroid as a common vertex, and then - /// subdividing the cell into tetrahedra with the cell - /// centroid as a common vertex. Then apply the tetrahedron - /// rule with the following 4 nodes (uniform weights): - /// a = 0.138196601125010515179541316563436 - /// x_i has all barycentric coordinates = a, except for - /// the i'th coordinate which is = 1 - 3a. - /// This rule is from http://nines.cs.kuleuven.be/ecf, - /// it is the second degree 2 4-point rule for tets, - /// referenced to Stroud(1971). - /// The tetrahedra are numbered T_{i,j}, and are given by the - /// cell centroid C, the face centroid FC_i, and two nodes - /// of face i: FN_{i,j}, FN_{i,j+1}. - class CellQuadrature - { - public: - CellQuadrature(const UnstructuredGrid& grid, - const int cell, - const int degree) - : grid_(grid), cell_(cell), degree_(degree) - { - if (grid.dimensions > 3) { - OPM_THROW(std::runtime_error, "CellQuadrature only implemented for up to 3 dimensions."); - } - if (degree > 2) { - OPM_THROW(std::runtime_error, "CellQuadrature exact for polynomial degrees > 1 not implemented."); - } - } - - int numQuadPts() const - { - if (degree_ < 2 || grid_.dimensions == 1) { - return 1; - } - // Degree 2 case. - if (grid_.dimensions == 2) { - return 3*(grid_.cell_facepos[cell_ + 1] - grid_.cell_facepos[cell_]); - } - assert(grid_.dimensions == 3); - int sumnodes = 0; - for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { - const int face = grid_.cell_faces[hf]; - sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; - } - return 4*sumnodes; - } - - void quadPtCoord(const int index, double* coord) const - { - const int dim = grid_.dimensions; - const double* cc = grid_.cell_centroids + dim*cell_; - if (degree_ < 2) { - std::copy(cc, cc + dim, coord); - return; - } - // Degree 2 case. - if (dim == 2) { - if (index % 3 == 0) { - // Boundary midpoint. This is the face centroid. - const int hface = grid_.cell_facepos[cell_] + index/3; - const int face = grid_.cell_faces[hface]; - const double* fc = grid_.face_centroids + dim*face; - std::copy(fc, fc + dim, coord); - } else { - // Interiour midpoint. This is the average of the - // cell centroid and a face node (they should - // always have two nodes in 2d). - const int hface = grid_.cell_facepos[cell_] + index/3; - const int face = grid_.cell_faces[hface]; - const int nodeoff = (index % 3) - 1; // == 0 or 1 - const int node = grid_.face_nodes[grid_.face_nodepos[face] + nodeoff]; - const double* nc = grid_.node_coordinates + dim*node; - for (int dd = 0; dd < dim; ++dd) { - coord[dd] = 0.5*(nc[dd] + cc[dd]); - } - } - return; - } - assert(dim == 3); - int tetindex = index / 4; - const int subindex = index % 4; - const double* nc = grid_.node_coordinates; - for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { - const int face = grid_.cell_faces[hf]; - const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; - if (nfn <= tetindex) { - // Our tet is not associated with this face. - tetindex -= nfn; - continue; - } - const double* fc = grid_.face_centroids + dim*face; - const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; - const int node0 = fnodes[tetindex]; - const int node1 = fnodes[(tetindex + 1) % nfn]; - const double* n0c = nc + dim*node0; - const double* n1c = nc + dim*node1; - const double a = 0.138196601125010515179541316563436; - // Barycentric coordinates of our point in the tet. - double baryc[4] = { a, a, a, a }; - baryc[subindex] = 1.0 - 3.0*a; - for (int dd = 0; dd < dim; ++dd) { - coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd]; - } - return; - } - OPM_THROW(std::runtime_error, "Should never reach this point."); - } - - double quadPtWeight(const int index) const - { - if (degree_ < 2) { - return grid_.cell_volumes[cell_]; - } - // Degree 2 case. - const int dim = grid_.dimensions; - const double* cc = grid_.cell_centroids + dim*cell_; - if (dim == 2) { - const int hface = grid_.cell_facepos[cell_] + index/3; - const int face = grid_.cell_faces[hface]; - const int* nptr = grid_.face_nodes + grid_.face_nodepos[face]; - const double* nc0 = grid_.node_coordinates + dim*nptr[0]; - const double* nc1 = grid_.node_coordinates + dim*nptr[1]; - return triangleArea2d(nc0, nc1, cc)/3.0; - } - assert(dim == 3); - int tetindex = index / 4; - const double* nc = grid_.node_coordinates; - for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { - const int face = grid_.cell_faces[hf]; - const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; - if (nfn <= tetindex) { - // Our tet is not associated with this face. - tetindex -= nfn; - continue; - } - const double* fc = grid_.face_centroids + dim*face; - const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; - const int node0 = fnodes[tetindex]; - const int node1 = fnodes[(tetindex + 1) % nfn]; - const double* n0c = nc + dim*node0; - const double* n1c = nc + dim*node1; - return 0.25*tetVolume(cc, fc, n0c, n1c); - } - OPM_THROW(std::runtime_error, "Should never reach this point."); - } - - private: - const UnstructuredGrid& grid_; - const int cell_; - const int degree_; - }; - -} // namespace Opm - -#endif // OPM_CELLQUADRATURE_HEADER_INCLUDED diff --git a/opm/core/grid/ColumnExtract.hpp b/opm/core/grid/ColumnExtract.hpp deleted file mode 100644 index 9bb4b76a..00000000 --- a/opm/core/grid/ColumnExtract.hpp +++ /dev/null @@ -1,122 +0,0 @@ -#include -#include -#include -#include - -namespace Opm { - - namespace { - - /// Helper struct for extractColumn - /// Compares the underlying k-index - struct ExtractColumnCompare - { - ExtractColumnCompare(const UnstructuredGrid& g) - : grid(g) - { - // empty - } - - bool operator()(const int i, const int j) - { - // Extract k-index - int index_i = grid.global_cell ? grid.global_cell[i] : i; - int k_i = index_i / grid.cartdims[0] / grid.cartdims[1]; - int index_j = grid.global_cell ? grid.global_cell[j] : j; - int k_j = index_j / grid.cartdims[0] / grid.cartdims[1]; - - return k_i < k_j; - } - - const UnstructuredGrid& grid; - }; - - - /// Neighbourhood query. - /// \return true if two cells are neighbours. - bool neighbours(const UnstructuredGrid& grid, const int c0, const int c1) - { - for (int hf = grid.cell_facepos[c0]; hf < grid.cell_facepos[c0 + 1]; ++hf) { - const int f = grid.cell_faces[hf]; - if (grid.face_cells[2*f] == c1 || grid.face_cells[2*f+1] == c1) { - return true; - } - } - return false; - } - - } // anonymous namespace - - -/// Extract each column of the grid. -/// \note Assumes the pillars of the grid are all vertically aligned. -/// \param grid The grid from which to extract the columns. -/// \param columns will for each (i, j) where (i, j) represents a non-empty column, -//// contain the cell indices contained in the column -/// centered at (i, j) in the second variable, and i+jN in the first variable. -inline void extractColumn( const UnstructuredGrid& grid, std::vector >& columns ) -{ - const int* dims = grid.cartdims; - - // Keeps track of column_index ---> index of vector - std::map global_to_local; - for (int cell = 0; cell < grid.number_of_cells; ++cell) { - // Extract Cartesian coordinates - int index = grid.global_cell ? grid.global_cell[cell] : cell; // If null, assume mapping is identity. - int i_cart = index % dims[0]; - int k_cart = index / dims[0] / dims[1]; - int j_cart = (index - k_cart*dims[0]*dims[1])/ dims[0]; - - int local_index; - std::map::iterator local_index_iterator = global_to_local.find(i_cart+j_cart*dims[0]); - if (local_index_iterator != global_to_local.end()) { - local_index = local_index_iterator->second; - } else { - local_index = columns.size(); - global_to_local[i_cart+j_cart*dims[0]] = local_index; - columns.push_back(std::vector()); - } - columns[local_index].push_back(cell); - } - - int num_cols = columns.size(); - for (int col = 0; col < num_cols; ++col) { - std::sort(columns[col].begin(), columns[col].end(), ExtractColumnCompare(grid)); - } - - // At this point, a column may contain multiple disjoint sets of cells. - // We must split these columns into connected parts. - std::vector< std::vector > new_columns; - for (int col = 0; col < num_cols; ++col) { - const int colsz = columns[col].size(); - int first_of_col = 0; - for (int k = 1; k < colsz; ++k) { - const int c0 = columns[col][k-1]; - const int c1 = columns[col][k]; - if (!neighbours(grid, c0, c1)) { - // Must split. Move the cells [first_of_col, ... , k-1] to - // a new column, known to be connected. - new_columns.push_back(std::vector()); - new_columns.back().assign(columns[col].begin() + first_of_col, columns[col].begin() + k); - // The working column now starts with index k. - first_of_col = k; - } - } - if (first_of_col != 0) { - // The column was split, the working part should be - // the entire column. We erase the cells before first_of_col. - // (Could be more efficient if we instead chop off end.) - columns[col].erase(columns[col].begin(), columns[col].begin() + first_of_col); - } - } - - // Must tack on the new columns to complete the set. - const int num_cols_all = num_cols + new_columns.size(); - columns.resize(num_cols_all); - for (int col = num_cols; col < num_cols_all; ++col) { - columns[col].swap(new_columns[col - num_cols]); - } - -} - -} // namespace Opm diff --git a/opm/core/grid/FaceQuadrature.hpp b/opm/core/grid/FaceQuadrature.hpp deleted file mode 100644 index 571b2f01..00000000 --- a/opm/core/grid/FaceQuadrature.hpp +++ /dev/null @@ -1,189 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_FACEQUADRATURE_HEADER_INCLUDED -#define OPM_FACEQUADRATURE_HEADER_INCLUDED - -#include -#include -#include - -namespace Opm -{ - - - namespace { - /// Calculates the cross product of two 3-vectors, represented as - /// 3-element arrays. Calculates res = a X b. The res array must - /// already be allocated with room for 3 elements. - inline void cross(const double* a, const double* b, double* res) - { - res[0] = a[1]*b[2] - a[2]*b[1]; - res[1] = a[2]*b[0] - a[0]*b[2]; - res[2] = a[0]*b[1] - a[1]*b[0]; - } - - /// Calculates the area of a triangle consisting of 3 vertices - /// with 3-dimensional coordinates - inline double triangleArea3d(const double* p0, - const double* p1, - const double* p2) - { - double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; - double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; - double cr[3]; - cross(a, b, cr); - return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]); - } - } // anonymous namespace - - - /// A class providing numerical quadrature for faces. - /// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). - /// Note that this class does multiply weights by face area, - /// so weights always sum to face area. - /// - /// Degree 1 method: - /// Midpoint (centroid) method. - /// n = 1, w_0 = face area, x_0 = face centroid - /// - /// Degree 2 method for 2d: - /// Simpson's method (actually this is degree 3). - /// - /// Degree 2 method for 3d: - /// Based on subdivision of the face into triangles, - /// with the centroid as a common vertex, and the triangle - /// edge midpoint rule. - /// Triangle i consists of the centroid C, nodes N_i and N_{i+1}. - /// Its area is A_i. - /// n = 2 * nn (nn = num nodes in face) - /// For i = 0..(nn-1): - /// w_i = 1/3 A_i. - /// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i - /// x_i = (N_i + N_{i+1})/2 - /// x_{nn+i} = (C + N_i)/2 - /// All N and A indices are interpreted cyclic, modulus nn. - class FaceQuadrature - { - public: - FaceQuadrature(const UnstructuredGrid& grid, - const int face, - const int degree) - : grid_(grid), face_(face), degree_(degree) - { - if (grid_.dimensions > 3) { - OPM_THROW(std::runtime_error, "FaceQuadrature only implemented for up to 3 dimensions."); - } - if (degree_ > 2) { - OPM_THROW(std::runtime_error, "FaceQuadrature exact for polynomial degrees > 2 not implemented."); - } - } - - int numQuadPts() const - { - if (degree_ < 2 || grid_.dimensions < 2) { - return 1; - } - // Degree 2 case. - if (grid_.dimensions == 2) { - return 3; - } else { - return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]); - } - } - - void quadPtCoord(const int index, double* coord) const - { - const int dim = grid_.dimensions; - const double* fc = grid_.face_centroids + dim*face_; - if (degree_ < 2 || dim < 2) { - std::copy(fc, fc + dim, coord); - return; - } - // Degree 2 case. - const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; - const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; - const double* nc = grid_.node_coordinates; - if (dim == 2) { - assert(nn == 2); - const double* pa[3] = { nc + dim*fnodes[0], fc, nc + dim*fnodes[1] }; - std::copy(pa[index], pa[index] + dim, coord); - return; - } - assert(dim == 3); - if (index < nn) { - // Boundary edge midpoint. - const int node0 = fnodes[index]; - const int node1 = fnodes[(index + 1)%nn]; - for (int dd = 0; dd < dim; ++dd) { - coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]); - } - } else { - // Interiour edge midpoint. - // Recall that index is now in [nn, 2*nn). - const int node = fnodes[index - nn]; - for (int dd = 0; dd < dim; ++dd) { - coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]); - } - } - } - - double quadPtWeight(const int index) const - { - if (degree_ < 2) { - return grid_.face_areas[face_]; - } - // Degree 2 case. - const int dim = grid_.dimensions; - if (dim == 2) { - const double simpsonw[3] = { 1.0/6.0, 4.0/6.0, 1.0/6.0 }; - return grid_.face_areas[face_]*simpsonw[index]; - } - assert(dim == 3); - const double* fc = grid_.face_centroids + dim*face_; - const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; - const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; - const double* nc = grid_.node_coordinates; - if (index < nn) { - // Boundary edge midpoint. - const int node0 = fnodes[index]; - const int node1 = fnodes[(index + 1)%nn]; - const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); - return area / 3.0; - } else { - // Interiour edge midpoint. - // Recall that index is now in [nn, 2*nn). - const int node0 = fnodes[(index - 1) % nn]; - const int node1 = fnodes[index - nn]; - const int node2 = fnodes[(index + 1) % nn]; - const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); - const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc); - return (area0 + area1) / 3.0; - } - } - - private: - const UnstructuredGrid& grid_; - const int face_; - const int degree_; - }; - -} // namespace Opm - -#endif // OPM_FACEQUADRATURE_HEADER_INCLUDED diff --git a/opm/core/grid/GridHelpers.cpp b/opm/core/grid/GridHelpers.cpp deleted file mode 100644 index 9b82507a..00000000 --- a/opm/core/grid/GridHelpers.cpp +++ /dev/null @@ -1,177 +0,0 @@ -/* - Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services. - Copyright 2014 Statoil AS - Copyright 2015 NTNU - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -#include "config.h" -#include -namespace Opm -{ -namespace UgGridHelpers -{ -int numCells(const UnstructuredGrid& grid) -{ - return grid.number_of_cells; -} - -int numFaces(const UnstructuredGrid& grid) -{ - return grid.number_of_faces; -} -int dimensions(const UnstructuredGrid& grid) -{ - return grid.dimensions; -} -int numCellFaces(const UnstructuredGrid& grid) -{ - return grid.cell_facepos[grid.number_of_cells]; -} - -const int* globalCell(const UnstructuredGrid& grid) -{ - return grid.global_cell; -} - -const int* cartDims(const UnstructuredGrid& grid) -{ - return grid.cartdims; -} - -const double* beginCellCentroids(const UnstructuredGrid& grid) -{ - return grid.cell_centroids; -} - -double cellCenterDepth(const UnstructuredGrid& grid, int cell_index) -{ - // This method is an alternative to the method cellCentroidCoordinate(...) below. - // The cell center depth is computed as a raw average of cell corner depths. - // For cornerpoint grids, this is likely to give slightly different depths that seem - // to agree with eclipse. - assert(grid.dimensions == 3); - const int nd = 3; // Assuming 3-dimensional grid ... - const int nv = 8; // Assuming 2*4 vertices ... - double zz = 0.0; - // Traverse the bottom and top cell-face - for (int i=grid.cell_facepos[cell_index+1]-2; i::const_iterator face) -{ - return grid.cell_facetag[face-cell2Faces(grid)[0].begin()]; -} - -SparseTableView cell2Faces(const UnstructuredGrid& grid) -{ - return SparseTableView(grid.cell_faces, grid.cell_facepos, numCells(grid)); -} - -SparseTableView face2Vertices(const UnstructuredGrid& grid) -{ - return SparseTableView(grid.face_nodes, grid.face_nodepos, numFaces(grid)); -} - -const double* vertexCoordinates(const UnstructuredGrid& grid, int index) -{ - return grid.node_coordinates+dimensions(grid)*index; -} - -double cellVolume(const UnstructuredGrid& grid, int cell_index) -{ - return grid.cell_volumes[cell_index]; -} - -FaceCellTraits::Type faceCells(const UnstructuredGrid& grid) -{ - return FaceCellsProxy(grid); -} - - -Opm::EclipseGrid createEclipseGrid(const UnstructuredGrid& grid, const Opm::EclipseGrid& inputGrid ) { - const int * dims = UgGridHelpers::cartDims( grid ); - - if ((inputGrid.getNX( ) == static_cast(dims[0])) && - (inputGrid.getNY( ) == static_cast(dims[1])) && - (inputGrid.getNZ( ) == static_cast(dims[2]))) { - std::vector updatedACTNUM; - const int* global_cell = UgGridHelpers::globalCell( grid ); - - if (global_cell) { - updatedACTNUM.assign( inputGrid.getCartesianSize( ) , 0 ); - for (int c = 0; c < numCells( grid ); c++) { - updatedACTNUM[global_cell[c]] = 1; - } - } - - return Opm::EclipseGrid( inputGrid, grid.zcorn, updatedACTNUM ); - } else { - throw std::invalid_argument("Size mismatch - dimensions of inputGrid argument and current UnstructuredGrid instance disagree"); - } -} - -} -} diff --git a/opm/core/grid/GridHelpers.hpp b/opm/core/grid/GridHelpers.hpp deleted file mode 100644 index cd89e416..00000000 --- a/opm/core/grid/GridHelpers.hpp +++ /dev/null @@ -1,381 +0,0 @@ -/* - Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services - Copyright 2014 Statoil AS - Copyright 2015 - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -#ifndef OPM_CORE_GRIDHELPERS_HEADER_INCLUDED -#define OPM_CORE_GRIDHELPERS_HEADER_INCLUDED - -#include -#include - -#include -#include -#include - - -namespace Opm -{ -namespace UgGridHelpers -{ - -/// \brief Allows viewing a sparse table consisting out of C-array -/// -/// This class can be used to convert two int array (like they are -/// in UnstructuredGrid for representing the cell to faces mapping -/// as a sparse table object. -class SparseTableView -{ -public: - class IntRange : public boost::iterator_range - { - public: - typedef boost::iterator_range BaseRowType; - typedef BaseRowType::size_type size_type; - typedef int value_type; - - IntRange(const int* start_arg, const int* end_arg) - : BaseRowType(start_arg, end_arg) - {} - }; - /// \brief The type of the roww. - typedef boost::iterator_range row_type; - - /// \brief Creates a sparse table view - /// \param data The array with data of the table. - /// \param offset The offsets of the rows. Row i starts - /// at offset[i] and ends a offset[i+1] - /// \param size The number of entries/rows of the table - SparseTableView(int* data, int *offset, std::size_t size_arg) - : data_(data), offset_(offset), size_(size_arg) - {} - - /// \brief Get a row of the the table. - /// \param row The row index. - /// \return The corresponding row. - row_type operator[](std::size_t row) const - { - assert(row<=size()); - return row_type(data_ + offset_[row], data_ + offset_[row+1]); - } - - /// \brief Get the size of the table. - /// \return the number rows. - std::size_t size() const - { - return size_; - } - - /// \brief Get the number of non-zero entries. - std::size_t noEntries() const - { - return offset_[size_]; - } - -private: - /// \brief The array with data of the table. - const int* data_; - /// \brief offset The offsets of the rows. - /// - /// Row i starts at offset[i] and ends a offset[i+1] - const int* offset_; - /// \brief The size, i.e. the number of rows. - std::size_t size_; -}; - -/// \brief Get the number of cells of a grid. -int numCells(const UnstructuredGrid& grid); - -/// \brief Get the number of faces of a grid. -int numFaces(const UnstructuredGrid& grid); - -/// \brief Get the dimensions of a grid -int dimensions(const UnstructuredGrid& grid); - -/// \brief Get the number of faces, where each face counts as many times as there are adjacent faces -int numCellFaces(const UnstructuredGrid& grid); - -/// \brief Get the cartesion dimension of the underlying structured grid. -const int* cartDims(const UnstructuredGrid& grid); - -/// \brief Get the local to global index mapping. -/// -/// The global index is the index of the active cell -/// in the underlying structured grid. -const int* globalCell(const UnstructuredGrid& grid); - -/// \brief Traits of the cell centroids of a grid. -/// -/// This class exports two types: IteratorType, the type of the iterator -/// over the cell centroids, and the ValueTpe, the type of the cell centroid. -/// \tpatam G The type of the grid. -template -struct CellCentroidTraits -{ -}; - -template<> -struct CellCentroidTraits -{ - typedef const double* IteratorType; - typedef const double* ValueType; -}; - -/// \brief Get an iterator over the cell centroids positioned at the first cell. -/// -/// The return type needs to be usable with the functions increment, and -/// getCoordinate. -CellCentroidTraits::IteratorType -beginCellCentroids(const UnstructuredGrid& grid); - - -/// \brief Get vertical position of cell center ("zcorn" average.) -/// \brief grid The grid. -/// \brief cell_index The index of the specific cell. -double cellCenterDepth(const UnstructuredGrid& grid, int cell_index); - - -/// \brief Get a coordinate of a specific cell centroid. -/// \brief grid The grid. -/// \brief cell_index The index of the specific cell. -/// \breif coordinate The coordinate index. -double cellCentroidCoordinate(const UnstructuredGrid& grid, int cell_index, - int coordinate); - - -/// \brief Get the centroid of a cell. -/// \param grid The grid whose cell centroid we query. -/// \param cell_index The index of the corresponding cell. -const double* cellCentroid(const UnstructuredGrid& grid, int cell_index); - - -/// \brief Get the volume of a cell. -/// \param grid The grid the cell belongs to. -/// \param cell_index The index of the cell. -double cellVolume(const UnstructuredGrid& grid, int cell_index); - -/// \brief The mapping of the grid type to type of the iterator over -/// the cell volumes. -/// -/// The value of the mapping is stored in nested type IteratorType -/// \tparam T The type of the grid. -template -struct CellVolumeIteratorTraits -{ -}; - -template<> -struct CellVolumeIteratorTraits -{ - typedef const double* IteratorType; -}; - - /** - Will create an EclipseGrid representation (i.e. based on - ZCORN and COORD) of the current UnstructuredGrid - instance. When creating the UnstructuredGrid the detailed - cornerpoint information is discarded, and it is difficult - to go backwards to recreated ZCORN and COORD. - - The current implementation is based on retaining a copy of the - zcorn keyword after the Minpvprocessor has modified it. - - We then create a new EclipseGrid instance based on the original - input grid, but we "replace" the ZCORN and ACTNUM keywords with - the updated versions. - - If the tolerance in the call to create_grid_cornerpoint( ) is - finite the grid processing code might collapse cells, the z - coordinate transformations from this process will *not* be - correctly represented in the EclipseGrid created by this - method. - */ - -/// \brief Construct an EclipseGrid instance based on the inputGrid, with modifications to -/// zcorn and actnum from the dune UnstructuredGrid. -Opm::EclipseGrid createEclipseGrid(const UnstructuredGrid& grid, const Opm::EclipseGrid& inputGrid ); - -/// \brief Get an iterator over the cell volumes of a grid positioned at the first cell. -const double* beginCellVolumes(const UnstructuredGrid& grid); - -/// \brief Get an iterator over the cell volumes of a grid positioned after the last cell. -const double* endCellVolumes(const UnstructuredGrid& grid); - - - -/// \brief Traits of the face centroids of a grid. -/// -/// This class exports two types: IteratorType, the type of the iterator -/// over the face centroids, and the ValueTpe, the type of the face centroid. -/// \tpatam G The type of the grid. -template -struct FaceCentroidTraits -{ -}; - -template<> -struct FaceCentroidTraits -{ - typedef const double* IteratorType; - typedef const double* ValueType; -}; - -/// \brief Get an iterator over the face centroids positioned at the first cell. -FaceCentroidTraits::IteratorType -beginFaceCentroids(const UnstructuredGrid& grid); - -/// \brief Get a coordinate of a specific face centroid. -/// \param grid The grid. -/// \param face_index The index of the specific face. -FaceCentroidTraits::ValueType -faceCentroid(const UnstructuredGrid& grid, int face_index); - -/// \brief Get the normal of a face. -/// \param grid The grid that the face is part of. -/// \param face_index The index of the face in the grid. -const double* faceNormal(const UnstructuredGrid& grid, int face_index); - -/// \brief Get the area of a face -/// \param grid The grid that the face is part of. -/// \param face_index The index of the face in the grid. -double faceArea(const UnstructuredGrid& grid, int face_index); - -/// \brief Get Eclipse Cartesian tag of a face -/// \param grid The grid that the face is part of. -/// \param cell_face The face attached to a cell as obtained from cell2Faces() -/// \return 0, 1, 2, 3, 4, 5 for I-, I+, J-, J+, K-, K+ -int faceTag(const UnstructuredGrid& grid, boost::iterator_range::const_iterator cell_face); - -/// \brief Maps the grid type to the associated type of the cell to faces mapping. -/// -/// Provides a type named Type. -/// \tparam T The type of the grid. -template -struct Cell2FacesTraits -{ -}; - -template<> -struct Cell2FacesTraits -{ - typedef SparseTableView Type; -}; - -/// \brief Maps the grid type to the associated type of the face to vertices mapping. -/// -/// Provides a type named Type. -/// \tparam T The type of the grid. -template -struct Face2VerticesTraits -{ -}; - -template<> -struct Face2VerticesTraits -{ - typedef SparseTableView Type; -}; - -/// \brief Get the cell to faces mapping of a grid. -Cell2FacesTraits::Type -cell2Faces(const UnstructuredGrid& grid); - -/// \brief Get the face to vertices mapping of a grid. -Face2VerticesTraits::Type -face2Vertices(const UnstructuredGrid& grid); - -/// \brief Get the coordinates of a vertex of the grid. -/// \param grid The grid the vertex is part of. -/// \param index The index identifying the vertex. -const double* vertexCoordinates(const UnstructuredGrid& grid, int index); - -class FaceCellsProxy -{ -public: - FaceCellsProxy(const UnstructuredGrid& grid) - : face_cells_(grid.face_cells) - {} - int operator()(int face_index, int local_index) const - { - return face_cells_[2*face_index+local_index]; - } -private: - const int* face_cells_; -}; - -/// \brief Traits of the face to attached cell mappping of a grid. -/// -/// Exports the type Type, the type of the mapping -/// \tparam T The type of the grid -template -struct FaceCellTraits -{}; - -template<> -struct FaceCellTraits -{ - typedef FaceCellsProxy Type; -}; - -/// \brief Get the face to cell mapping of a grid. -FaceCellTraits::Type faceCells(const UnstructuredGrid& grid); - -/// \brief Increment an iterator over an array that reresents a dense row-major -/// matrix with dims columns -/// \param cc The iterator. -/// \param i The nzumber of rows to increment -/// \param dim The number of columns of the matrix. -template -T* increment(T* cc, int i, int dim) -{ - return cc+(i*dim); -} -/// \brief Increment an iterator over an array that reresents a dense row-major -/// matrix with dims columns -/// \param cc The iterator. -/// \param i The nzumber of rows to increment -template -T increment(const T& t, int i, int) -{ - return t+i; -} - -/// \brief Get the i-th corrdinate of a centroid. -/// \param cc The array with the coordinates. -/// \param i The index of the coordinate. -/// \tparam T The type of the coordinate of the centroid. -template -double getCoordinate(T* cc, int i) -{ - return cc[i]; -} - -/// \brief Get the i-th corrdinate of an array. -/// \param t The iterator over the centroids -/// \brief i The index of the coordinate. -/// \tparam T The type of the iterator representing the centroid. -/// Its value_type has to provide an operator[] to access the coordinates. -template -double getCoordinate(T t, int i) -{ - return (*t)[i]; -} - -} // end namespace UGGridHelpers -} // end namespace OPM -#endif diff --git a/opm/core/grid/GridManager.cpp b/opm/core/grid/GridManager.cpp deleted file mode 100644 index 02019b0d..00000000 --- a/opm/core/grid/GridManager.cpp +++ /dev/null @@ -1,233 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -namespace Opm -{ - - /// Construct a 3d corner-point grid from a deck. - GridManager::GridManager(const Opm::EclipseGrid& inputGrid) - : ug_(0) - { - initFromEclipseGrid(inputGrid, std::vector()); - } - - - GridManager::GridManager(const Opm::EclipseGrid& inputGrid, - const std::vector& poreVolumes) - : ug_(0) - { - initFromEclipseGrid(inputGrid, poreVolumes); - } - - - /// Construct a 2d cartesian grid with cells of unit size. - GridManager::GridManager(int nx, int ny) - { - ug_ = create_grid_cart2d(nx, ny, 1.0, 1.0); - if (!ug_) { - OPM_THROW(std::runtime_error, "Failed to construct grid."); - } - } - - GridManager::GridManager(int nx, int ny,double dx, double dy) - { - ug_ = create_grid_cart2d(nx, ny, dx, dy); - if (!ug_) { - OPM_THROW(std::runtime_error, "Failed to construct grid."); - } - } - - - /// Construct a 3d cartesian grid with cells of unit size. - GridManager::GridManager(int nx, int ny, int nz) - { - ug_ = create_grid_cart3d(nx, ny, nz); - if (!ug_) { - OPM_THROW(std::runtime_error, "Failed to construct grid."); - } - } - - - - - /// Construct a 3d cartesian grid with cells of size [dx, dy, dz]. - GridManager::GridManager(int nx, int ny, int nz, - double dx, double dy, double dz) - { - ug_ = create_grid_hexa3d(nx, ny, nz, dx, dy, dz); - if (!ug_) { - OPM_THROW(std::runtime_error, "Failed to construct grid."); - } - } - - - - - /// Construct a grid from an input file. - /// The file format used is currently undocumented, - /// and is therefore only suited for internal use. - GridManager::GridManager(const std::string& input_filename) - { - ug_ = read_grid(input_filename.c_str()); - if (!ug_) { - OPM_THROW(std::runtime_error, "Failed to read grid from file " << input_filename); - } - } - - /// Destructor. - GridManager::~GridManager() - { - destroy_grid(ug_); - } - - - - - /// Access the managed UnstructuredGrid. - /// The method is named similarly to c_str() in std::string, - /// to make it clear that we are returning a C-compatible struct. - const UnstructuredGrid* GridManager::c_grid() const - { - return ug_; - } - - - - - // Construct corner-point grid from EclipseGrid. - void GridManager::initFromEclipseGrid(const Opm::EclipseGrid& inputGrid, - const std::vector& poreVolumes) - { - struct grdecl g; - int cells_modified = 0; - std::vector actnum; - std::vector coord; - std::vector zcorn; - std::vector mapaxes; - - g.dims[0] = inputGrid.getNX(); - g.dims[1] = inputGrid.getNY(); - g.dims[2] = inputGrid.getNZ(); - - inputGrid.exportMAPAXES( mapaxes ); - inputGrid.exportCOORD( coord ); - inputGrid.exportZCORN( zcorn ); - inputGrid.exportACTNUM( actnum ); - - g.coord = coord.data(); - g.zcorn = zcorn.data(); - g.actnum = actnum.data(); - g.mapaxes = mapaxes.data(); - - if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != MinpvMode::ModeEnum::Inactive)) { - MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]); - const double minpv_value = inputGrid.getMinpvValue(); - // Currently the pinchProcessor is not used and only opmfil is supported - //bool opmfil = inputGrid.getMinpvMode() == MinpvMode::OpmFIL; - bool opmfil = true; - cells_modified = mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data()); - } - - const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0; - ug_ = create_grid_cornerpoint(&g, z_tolerance); - if (!ug_) { - OPM_THROW(std::runtime_error, "Failed to construct grid."); - } - - if (cells_modified > 0) { - attach_zcorn_copy( ug_ , zcorn.data() ); - } - } - - - void GridManager::createGrdecl(const Opm::Deck& deck, struct grdecl &grdecl) - { - // Extract data from deck. - const std::vector& zcorn = deck.getKeyword("ZCORN").getSIDoubleData(); - const std::vector& coord = deck.getKeyword("COORD").getSIDoubleData(); - const int* actnum = NULL; - if (deck.hasKeyword("ACTNUM")) { - actnum = &(deck.getKeyword("ACTNUM").getIntData()[0]); - } - - std::array dims; - if (deck.hasKeyword("DIMENS")) { - const auto& dimensKeyword = deck.getKeyword("DIMENS"); - dims[0] = dimensKeyword.getRecord(0).getItem(0).get< int >(0); - dims[1] = dimensKeyword.getRecord(0).getItem(1).get< int >(0); - dims[2] = dimensKeyword.getRecord(0).getItem(2).get< int >(0); - } else if (deck.hasKeyword("SPECGRID")) { - const auto& specgridKeyword = deck.getKeyword("SPECGRID"); - dims[0] = specgridKeyword.getRecord(0).getItem(0).get< int >(0); - dims[1] = specgridKeyword.getRecord(0).getItem(1).get< int >(0); - dims[2] = specgridKeyword.getRecord(0).getItem(2).get< int >(0); - } else { - OPM_THROW(std::runtime_error, "Deck must have either DIMENS or SPECGRID."); - } - - // Collect in input struct for preprocessing. - - grdecl.zcorn = &zcorn[0]; - grdecl.coord = &coord[0]; - grdecl.actnum = actnum; - grdecl.dims[0] = dims[0]; - grdecl.dims[1] = dims[1]; - grdecl.dims[2] = dims[2]; - - if (deck.hasKeyword("MAPAXES")) { - const auto& mapaxesKeyword = deck.getKeyword("MAPAXES"); - const auto& mapaxesRecord = mapaxesKeyword.getRecord(0); - - // memleak alert: here we need to make sure that C code - // can properly take ownership of the grdecl.mapaxes - // object. if it is not freed, it will result in a - // memleak... - double *cWtfMapaxes = static_cast(malloc(sizeof(double)*mapaxesRecord.size())); - for (unsigned i = 0; i < mapaxesRecord.size(); ++i) - cWtfMapaxes[i] = mapaxesRecord.getItem(i).getSIDouble(0); - grdecl.mapaxes = cWtfMapaxes; - } else - grdecl.mapaxes = NULL; - - } - - - - - - -} // namespace Opm diff --git a/opm/core/grid/GridManager.hpp b/opm/core/grid/GridManager.hpp deleted file mode 100644 index 41ae6455..00000000 --- a/opm/core/grid/GridManager.hpp +++ /dev/null @@ -1,98 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_GRIDMANAGER_HEADER_INCLUDED -#define OPM_GRIDMANAGER_HEADER_INCLUDED - -#include -#include - -#include - -struct UnstructuredGrid; -struct grdecl; - -namespace Opm -{ - /// This class manages an Opm::UnstructuredGrid in the sense that it - /// encapsulates creation and destruction of the grid. - /// The following grid types can be constructed: - /// - 3d corner-point grids (from deck input) - /// - 3d tensor grids (from deck input) - /// - 2d cartesian grids - /// - 3d cartesian grids - /// The resulting UnstructuredGrid is available through the c_grid() method. - class GridManager - { - public: - /// Construct a grid from an EclipseState::EclipseGrid instance. - explicit GridManager(const Opm::EclipseGrid& inputGrid); - - /// Construct a grid from an EclipseState::EclipseGrid instance, - /// giving an explicit set of pore volumes to be used for MINPV - /// considerations. - /// \input[in] eclipseGrid encapsulates a corner-point grid given from a deck - /// \input[in] poreVolumes one element per logical cartesian grid element - GridManager(const Opm::EclipseGrid& inputGrid, - const std::vector& poreVolumes); - - /// Construct a 2d cartesian grid with cells of unit size. - GridManager(int nx, int ny); - - /// Construct a 2d cartesian grid with cells of size [dx, dy]. - GridManager(int nx, int ny, double dx, double dy); - - /// Construct a 3d cartesian grid with cells of unit size. - GridManager(int nx, int ny, int nz); - - /// Construct a 3d cartesian grid with cells of size [dx, dy, dz]. - GridManager(int nx, int ny, int nz, - double dx, double dy, double dz); - - /// Construct a grid from an input file. - /// The file format used is currently undocumented, - /// and is therefore only suited for internal use. - explicit GridManager(const std::string& input_filename); - - /// Destructor. - ~GridManager(); - - /// Access the managed UnstructuredGrid. - /// The method is named similarly to c_str() in std::string, - /// to make it clear that we are returning a C-compatible struct. - const UnstructuredGrid* c_grid() const; - - static void createGrdecl(const Opm::Deck& deck, struct grdecl &grdecl); - - private: - // Disable copying and assignment. - GridManager(const GridManager& other); - GridManager& operator=(const GridManager& other); - - // Construct corner-point grid from EclipseGrid. - void initFromEclipseGrid(const Opm::EclipseGrid& inputGrid, - const std::vector& poreVolumes); - - // The managed UnstructuredGrid. - UnstructuredGrid* ug_; - }; - -} // namespace Opm - -#endif // OPM_GRIDMANAGER_HEADER_INCLUDED diff --git a/opm/core/grid/GridUtilities.cpp b/opm/core/grid/GridUtilities.cpp deleted file mode 100644 index b653c3a2..00000000 --- a/opm/core/grid/GridUtilities.cpp +++ /dev/null @@ -1,133 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -namespace Opm -{ - /// For each cell, find indices of all other cells sharing a vertex with it. - /// \param[in] grid A grid object. - /// \return A table of neighbour cell-indices by cell. - SparseTable cellNeighboursAcrossVertices(const UnstructuredGrid& grid) - { - // 1. Create vertex->cell mapping. We do this by iterating - // over all faces, and adding both its cell neighbours - // to each of its vertices' data. - using namespace UgGridHelpers; - const int num_vertices = grid.number_of_nodes; - std::vector> v2c(num_vertices); - const int num_faces = numFaces(grid); - const auto fc = faceCells(grid); - for (int face = 0; face < num_faces; ++face) { - for (int nodepos = grid.face_nodepos[face]; nodepos < grid.face_nodepos[face + 1]; ++nodepos) { - const int vertex = grid.face_nodes[nodepos]; - for (int face_nb = 0; face_nb < 2; ++face_nb) { - const int face_nb_cell = fc(face, face_nb); - if (face_nb_cell >= 0) { - v2c[vertex].insert(face_nb_cell); - } - } - } - } - - // 2. For each cell, iterate over its faces, iterate over - // their vertices, and collect all those vertices' cell - // neighbours. Add as row to sparse table. - SparseTable cell_nb; - const int num_cells = numCells(grid); - const auto c2f = cell2Faces(grid); - // Reserve sufficient room for cartesian grids in 2 and 3 - // dimensions. Note that this is not a limit, just an - // optimization similar to std::vector. - cell_nb.reserve(num_cells, (dimensions(grid) == 2 ? 8 : 26) * num_cells); - std::set nb; - for (int cell = 0; cell < num_cells; ++cell) { - nb.clear(); - const auto cell_faces = c2f[cell]; - const int num_cell_faces = cell_faces.size(); - for (int local_face = 0; local_face < num_cell_faces; ++local_face) { - const int face = cell_faces[local_face]; - for (int nodepos = grid.face_nodepos[face]; nodepos < grid.face_nodepos[face + 1]; ++nodepos) { - const int vertex = grid.face_nodes[nodepos]; - nb.insert(v2c[vertex].begin(), v2c[vertex].end()); - } - } - nb.erase(cell); - cell_nb.appendRow(nb.begin(), nb.end()); - } - - // 3. Done. Return. - return cell_nb; - } - - - - - - - /// For each cell, order the (cell) neighbours counterclockwise. - /// \param[in] grid A 2d grid object. - /// \param[in, out] nb A cell-cell neighbourhood table, such as from cellNeighboursAcrossVertices(). - void orderCounterClockwise(const UnstructuredGrid& grid, - SparseTable& nb) - { - if (grid.dimensions != 2) { - OPM_THROW(std::logic_error, "Cannot use orderCounterClockwise in " << grid.dimensions << " dimensions."); - } - const int num_cells = grid.number_of_cells; - if (nb.size() != num_cells) { - OPM_THROW(std::logic_error, "Inconsistent arguments for orderCounterClockwise()."); - } - - // For each cell, compute each neighbour's angle with the x axis, - // sort that to find the correct permutation of the neighbours. - typedef std::pair AngleAndPos; - std::vector angle_and_pos; - std::vector original; - for (int cell = 0; cell < num_cells; ++cell) { - const int num_nb = nb[cell].size(); - angle_and_pos.clear(); - angle_and_pos.resize(num_nb); - for (int ii = 0; ii < num_nb; ++ii) { - const int cell2 = nb[cell][ii]; - const double v[2] = { grid.cell_centroids[2*cell2] - grid.cell_centroids[2*cell], - grid.cell_centroids[2*cell2 + 1] - grid.cell_centroids[2*cell + 1] }; - // The formula below gives an angle in [0, 2*pi] with the positive x axis. - const double angle = boost::math::constants::pi() - std::atan2(v[1], -v[0]); - angle_and_pos[ii] = std::make_pair(angle, ii); - } - original.assign(nb[cell].begin(), nb[cell].end()); - std::sort(angle_and_pos.begin(), angle_and_pos.end()); - for (int ii = 0; ii < num_nb; ++ii) { - nb[cell][ii] = original[angle_and_pos[ii].second]; - } - } - } - -} // namespace Opm diff --git a/opm/core/grid/GridUtilities.hpp b/opm/core/grid/GridUtilities.hpp deleted file mode 100644 index 78363b40..00000000 --- a/opm/core/grid/GridUtilities.hpp +++ /dev/null @@ -1,42 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_GRIDUTILITIES_HEADER_INCLUDED -#define OPM_GRIDUTILITIES_HEADER_INCLUDED - -#include -#include - -namespace Opm -{ - - /// For each cell, find indices of all cells sharing a vertex with it. - /// \param[in] grid A grid object. - /// \return A table of neighbour cell-indices by cell. - SparseTable cellNeighboursAcrossVertices(const UnstructuredGrid& grid); - - /// For each cell, order the (cell) neighbours counterclockwise. - /// \param[in] grid A 2d grid object. - /// \param[in, out] nb A cell-cell neighbourhood table, such as from vertexNeighbours(). - void orderCounterClockwise(const UnstructuredGrid& grid, - SparseTable& nb); - -} // namespace Opm - -#endif // OPM_GRIDUTILITIES_HEADER_INCLUDED diff --git a/opm/core/grid/MinpvProcessor.hpp b/opm/core/grid/MinpvProcessor.hpp deleted file mode 100644 index 64513489..00000000 --- a/opm/core/grid/MinpvProcessor.hpp +++ /dev/null @@ -1,176 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_MINPVPROCESSOR_HEADER_INCLUDED -#define OPM_MINPVPROCESSOR_HEADER_INCLUDED - - -#include -#include - - -namespace Opm -{ - - /// \brief Transform a corner-point grid ZCORN field to account for MINPV processing. - class MinpvProcessor - { - public: - /// \brief Create a processor. - /// \param[in] nx logical cartesian number of cells in I-direction - /// \param[in] ny logical cartesian number of cells in J-direction - /// \param[in] nz logical cartesian number of cells in K-direction - MinpvProcessor(const int nx, const int ny, const int nz); - /// Change zcorn so that it respects the minpv property. - /// \param[in] pv pore volumes of all logical cartesian cells - /// \param[in] minpv minimum pore volume to accept a cell - /// \param[in] actnum active cells, inactive cells are not considered - /// \param[in] mergeMinPVCells flag to determine whether cells below minpv - /// should be included in the cell below - /// \param[in, out] zcorn ZCORN array to be manipulated - /// After processing, all cells that have lower pore volume than minpv - /// will have the zcorn numbers changed so they are zero-thickness. Any - /// cell below will be changed to include the deleted volume if mergeMinPCCells is true - /// els the volume will be lost - int process(const std::vector& pv, - const double minpv, - const std::vector& actnum, - const bool mergeMinPVCells, - double* zcorn) const; - private: - std::array cornerIndices(const int i, const int j, const int k) const; - std::array getCellZcorn(const int i, const int j, const int k, const double* z) const; - void setCellZcorn(const int i, const int j, const int k, const std::array& cellz, double* z) const; - std::array dims_; - std::array delta_; - }; - - inline MinpvProcessor::MinpvProcessor(const int nx, const int ny, const int nz) : - dims_( {nx,ny,nz} ), - delta_( {1 , 2*nx , 4*nx*ny} ) - { } - - - - inline int MinpvProcessor::process(const std::vector& pv, - const double minpv, - const std::vector& actnum, - const bool mergeMinPVCells, - double* zcorn) const - { - // Algorithm: - // 1. Process each column of cells (with same i and j - // coordinates) from top (low k) to bottom (high k). - // 2. For each cell 'c' visited, check if its pore volume - // pv[c] is less than minpv. - // 3. If below the minpv threshold, move the lower four - // zcorn associated with the cell c to coincide with - // the upper four (so it becomes degenerate). If mergeMinPVcells - // is true, the higher four zcorn associated with the cell below - // is moved to these values (so it gains the deleted volume). - - int cells_modified = 0; - - // Check for sane input sizes. - const size_t log_size = dims_[0] * dims_[1] * dims_[2]; - if (pv.size() != log_size) { - OPM_THROW(std::runtime_error, "Wrong size of PORV input, must have one element per logical cartesian cell."); - } - if (!actnum.empty() && actnum.size() != log_size) { - OPM_THROW(std::runtime_error, "Wrong size of ACTNUM input, must have one element per logical cartesian cell."); - } - - // Main loop. - for (int kk = 0; kk < dims_[2]; ++kk) { - for (int jj = 0; jj < dims_[1]; ++jj) { - for (int ii = 0; ii < dims_[0]; ++ii) { - const int c = ii + dims_[0] * (jj + dims_[1] * kk); - if (pv[c] < minpv && (actnum.empty() || actnum[c])) { - // Move deeper (higher k) coordinates to lower k coordinates. - std::array cz = getCellZcorn(ii, jj, kk, zcorn); - for (int count = 0; count < 4; ++count) { - cz[count + 4] = cz[count]; - } - setCellZcorn(ii, jj, kk, cz, zcorn); - - - // optionally add removed volume to the cell below. - if (mergeMinPVCells) { - // Check if there is a cell below. - if (pv[c] > 0.0 && kk < dims_[2] - 1) { - // Set lower k coordinates of cell below to upper cells's coordinates. - std::array cz_below = getCellZcorn(ii, jj, kk + 1, zcorn); - for (int count = 0; count < 4; ++count) { - cz_below[count] = cz[count]; - } - setCellZcorn(ii, jj, kk + 1, cz_below, zcorn); - } - } - ++cells_modified; - } - } - } - } - return cells_modified; - } - - - - inline std::array MinpvProcessor::cornerIndices(const int i, const int j, const int k) const - { - const int ix = 2*(i*delta_[0] + j*delta_[1] + k*delta_[2]); - std::array ixs = {{ ix, ix + delta_[0], - ix + delta_[1], ix + delta_[1] + delta_[0], - ix + delta_[2], ix + delta_[2] + delta_[0], - ix + delta_[2] + delta_[1], ix + delta_[2] + delta_[1] + delta_[0] }}; - - return ixs; - } - - - - // Returns the eight z-values associated with a given cell. - // The ordering is such that i runs fastest. That is, with - // L = low and H = high: - // {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }. - inline std::array MinpvProcessor::getCellZcorn(const int i, const int j, const int k, const double* z) const - { - const std::array ixs = cornerIndices(i, j, k); - std::array cellz; - for (int count = 0; count < 8; ++count) { - cellz[count] = z[ixs[count]]; - } - return cellz; - } - - - - inline void MinpvProcessor::setCellZcorn(const int i, const int j, const int k, const std::array& cellz, double* z) const - { - const std::array ixs = cornerIndices(i, j, k); - for (int count = 0; count < 8; ++count) { - z[ixs[count]] = cellz[count]; - } - } - - - -} // namespace Opm - -#endif // OPM_MINPVPROCESSOR_HEADER_INCLUDED diff --git a/opm/core/grid/PinchProcessor.hpp b/opm/core/grid/PinchProcessor.hpp deleted file mode 100755 index 25d5252c..00000000 --- a/opm/core/grid/PinchProcessor.hpp +++ /dev/null @@ -1,488 +0,0 @@ -/* - Copyright 2015 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_PINCHPROCESSOR_HEADER_INCLUDED -#define OPM_PINCHPROCESSOR_HEADER_INCLUDED - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace Opm -{ - - template - class PinchProcessor - { - public: - /// \brief Create a Pinch processor. - /// \param[in] minpvValue value in MINPV keyword - /// \param[in] thickness item 2 in PINCH keyword - /// \param[in] transMode item 4 in PINCH keyword - /// \param[in] multzMode item 5 in PINCH keyword - PinchProcessor(const double minpvValue, - const double thickness, - const PinchMode::ModeEnum transMode, - const PinchMode::ModeEnum multzMode); - /// Generate NNCs for cells which pv is less than MINPV. - /// \param[in] Grid cpgrid or unstructured grid - /// \param[in] htrans half cell transmissibility, size is number of cellfaces. - /// \param[in] multz Z+ transmissibility multiplier for all active cells - /// \param[in] pv pore volume for all the cartesian cells - /// \param[in] nnc non-neighbor connection class - /// Algorithm: - /// 1. Mark all the cells which pv less than minpvValue. - /// 2. Find out proper pinchouts column and associate top and bottom cells. - /// 3. Compute transmissibility for nncs. - /// 4. Apply multz due to different multz options. - void process(const Grid& grid, - const std::vector& htrans, - const std::vector& actnum, - const std::vector& multz, - const std::vector& pv, - NNC& nnc); - - private: - double minpvValue_; - double thickness_; - PinchMode::ModeEnum transMode_; - PinchMode::ModeEnum multzMode_; - - /// Mark minpved cells. - std::vector getMinpvCells_(const std::vector& actnum, - const std::vector& pv); - - /// Get the interface for two cells. - int interface_(const Grid& grid, - const int cellIdx1, - const int cellIdx2); - - /// Get the proper face for one cell. - int interface_(const Grid& grid, - const int cellIdx, - const Opm::FaceDir::DirEnum& faceDir); - - /// Get pinchouts column. - std::vector > - getPinchoutsColumn_(const Grid& grid, - const std::vector& actnum, - const std::vector& pv); - - /// Get global cell index. - int getGlobalIndex_(const int i, const int j, const int k, const int* dims); - - /// Get cartesian index. - std::array getCartIndex_(const int idx, - const int* dims); - - /// Compute transmissibility for nnc. - std::vector transCompute_(const Grid& grid, - const std::vector& htrans, - const std::vector& pinCells, - const std::vector& pinFaces); - - /// Get map between half-trans index and the pair of face index and cell index. - std::vector getHfIdxMap_(const Grid& grid); - - /// Get active cell index. - int getActiveCellIdx_(const Grid& grid, - const int globalIdx); - - /// Item 4 in PINCH keyword. - void transTopbot_(const Grid& grid, - const std::vector& htrans, - const std::vector& actnum, - const std::vector& multz, - const std::vector& pv, - NNC& nnc); - - /// Item 5 in PINCH keyword. - std::unordered_multimap multzOptions_(const Grid& grid, - const std::vector& pinCells, - const std::vector& pinFaces, - const std::vector& multz, - const std::vector >& seg); - - /// Apply multz vector to face transmissibility. - void applyMultz_(std::vector& trans, - const std::unordered_multimap& multzmap); - - }; - - - - template - inline PinchProcessor::PinchProcessor(const double minpv, - const double thickness, - const PinchMode::ModeEnum transMode, - const PinchMode::ModeEnum multzMode) - { - minpvValue_ = minpv; - thickness_ = thickness; - transMode_ = transMode; - multzMode_ = multzMode; - } - - - - template - inline int PinchProcessor::getGlobalIndex_(const int i, const int j, const int k, const int* dims) - { - return i + dims[0] * (j + dims[1] * k); - } - - - - template - inline std::array PinchProcessor::getCartIndex_(const int idx, - const int* dims) - { - std::array ijk; - ijk[0] = (idx % dims[0]); - ijk[1] = ((idx / dims[0]) % dims[1]); - ijk[2] = ((idx / dims[0]) / dims[1]); - - return ijk; - } - - - - template - inline int PinchProcessor::interface_(const Grid& grid, - const int cellIdx1, - const int cellIdx2) - { - const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid); - int commonFace = -1; - const int actCellIdx1 = getActiveCellIdx_(grid, cellIdx1); - const int actCellIdx2 = getActiveCellIdx_(grid, cellIdx2); - const auto cellFacesRange1 = cell_faces[actCellIdx1]; - const auto cellFacesRange2 = cell_faces[actCellIdx2]; - for (const auto& f1 : cellFacesRange1) { - for (const auto& f2 : cellFacesRange2) { - if (f1 == f2) { - commonFace = f1; - break; - } - } - } - - if (commonFace == -1) { - const auto dims = Opm::UgGridHelpers::cartDims(grid); - const auto ijk1 = getCartIndex_(cellIdx1, dims); - const auto ijk2 = getCartIndex_(cellIdx2, dims); - - OPM_THROW(std::logic_error, "Couldn't find the common face for cell " - << cellIdx1<< "("< - inline int PinchProcessor::interface_(const Grid& grid, - const int cellIdx, - const Opm::FaceDir::DirEnum& faceDir) - { - const auto actCellIdx = getActiveCellIdx_(grid, cellIdx); - const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid); - const auto cellFacesRange = cell_faces[actCellIdx]; - int faceIdx = -1; - for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) { - int tag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter); - if ( (faceDir == Opm::FaceDir::ZMinus && tag == 4) || (faceDir == Opm::FaceDir::ZPlus && tag == 5) ) { - faceIdx = *cellFaceIter; - } - } - - if (faceIdx == -1) { - OPM_THROW(std::logic_error, "Couldn't find the face for cell ." << cellIdx); - } - - return faceIdx; - } - - - - template - inline std::vector PinchProcessor::getMinpvCells_(const std::vector& actnum, - const std::vector& pv) - { - std::vector minpvCells(pv.size(), 0); - for (int idx = 0; idx < static_cast(pv.size()); ++idx) { - if (actnum[idx]) { - if (pv[idx] < minpvValue_) { - minpvCells[idx] = 1; - } - } - } - return minpvCells; - } - - - - template - inline std::vector PinchProcessor::getHfIdxMap_(const Grid& grid) - { - std::vector hf_ix(2*Opm::UgGridHelpers::numFaces(grid), -1); - const auto& f2c = Opm::UgGridHelpers::faceCells(grid); - const auto& cf = Opm::UgGridHelpers::cell2Faces(grid); - - for (int c = 0, i = 0; c < Opm::UgGridHelpers::numCells(grid); ++c) { - for (const auto& f: cf[c]) { - const auto off = 0 + (f2c(f, 0) != c); - hf_ix[2*f + off] = i++; - } - } - return hf_ix; - } - - - - template - inline int PinchProcessor::getActiveCellIdx_(const Grid& grid, - const int globalIdx) - { - const int nc = Opm::UgGridHelpers::numCells(grid); - const int* global_cell = Opm::UgGridHelpers::globalCell(grid); - int idx = -1; - for (int i = 0; i < nc; ++i) { - if (global_cell[i] == globalIdx) { - idx = i; - break; - } - } - return idx; - } - - - - template - inline std::vector PinchProcessor::transCompute_(const Grid& grid, - const std::vector& htrans, - const std::vector& pinCells, - const std::vector& pinFaces) - { - const int nc = Opm::UgGridHelpers::numCells(grid); - const int nf = Opm::UgGridHelpers::numFaces(grid); - std::vector trans(nf, 0); - int cellFaceIdx = 0; - auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid); - const auto& hfmap = getHfIdxMap_(grid); - const auto& f2c = Opm::UgGridHelpers::faceCells(grid); - for (int cellIdx = 0; cellIdx < nc; ++cellIdx) { - auto cellFacesRange = cell_faces[cellIdx]; - for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter, ++cellFaceIdx) { - const int faceIdx = *cellFaceIter; - const auto pos = std::find(pinFaces.begin(), pinFaces.end(), faceIdx); - if (pos == pinFaces.end()) { - trans[faceIdx] += 1. / htrans[cellFaceIdx]; - } else { - const int idx1 = std::distance(std::begin(pinFaces), pos); - int idx2; - if (idx1 % 2 == 0) { - idx2 = idx1 + 1; - } else { - idx2 = idx1 - 1; - } - const int f1 = hfmap[2*pinFaces[idx1] + (f2c(pinFaces[idx1], 0) != getActiveCellIdx_(grid, pinCells[idx1]))]; - const int f2 = hfmap[2*pinFaces[idx2] + (f2c(pinFaces[idx2], 0) != getActiveCellIdx_(grid, pinCells[idx2]))]; - trans[faceIdx] = (1. / htrans[f1] + 1. / htrans[f2]); - trans[pinFaces[idx2]] = trans[faceIdx]; - } - } - } - - for (auto f = 0; f < nf; ++f) { - trans[f] = 1. / trans[f]; - } - - return trans; - } - - - - template - inline std::vector> PinchProcessor::getPinchoutsColumn_(const Grid& grid, - const std::vector& actnum, - const std::vector& pv) - { - const int* dims = Opm::UgGridHelpers::cartDims(grid); - std::vector minpvCells = getMinpvCells_(actnum, pv); - std::vector> segment; - for (int z = 0; z < dims[2]; ++z) { - for (int y = 0; y < dims[1]; ++y) { - for (int x = 0; x < dims[0]; ++x) { - const int c = getGlobalIndex_(x, y, z, dims); - std::vector seg; - if (minpvCells[c]) { - seg.push_back(c); - minpvCells[c] = 0; - for (int zz = z+1; zz < dims[2]; ++zz) { - const int cc = getGlobalIndex_(x, y, zz, dims); - if (minpvCells[cc]) { - seg.push_back(cc); - minpvCells[cc] = 0; - } else { - break; - } - } - segment.push_back(seg); - } - } - } - } - - return segment; - } - - - - template - inline void PinchProcessor::transTopbot_(const Grid& grid, - const std::vector& htrans, - const std::vector& actnum, - const std::vector& multz, - const std::vector& pv, - NNC& nnc) - { - const int* dims = Opm::UgGridHelpers::cartDims(grid); - std::vector pinFaces; - std::vector pinCells; - std::vector > newSeg; - auto minpvSeg = getPinchoutsColumn_(grid, actnum, pv); - for (auto& seg : minpvSeg) { - std::array ijk1 = getCartIndex_(seg.front(), dims); - std::array ijk2 = getCartIndex_(seg.back(), dims); - auto tmp = seg; - if ((ijk1[2]-1) >= 0 && (ijk2[2]+1) < dims[2]) { - int topCell = getGlobalIndex_(ijk1[0], ijk1[1], ijk1[2]-1, dims); - int botCell = getGlobalIndex_(ijk2[0], ijk2[1], ijk2[2]+1, dims); - /// for any segments, we need to find the active top and bottom cells. - /// if the original segment's top and bottom is inactive, we need to lookup - /// the column until they're found otherwise just ignore this segment. - if (!actnum[topCell]) { - for (int topk = ijk1[2]-2; topk > 0; --topk) { - topCell = getGlobalIndex_(ijk1[0], ijk1[1], topk, dims); - if (actnum[topCell]) { - break; - } - } - } - pinFaces.push_back(interface_(grid, topCell, Opm::FaceDir::ZPlus)); - pinCells.push_back(topCell); - - tmp.insert(tmp.begin(), topCell); - newSeg.push_back(tmp); - if (!actnum[botCell]) { - for (int botk = ijk2[2]+2; botk < dims[2]; ++botk) { - botCell = getGlobalIndex_(ijk2[0], ijk2[1], botk, dims); - if (actnum[botCell]) { - break; - - } - } - } - pinFaces.push_back(interface_(grid, botCell, Opm::FaceDir::ZMinus)); - pinCells.push_back(botCell); - - } - } - - auto faceTrans = transCompute_(grid, htrans, pinCells, pinFaces); - auto multzmap = multzOptions_(grid, pinCells, pinFaces, multz, newSeg); - applyMultz_(faceTrans, multzmap); - for (int i = 0; i < static_cast(pinCells.size())/2; ++i) { - nnc.addNNC(static_cast(pinCells[2*i]), static_cast(pinCells[2*i+1]), faceTrans[pinFaces[2*i]]); - } - } - - - - - template - inline std::unordered_multimap PinchProcessor::multzOptions_(const Grid& grid, - const std::vector& pinCells, - const std::vector& pinFaces, - const std::vector& multz, - const std::vector >& segs) - { - std::unordered_multimap multzmap; - if (multzMode_ == PinchMode::ModeEnum::TOP) { - for (int i = 0; i < static_cast(pinFaces.size())/2; ++i) { - multzmap.insert(std::make_pair(pinFaces[2*i], multz[getActiveCellIdx_(grid, pinCells[2*i])])); - multzmap.insert(std::make_pair(pinFaces[2*i+1],multz[getActiveCellIdx_(grid, pinCells[2*i])])); - } - } else if (multzMode_ == PinchMode::ModeEnum::ALL) { - for (auto& seg : segs) { - //find the min multz in seg cells. - auto multzValue = std::numeric_limits::max(); - for (auto& cellIdx : seg) { - auto activeIdx = getActiveCellIdx_(grid, cellIdx); - if (activeIdx != -1) { - multzValue = std::min(multzValue, multz[activeIdx]); - } - } - //find the right face. - auto index = std::distance(std::begin(pinCells), std::find(pinCells.begin(), pinCells.end(), seg.front())); - multzmap.insert(std::make_pair(pinFaces[index], multzValue)); - multzmap.insert(std::make_pair(pinFaces[index+1], multzValue)); - } - } - - return multzmap; - } - - - - template - inline void PinchProcessor::applyMultz_(std::vector& trans, - const std::unordered_multimap& multzmap) - { - for (auto& x : multzmap) { - trans[x.first] *= x.second; - } - } - - - - template - inline void PinchProcessor::process(const Grid& grid, - const std::vector& htrans, - const std::vector& actnum, - const std::vector& multz, - const std::vector& pv, - NNC& nnc) - { - transTopbot_(grid, htrans, actnum, multz, pv, nnc); - } - -} // namespace Opm -#endif // OPM_PINCHPROCESSOR_HEADER_INCLUDED diff --git a/opm/core/grid/cart_grid.c b/opm/core/grid/cart_grid.c deleted file mode 100644 index c6a0a8cc..00000000 --- a/opm/core/grid/cart_grid.c +++ /dev/null @@ -1,708 +0,0 @@ -/*=========================================================================== -// -// File: cart_grid.c -// -// Author: Jostein R. Natvig -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include - -#include -#include - -#include - -static struct UnstructuredGrid *allocate_cart_grid_3d(int nx, int ny, int nz); -static void fill_cart_topology_3d(struct UnstructuredGrid *G); -static void fill_cart_geometry_3d(struct UnstructuredGrid *G, - const double *x, - const double *y, - const double *z); -static void -fill_layered_geometry_3d(struct UnstructuredGrid *G, - const double *x, - const double *y, - const double *z, - const double *depthz); - -struct UnstructuredGrid * -create_grid_cart3d(int nx, int ny, int nz) -{ - return create_grid_hexa3d(nx, ny, nz, 1.0, 1.0, 1.0); -} - -struct UnstructuredGrid * -create_grid_hexa3d(int nx, int ny, int nz, - double dx, double dy, double dz) -{ - int i; - double *x, *y, *z; - struct UnstructuredGrid *G; - - x = malloc((nx + 1) * sizeof *x); - y = malloc((ny + 1) * sizeof *y); - z = malloc((nz + 1) * sizeof *z); - - if ((x == NULL) || (y == NULL) || (z == NULL)) { - G = NULL; - } else { - for (i = 0; i < nx + 1; i++) { x[i] = i * dx; } - for (i = 0; i < ny + 1; i++) { y[i] = i * dy; } - for (i = 0; i < nz + 1; i++) { z[i] = i * dz; } - - G = create_grid_tensor3d(nx, ny, nz, x, y, z, - (const double *) NULL); - } - - free(z); free(y); free(x); - - return G; -} - -/* --------------------------------------------------------------------- */ - -static struct UnstructuredGrid *allocate_cart_grid_2d(int nx, int ny); -static void fill_cart_topology_2d(struct UnstructuredGrid *G); -static void fill_cart_geometry_2d(struct UnstructuredGrid *G, - const double *x, - const double *y); - -struct UnstructuredGrid* -create_grid_cart2d(int nx, int ny, double dx, double dy) -{ - int i; - double *x, *y; - struct UnstructuredGrid *G; - - x = malloc((nx + 1) * sizeof *x); - y = malloc((ny + 1) * sizeof *y); - - if ((x == NULL) || (y == NULL)) { - G = NULL; - } else { - - for (i = 0; i < nx + 1; i++) { x[i] = i*dx; } - for (i = 0; i < ny + 1; i++) { y[i] = i*dy; } - - G = create_grid_tensor2d(nx, ny, x, y); - } - - free(y); free(x); - - return G; -} - -/* --------------------------------------------------------------------- */ - -struct UnstructuredGrid * -create_grid_tensor2d(int nx, int ny, const double *x, const double *y) -{ - struct UnstructuredGrid *G; - - G = allocate_cart_grid_2d(nx, ny); - - if (G != NULL) - { - fill_cart_topology_2d(G); - fill_cart_geometry_2d(G, x, y); - } - - return G; -} - -/* --------------------------------------------------------------------- */ - -struct UnstructuredGrid * -create_grid_tensor3d(int nx , - int ny , - int nz , - const double *x , - const double *y , - const double *z , - const double *depthz) -{ - struct UnstructuredGrid *G; - - G = allocate_cart_grid_3d(nx, ny, nz); - - if (G != NULL) - { - fill_cart_topology_3d(G); - - if (depthz == NULL) { - fill_cart_geometry_3d(G, x, y, z); - } - else { - fill_layered_geometry_3d(G, x, y, z, depthz); - } - } - - return G; -} - -/* --------------------------------------------------------------------- */ -/* Static functions follow: */ -/* --------------------------------------------------------------------- */ - - -static struct UnstructuredGrid * -allocate_cart_grid(size_t ndims , - size_t ncells, - size_t nfaces, - size_t nnodes) -{ - size_t nfacenodes, ncellfaces; - - nfacenodes = nfaces * (2 * (ndims - 1)); - ncellfaces = ncells * (2 * ndims); - - return allocate_grid(ndims, ncells, nfaces, - nfacenodes, ncellfaces, nnodes); -} - - -static struct UnstructuredGrid* -allocate_cart_grid_3d(int nx, int ny, int nz) -{ - struct UnstructuredGrid *G; - int Nx, Ny, Nz; - int nxf, nyf, nzf; - - int ncells, nfaces, nnodes; - - Nx = nx + 1; - Ny = ny + 1; - Nz = nz +1; - - nxf = Nx * ny * nz; - nyf = nx * Ny * nz; - nzf = nx * ny * Nz; - - ncells = nx * ny * nz ; - nfaces = nxf + nyf + nzf; - nnodes = Nx * Ny * Nz ; - - G = allocate_cart_grid(3, ncells, nfaces, nnodes); - - if (G != NULL) - { - G->dimensions = 3 ; - G->cartdims[0] = nx; - G->cartdims[1] = ny; - G->cartdims[2] = nz; - - G->number_of_cells = ncells; - G->number_of_faces = nfaces; - G->number_of_nodes = nnodes; - } - - return G; -} - -/* --------------------------------------------------------------------- */ - -static void -fill_cart_topology_3d(struct UnstructuredGrid *G) -{ - int nx, ny, nz; - int Nx, Ny; - int nxf, nyf; - int i,j,k; - - int *cfaces, *cfacepos, *fnodes, *fnodepos, *fcells; - - nx = G->cartdims[0]; - ny = G->cartdims[1]; - nz = G->cartdims[2]; - - Nx = nx+1; - Ny = ny+1; - - nxf = Nx*ny*nz; - nyf = nx*Ny*nz; - - cfaces = G->cell_faces; - cfacepos = G->cell_facepos; - cfacepos[0] = 0; - for (k=0; kcell_facetag[k*6 + i] = i; - } - } - - fnodes = G->face_nodes; - fnodepos = G->face_nodepos; - fcells = G->face_cells; - fnodepos[0] = 0; - - /* Faces with x-normal */ - for (k=0; kcartdims[0]; - ny = G->cartdims[1]; - nz = G->cartdims[2]; - - ccentroids = G->cell_centroids; - cvolumes = G->cell_volumes; - for (k=0; kface_normals; - fcentroids = G->face_centroids; - fareas = G->face_areas; - - /* Faces with x-normal */ - for (k=0; knode_coordinates; - for (k=0; kcartdims[0]; ny = G->cartdims[1]; nz = G->cartdims[2]; - - coord = G->node_coordinates; - for (k = 0; k < nz + 1; k++) { - - depth = depthz; - - for (j = 0; j < ny + 1; j++) { - for (i = 0; i < nx + 1; i++) { - *coord++ = x[i]; - *coord++ = y[j]; - *coord++ = z[k] + *depth++; - } - } - } - - compute_geometry(G); -} - -/* --------------------------------------------------------------------- */ - -static struct UnstructuredGrid* -allocate_cart_grid_2d(int nx, int ny) -{ - int nxf, nyf; - int Nx , Ny ; - - int ncells, nfaces, nnodes; - - struct UnstructuredGrid *G; - - Nx = nx + 1; - Ny = ny + 1; - - nxf = Nx * ny; - nyf = nx * Ny; - - ncells = nx * ny ; - nfaces = nxf + nyf; - nnodes = Nx * Ny ; - - G = allocate_cart_grid(2, ncells, nfaces, nnodes); - - if (G != NULL) - { - G->dimensions = 2 ; - G->cartdims[0] = nx; - G->cartdims[1] = ny; - G->cartdims[2] = 1 ; - - G->number_of_cells = ncells; - G->number_of_faces = nfaces; - G->number_of_nodes = nnodes; - } - - return G; -} - -/* --------------------------------------------------------------------- */ - -static void -fill_cart_topology_2d(struct UnstructuredGrid *G) -{ - int i,j; - int nx, ny; - int nxf; - int Nx; - - int *fnodes, *fnodepos, *fcells, *cfaces, *cfacepos; - - cfaces = G->cell_faces; - cfacepos = G->cell_facepos; - - nx = G->cartdims[0]; - ny = G->cartdims[1]; - Nx = nx + 1; - nxf = Nx * ny; - - cfacepos[0] = 0; - for (j=0; jcell_facetag[j*4 + 0] = 0; - G->cell_facetag[j*4 + 1] = 2; - G->cell_facetag[j*4 + 2] = 1; - G->cell_facetag[j*4 + 3] = 3; - } - - - fnodes = G->face_nodes; - fnodepos = G->face_nodepos; - fcells = G->face_cells; - fnodepos[0] = 0; - - /* Faces with x-normal */ - for (j=0; jcartdims[0]; - ny = G->cartdims[1]; - - ccentroids = G->cell_centroids; - cvolumes = G->cell_volumes; - - for (j=0; jface_normals; - fcentroids = G->face_centroids; - fareas = G->face_areas; - - /* Faces with x-normal */ - for (j=0; jnode_coordinates; - for (j=0; j -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. - Copyright 2011 Statoil ASA. - - This file is part of the Open Porous Media Project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_CART_GRID_H_HEADER -#define OPM_CART_GRID_H_HEADER - -/** - * \file - * Routines to construct fully formed grid structures from a simple Cartesian - * (i.e., tensor product) description. - * - * The cells are lexicographically ordered with the @c i index cycling the most - * rapidly, followed by the @c j index and then, in three space dimensions, the - * @c k (`layer') index as the least rapidly cycling index. - */ - -#ifdef __cplusplus -extern "C" { -#endif - - - -/** - * Form geometrically Cartesian grid in two space dimensions with equally - * sized cells. - * - * @param[in] nx Number of cells in @c x direction. - * @param[in] ny Number of cells in @c y direction. - * @param[in] dx Length, in meters, of each cell's @c x extent. - * @param[in] dy Length, in meters, of each cell's @c y extent. - * - * @return Fully formed grid structure containing valid geometric primitives. - * Must be destroyed using function destroy_grid(). - */ -struct UnstructuredGrid * -create_grid_cart2d(int nx, int ny, double dx, double dy); - - -/** - * Form geometrically Cartesian grid in three space dimensions with unit-sized - * cells. - * - * @param[in] nx Number of cells in @c x direction. - * @param[in] ny Number of cells in @c y direction. - * @param[in] nz Number of cells in @c z direction. - * - * @return Fully formed grid structure containing valid geometric primitives. - * Must be destroyed using function destroy_grid(). - */ -struct UnstructuredGrid * -create_grid_cart3d(int nx, int ny, int nz); - - -/** - * Form geometrically Cartesian grid in three space dimensions with equally - * sized cells. - * - * Each cell has physical size (volume) \f$\mathit{dx}\times \mathit{dy}\times - * \mathit{dz}\f$. - * - * @param[in] nx Number of cells in @c x direction. - * @param[in] ny Number of cells in @c y direction. - * @param[in] nz Number of cells in @c z direction. - * - * @param[in] dx Length, in meters, of each cell's @c x extent. - * @param[in] dy Length, in meters, of each cell's @c y extent. - * @param[in] dz Length, in meters, of each cell's @c z extent. - * - * @return Fully formed grid structure containing valid geometric primitives. - * Must be destroyed using function destroy_grid(). - */ -struct UnstructuredGrid * -create_grid_hexa3d(int nx, int ny, int nz, - double dx, double dy, double dz); - - -/** - * Form tensor product (Cartesian) grid in two space dimensions. - * - * The size (volume) of cell \f$(i,j)\f$ is - * \f[ - * v_{ij} = (x_{i+1} - x_i)\cdot (y_{j+1} - y_j) - * \f] - * Similar relations hold for the cell and interface centroids as well as the - * interface areas and normal vectors. In other words, cell \f$(i,j)\f$ is the - * convex hull bounded by the tensor product of nodes \f$x_i\f$, \f$x_{i+1}\f$, - * \f$y_j\f$, and \f$y_{j+1}\f$. - * - * @param[in] nx Number of cells in @c x direction. - * @param[in] ny Number of cells in @c y direction. - * - * @param[in] x Position along @c x axis of each grid line with constant @c x - * coordinate. Array of size nx + 1. - * @param[in] y Position along @c y axis of each grid line with constant @c y - * coordinate. Array of size ny + 1. - * - * @return Fully formed grid structure containing valid geometric primitives. - * Must be destroyed using function destroy_grid(). - */ -struct UnstructuredGrid * -create_grid_tensor2d(int nx, int ny, - const double *x , const double *y ); - - -/** - * Form tensor product (i.e., topologically Cartesian) grid in three space - * dimensions--possibly with a variable top-layer topography. - * - * If @c depthz is @c NULL, then geometric information such as volumes and - * centroids is calculated from analytic expressions. Otherwise, these values - * are computed using function compute_geometry(). - * - * @param[in] nx Number of cells in @c x direction. - * @param[in] ny Number of cells in @c y direction. - * @param[in] nz Number of cells in @c z direction. - * - * @param[in] x Position along @c x axis of each grid line with constant @c x - * coordinate. Array of size nx + 1. - * @param[in] y Position along @c y axis of each grid line with constant @c y - * coordinate. Array of size ny + 1. - * @param[in] z Distance (depth) from top-layer measured along the @c z axis of - * each grid line with constant @c z coordinate. Array of size - * nz + 1. - * - * @param[in] depthz - * Top-layer topography specification. If @c NULL, interpreted as - * horizontal top-layer at z=0. Otherwise, must be - * an array of size (nx + 1) * (ny + 1), ordered - * lexicographically, that defines the depth of each top-layer - * pillar vertex. - * - * @return Fully formed grid structure containing valid geometric primitives. - * Must be destroyed using function destroy_grid(). - */ -struct UnstructuredGrid * -create_grid_tensor3d(int nx , - int ny , - int nz , - const double *x , - const double *y , - const double *z , - const double *depthz); -#ifdef __cplusplus -} -#endif -#endif /* OPM_CART_GRID_H_HEADER */ diff --git a/opm/core/grid/cornerpoint_grid.c b/opm/core/grid/cornerpoint_grid.c deleted file mode 100644 index 8ae8d885..00000000 --- a/opm/core/grid/cornerpoint_grid.c +++ /dev/null @@ -1,235 +0,0 @@ -/*=========================================================================== -// -// File: cgridinterface.c -// -// Author: Jostein R. Natvig -// -//==========================================================================*/ - - -/* - Copyright 2011 SINTEF ICT, Applied Mathematics. -*/ - -#include "config.h" -#include -#include - -#include -#include -#include -#include - - -static int -fill_cell_topology(struct processed_grid *pg, - struct UnstructuredGrid *g ) -{ - int f, c1, c2, tag; - size_t c, nc, nhf; - - nc = g->number_of_cells; - - g->cell_facepos = malloc((nc + 1) * sizeof *g->cell_facepos); - - if (g->cell_facepos != NULL) { - /* Allocate and initialise compressed cell-to-face topology. */ - - for (c = 0; c < nc + 1; c++) { g->cell_facepos[c] = 0; } - - for (f = 0; f < g->number_of_faces; f++) { - c1 = g->face_cells[2*f + 0]; - c2 = g->face_cells[2*f + 1]; - - if (c1 >= 0) { g->cell_facepos[c1 + 1] += 1; } - if (c2 >= 0) { g->cell_facepos[c2 + 1] += 1; } - } - - for (c = 1; c <= nc; c++) { - g->cell_facepos[0] += g->cell_facepos[c]; - g->cell_facepos[c] = g->cell_facepos[0] - g->cell_facepos[c]; - } - - nhf = g->cell_facepos[0]; - g->cell_facepos[0] = 0; - - g->cell_faces = malloc(nhf * sizeof *g->cell_faces ); - g->cell_facetag = malloc(nhf * sizeof *g->cell_facetag); - - if ((g->cell_faces == NULL) || (g->cell_facetag == NULL)) { - free(g->cell_facetag); g->cell_facetag = NULL; - free(g->cell_faces); g->cell_faces = NULL; - free(g->cell_facepos); g->cell_facepos = NULL; - } - } - - if (g->cell_facepos != NULL) { - /* Compute final cell-to-face mapping and half-face tags. - * - * Process relies on preprocess() producing neighbourship - * definitions for which the normals point (logically) in the - * positive I,J,K directions *and* from ->face_cells[2*f+0] to - * ->face_cells[2*f+1] (i.e., from first to second cell of - * interface 'f'--be it internal or outer). - * - * For instance, a "LEFT" connection (pg->face_tag==LEFT==0) - * for which the normal points into the cell (i.e., when - * ->face_cells[2*f+1] >= 0), is a half-face of type 0. - * - * Simlarly, a "TOP" connection (pg->face_tag==TOP==2) for - * which the normal points out of the cell (i.e., when - * ->face_cells[2*f+0] >= 0), is a half-face of type 5. */ - - for (f = 0; f < g->number_of_faces; f++) { - tag = 2 * pg->face_tag[f]; /* [0, 2, 4] */ - c1 = g->face_cells[2*f + 0]; - c2 = g->face_cells[2*f + 1]; - - if (c1 >= 0) { - g->cell_faces [ g->cell_facepos[c1 + 1] ] = f; - g->cell_facetag [ g->cell_facepos[c1 + 1] ] = tag + 1; - - g->cell_facepos[c1 + 1] += 1; - } - - if (c2 >= 0) { - g->cell_faces [ g->cell_facepos[c2 + 1] ] = f; - g->cell_facetag [ g->cell_facepos[c2 + 1] ] = tag + 0; - - g->cell_facepos[c2 + 1] += 1; - } - } - } - - return g->cell_facepos != NULL; -} - -static int -allocate_geometry(struct UnstructuredGrid *g) -{ - int ok; - size_t nc, nf, nd; - - assert (g->dimensions == 3); - - nc = g->number_of_cells; - nf = g->number_of_faces; - nd = 3; - - g->face_areas = malloc(nf * 1 * sizeof *g->face_areas); - g->face_centroids = malloc(nf * nd * sizeof *g->face_centroids); - g->face_normals = malloc(nf * nd * sizeof *g->face_normals); - - g->cell_volumes = malloc(nc * 1 * sizeof *g->cell_volumes); - g->cell_centroids = malloc(nc * nd * sizeof *g->cell_centroids); - - ok = g->face_areas != NULL; - ok += g->face_centroids != NULL; - ok += g->face_normals != NULL; - - ok += g->cell_volumes != NULL; - ok += g->cell_centroids != NULL; - - return ok == 5; -} - - -void compute_geometry(struct UnstructuredGrid *g) -{ - assert (g != NULL); - if (g!=NULL) - { - assert (g->face_centroids != NULL); - assert (g->face_normals != NULL); - assert (g->face_areas != NULL); - assert (g->cell_centroids != NULL); - assert (g->cell_volumes != NULL); - - compute_face_geometry(g->dimensions , g->node_coordinates, - g->number_of_faces, g->face_nodepos, - g->face_nodes, g->face_normals, - g->face_centroids, g->face_areas); - - compute_cell_geometry(g->dimensions, g->node_coordinates, - g->face_nodepos, g->face_nodes, - g->face_cells, g->face_normals, - g->face_centroids, g->number_of_cells, - g->cell_facepos, g->cell_faces, - g->cell_centroids, g->cell_volumes); - } -} - - -struct UnstructuredGrid * -create_grid_cornerpoint(const struct grdecl *in, double tol) -{ - struct UnstructuredGrid *g; - int ok; - struct processed_grid pg; - - g = create_grid_empty(); - if (g == NULL) - { - return NULL; - } - - process_grdecl(in, tol, &pg); - - /* - * Convert "struct processed_grid" to "struct UnstructuredGrid". - * - * In particular, convey resource ownership from 'pg' to 'g'. - * Consequently, memory resources obtained in process_grdecl() - * will be released in destroy_grid(). - */ - g->dimensions = 3; - - g->number_of_nodes = pg.number_of_nodes; - g->number_of_faces = pg.number_of_faces; - g->number_of_cells = pg.number_of_cells; - - g->node_coordinates = pg.node_coordinates; - - g->face_nodes = pg.face_nodes; - g->face_nodepos = pg.face_ptr; - g->face_cells = pg.face_neighbors; - - /* Explicitly relinquish resource references conveyed to 'g'. This - * is needed to avoid creating dangling references in the - * free_processed_grid() call. */ - pg.node_coordinates = NULL; - pg.face_nodes = NULL; - pg.face_ptr = NULL; - pg.face_neighbors = NULL; - - /* allocate and fill g->cell_faces/g->cell_facepos and - * g->cell_facetag as well as the geometry-related fields. */ - ok = fill_cell_topology(&pg, g); - ok = ok && allocate_geometry(g); - - if (!ok) - { - destroy_grid(g); - g = NULL; - } - else - { - - compute_geometry(g); - - g->cartdims[0] = pg.dimensions[0]; - g->cartdims[1] = pg.dimensions[1]; - g->cartdims[2] = pg.dimensions[2]; - - g->global_cell = pg.local_cell_index; - - /* Explicitly relinquish resource references conveyed to 'g'. - * This is needed to avoid creating dangling references in the - * free_processed_grid() call. */ - pg.local_cell_index = NULL; - } - - free_processed_grid(&pg); - - return g; -} diff --git a/opm/core/grid/cornerpoint_grid.h b/opm/core/grid/cornerpoint_grid.h deleted file mode 100644 index 3f371d8f..00000000 --- a/opm/core/grid/cornerpoint_grid.h +++ /dev/null @@ -1,103 +0,0 @@ -/*=========================================================================== -// -// File: preprocess.h -// -// Created: Fri Jun 19 08:43:04 2009 -// -// Author: Jostein R. Natvig -// -// $Date: 2010-08-27 19:12:16 +0200 (Fri, 27 Aug 2010) $ -// -// $Revision: 930 $ -// -//==========================================================================*/ - -/* - Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_CORNERPOINT_GRID_HEADER_INCLUDED -#define OPM_CORNERPOINT_GRID_HEADER_INCLUDED - -/** - * \file - * Routines to form a complete UnstructuredGrid from a corner-point - * specification. - */ - -#include -#include - -#ifdef __cplusplus -extern "C" { -#endif - - /** - * Construct grid representation from corner-point specification of a - * particular geological model. - * - * Pinched cells will be removed irrespective of any explicit "active" map - * in the geological model input specification. Geometric primitives such - * as cell barycenters (i.e., centroids), volumes and interface areas are - * computed internally using function compute_geometry(). The caller does - * not need to compute this information separately. - * - * @param[in] in Corner-point specification. If "actnum" is NULL, then the - * specification is interpreted as if all cells are initially - * active. - * - * @param[in] tol Absolute tolerance of node-coincidence. - * - * @return Fully formed grid data structure that manages the grid defined by - * the input corner-point specification. Must be destroyed using function - * destroy_grid(). - */ - struct UnstructuredGrid * - create_grid_cornerpoint(const struct grdecl *in, double tol); - - - /** - * Compute derived geometric primitives in a grid. - * - * This function computes values for each of the following quantities - * - Quantities pertaining to interfaces (connections, faces) - * -# Barycenters (centroids), g->dimensions scalars per face - * stored sequentially in g->face_centroids. - * -# Areas, one scalar per face stored sequentially in - * g->face_areas. - * -# Normals, g->dimensions scalars per face stored - * sequentially in g->face_normals. The Euclidian norm of - * each normal is equal to the corresponding face's area. - * - * - Quantities pertaining to cells (volumes) - * -# Barycenters (centroids), g->dimensions scalars per cell - * stored sequentially in g->cell_centroids. - * -# Volumes, one scalar per cell stored sequentially in - * g->cell_volumes. - * - * These fields must be allocated prior to calling compute_geometry(). - * - * @param[in,out] g Grid structure. - */ - void compute_geometry(struct UnstructuredGrid *g); - -#ifdef __cplusplus -} -#endif - -#endif /* OPM_CORNERPOINT_GRID_HEADER_INCLUDED */ diff --git a/opm/core/grid/cpgpreprocess/facetopology.c b/opm/core/grid/cpgpreprocess/facetopology.c deleted file mode 100644 index 29561556..00000000 --- a/opm/core/grid/cpgpreprocess/facetopology.c +++ /dev/null @@ -1,369 +0,0 @@ -/*=========================================================================== -// -// File: facetopology.c -// -// Created: Fri Jun 19 08:46:53 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//===========================================================================*/ - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include - -#include "preprocess.h" -#include "facetopology.h" - -/* No checking of input arguments in this code! */ -#define MIN(i,j) ((i)<(j) ? (i) : (j)) -#define MAX(i,j) ((i)>(j) ? (i) : (j)) - -#define DEBUG 1 - -/*------------------------------------------------------*/ -/* */ -/* */ -/* Find connections for each pair of pillars */ -/* */ -/* */ -/* */ -/*------------------------------------------------------*/ - - -/* Determine face topology first, then compute intersection. */ -/* All intersections that occur are present in the final face geometry.*/ -static int * -computeFaceTopology(const int *a1, const int *a2, - const int *b1, const int *b2, - int intersect[4], int *faces) -{ - int mask[8]; - int k; - int *f; - - for (k = 0; k < 8; k++) { mask[k] = -1; } - - /* Which pillar points should we use? */ - if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; } - if (a2[1] > b2[1]){ mask[2] = b2[1]; } else { mask[2] = a2[1]; } - if (a2[0] > b2[0]){ mask[4] = a2[0]; } else { mask[4] = b2[0]; } - if (a1[0] > b1[0]){ mask[6] = a1[0]; } else { mask[6] = b1[0]; } - -#if DEBUG - /* Illegal situations */ - if (mask [0] == mask[2] || - mask [0] == mask[4] || - mask [2] == mask[6] || - mask [4] == mask[6]){ - fprintf(stderr, "Illegal Partial pinch!\n"); - } -#endif - - - /* Partial pinch of face */ - if (mask[0] == mask[6]){ - mask[6] = -1; - } - - - if (mask[2] == mask[4]){ - mask[4] = -1; - } - - - - /* Get shape of face: */ - /* each new intersection will be part of the new face, */ - /* but not all pillar points. This is encoded in mask. */ - - - mask[1] = intersect[3]; /* top-top */ - mask[3] = -1; - mask[5] = intersect[0]; /* bottom-bottom*/ - mask[7] = -1; - - /* bottom-top */ - if (intersect[1] != -1){ - if(a1[0] > b1[1]){ /* intersection[1] left of (any) intersection[0] */ - mask[0] = -1; - mask[6] = -1; - mask[7] = intersect[1]; - } - else{ - mask[2] = -1; - mask[4] = -1; - mask[3] = intersect[1]; - } - } - - /* top-bottom */ - if (intersect[2] != -1){ - if(a1[1] < b1[0]){ /* intersection[2] left of (any) intersection[3] */ - mask[0] = -1; - mask[6] = -1; - mask[7] = intersect[2]; - } - else{ - mask[2] = -1; - mask[4] = -1; - mask[3] = intersect[2]; - } - } - f = faces; - for (k=7; k>=0; --k){ - if(mask[k] != -1){ - *f++ = mask[k]; - } - } - -#if DEBUG>1 - /* Check for repeated nodes:*/ - int i; - fprintf(stderr, "face: "); - for (i=0; i<8; ++i){ - fprintf(stderr, "%d ", mask[i]); - for (k=0; k<8; ++k){ - if (i!=k && mask[i] != -1 && mask[i] == mask[k]){ - fprintf(stderr, "Repeated node in faulted face\n"); - } - } - } - fprintf(stderr, "\n"); -#endif - - - return f; - - - -} - - - - -/* a) If we assume that the index increase when z increase for each - pillar (but only separately), we can use only the point indices, - since we only need to compare z-values on one pillar at a time. - - b) We assume input is preprocessed such that no intersections occur - on the first and last lines, for instance by padding the grid with - extra cells. This is convenient in the identification of (unique) - intersections. - -*/ - -#define LINE_INTERSECTION(a1, a2, b1, b2) \ - (((a1 > b1) && (a2 < b2)) || \ - ((a1 < b1) && (a2 > b2))) - -static int -faceintersection(const int *a1, const int *a2, - const int *b1, const int *b2) -{ - return - MAX(a1[0], b1[0]) < MIN(a1[1], b1[1]) || - MAX(a2[0], b2[0]) < MIN(a2[1], b2[1]) || - LINE_INTERSECTION(a1[0], a2[0], b1[0], b2[0]) || - LINE_INTERSECTION(a1[1], a2[1], b1[1], b2[1]); -} - - -#define MEANINGFUL_FACE(i, j) \ - (! ((a1[i] == INT_MIN) && (b1[j] == INT_MIN)) && \ - ! ((a1[i+1] == INT_MAX) && (b1[j+1] == INT_MAX))) - - -/* work should be pointer to 2n ints initialised to zero . */ -void findconnections(int n, int *pts[4], - int *intersectionlist, - int *work, - struct processed_grid *out) -{ - /* vectors of point numbers for faces a(b) on pillar 1(2) */ - int *a1 = pts[0]; - int *a2 = pts[1]; - int *b1 = pts[2]; - int *b2 = pts[3]; - - /* Intersection record for top line and bottomline of a */ - int *itop = work; - int *ibottom = work + n; - int *f = out->face_nodes + out->face_ptr[out->number_of_faces]; - int *c = out->face_neighbors + 2*out->number_of_faces; - - int k1 = 0; - int k2 = 0; - - int i,j=0; - int intersect[4]; - int *tmp; - /* for (i=0; i<2*n; work[i++]=-1); */ - - for (i = 0; i < 4; i++) { intersect[i] = -1; } - - for (i = 0; i < n - 1; ++i) { - - /* pinched a-cell */ - if ((a1[i] == a1[i + 1]) && - (a2[i] == a2[i + 1])) { - continue; - } - - - while ((j < n-1) && - ((b1[j] < a1[i + 1]) || - (b2[j] < a2[i + 1]))) - { - /* pinched b-cell */ - if ((b1[j] == b1[j + 1]) && - (b2[j] == b2[j + 1])) { - - itop[j+1] = itop[j]; - ++j; - continue; - } - - - /* --------------------------------------------------------- */ - /* face a(i,i+1) and face b(j,j+1) have nonzero intersection */ - /* --------------------------------------------------------- */ - if (faceintersection(a1+i, a2+i, b1+j, b2+j)){ - - - /* Completely matching faces */ - if (((a1[i] == b1[j]) && (a1[i+1] == b1[j+1])) && - ((a2[i] == b2[j]) && (a2[i+1] == b2[j+1]))) { - - if (MEANINGFUL_FACE(i, j)) { - - int cell_a = i%2 != 0 ? (i-1)/2 : -1; - int cell_b = j%2 != 0 ? (j-1)/2 : -1; - - if (cell_a != -1 || cell_b != -1){ - *c++ = cell_a; - *c++ = cell_b; - - /* face */ - *f++ = a1[i]; - *f++ = a2[i]; - - /* avoid duplicating nodes in pinched faces */ - if (a2[i+1] != a2[i]) { *f++ = a2[i+1]; } - if (a1[i+1] != a1[i]) { *f++ = a1[i+1]; } - - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - - } - else{ - ; - } - } - } - - /* Non-matching faces */ - else{ - - /* Find new intersection */ - if (LINE_INTERSECTION(a1[i+1], a2[i+1], - b1[j+1], b2[j+1])) { - itop[j+1] = out->number_of_nodes++; - - /* store point numbers of intersecting lines */ - *intersectionlist++ = a1[i+1]; - *intersectionlist++ = a2[i+1]; - *intersectionlist++ = b1[j+1]; - *intersectionlist++ = b2[j+1]; - - - }else{ - itop[j+1] = -1; - } - - /* Update intersection record */ - intersect[0] = ibottom[j ]; /* i x j */ - intersect[1] = ibottom[j+1]; /* i x j+1 */ - intersect[2] = itop[j ]; /* i+1 x j */ - intersect[3] = itop[j+1]; /* i+1 x j+1 */ - - - /* Add face to list of faces if no INT_MIN or - * INT_MAX appear in a or b. */ - if (MEANINGFUL_FACE(i,j)) { - - /* - Even indices refer to space between cells, - odd indices refer to cells - */ - int cell_a = i%2 != 0 ? (i-1)/2 : -1; - int cell_b = j%2 != 0 ? (j-1)/2 : -1; - - - - if (cell_a != -1 || cell_b != -1){ - *c++ = cell_a; - *c++ = cell_b; - - f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f); - - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - } - else{ - ; - } - } - } - } - - /* Update candidates for restart of j for in next i-iteration */ - if (b1[j] < a1[i+1]) { k1 = j; } - if (b2[j] < a2[i+1]) { k2 = j; } - - j = j+1; - } - - - - /* Swap intersection records: top line of a[i,i+1] is bottom - * line of a[i+1,i+2] */ - tmp = itop; itop = ibottom; ibottom = tmp; - - /* Zero out the "new" itop */ - for(j=0;j -// -// $Date$ -// -// $Revision$ -// -//===========================================================================*/ - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_FACETOPOLOGY_HEADER -#define OPM_FACETOPOLOGY_HEADER - - -void findconnections(int n, int *pts[4], - int *intersectionlist, - int *work, - struct processed_grid *out); - -#endif /* OPM_FACETOPOLOGY_HEADER */ - -/* Local Variables: */ -/* c-basic-offset:4 */ -/* End: */ diff --git a/opm/core/grid/cpgpreprocess/geometry.c b/opm/core/grid/cpgpreprocess/geometry.c deleted file mode 100644 index e010935c..00000000 --- a/opm/core/grid/cpgpreprocess/geometry.c +++ /dev/null @@ -1,460 +0,0 @@ -/* - * Copyright 2010 (c) SINTEF ICT, Applied Mathematics. - * Jostein R. Natvig - */ -#include "config.h" -#include -#include -#include "geometry.h" -#include - -/* ------------------------------------------------------------------ */ -static void -cross(const double u[3], const double v[3], double w[3]) -/* ------------------------------------------------------------------ */ -{ - w[0] = u[1]*v[2]-u[2]*v[1]; - w[1] = u[2]*v[0]-u[0]*v[2]; - w[2] = u[0]*v[1]-u[1]*v[0]; -} - -/* ------------------------------------------------------------------ */ -static double -norm(const double w[3]) -/* ------------------------------------------------------------------ */ -{ - return sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]); -} - - -/* ------------------------------------------------------------------ */ -static void -compute_face_geometry_3d(double *coords, int nfaces, - int *nodepos, int *facenodes, double *fnormals, - double *fcentroids, double *fareas) -/* ------------------------------------------------------------------ */ -{ - - /* Assume 3D for now */ - const int ndims = 3; - int f; - double x[3]; - double u[3]; - double v[3]; - double w[3]; - - int i,k; - int node; - - double cface[3] = {0}; - double n[3] = {0}; - double twothirds = 0.666666666666666666666666666667; - double a; - int num_face_nodes; - double area; - /*#pragma omp parallel for */ - - /*#pragma omp parallel for shared(fnormals,fcentroids,fareas)*/ -#pragma omp parallel for default(none) \ - private(f,x,u,v,w,i,k,node,cface,n,a,num_face_nodes,area) \ - shared(fnormals,fcentroids,fareas \ - ,coords, nfaces, nodepos, facenodes, twothirds) - for (f=0; f0)) - { - fprintf(stderr, "Internal error in compute_face_geometry."); - } - */ - /* face normal */ - for (i=0; i - */ -#ifndef MRST_GEOMETRY_H_INCLUDED -#define MRST_GEOMETRY_H_INCLUDED - -void compute_face_geometry(int ndims, double *coords, int nfaces, - int *nodepos, int *facenodes, - double *fnormals, double *fcentroids, - double *fareas); -void compute_cell_geometry(int ndims, double *coords, - int *nodepos, int *facenodes, int *neighbours, - double *fnormals, - double *fcentroids, int ncells, - int *facepos, int *cellfaces, - double *ccentroids, double *cvolumes); - -#endif /* MRST_GEOMETRY_H_INCLUDED */ diff --git a/opm/core/grid/cpgpreprocess/preprocess.c b/opm/core/grid/cpgpreprocess/preprocess.c deleted file mode 100644 index 3b465c2a..00000000 --- a/opm/core/grid/cpgpreprocess/preprocess.c +++ /dev/null @@ -1,956 +0,0 @@ -/*=========================================================================== -// -// File: preprocess.c -// -// Created: Fri Jun 19 08:42:39 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//==========================================================================*/ - -/* - Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010, 2011, 2012 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include - -#include "preprocess.h" -#include "uniquepoints.h" -#include "facetopology.h" - - -#define MIN(i,j) ((i)<(j) ? (i) : (j)) -#define MAX(i,j) ((i)>(j) ? (i) : (j)) - -static void -compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len); - -static int -checkmemory(int nz, struct processed_grid *out, int **intersections); - -static void -process_vertical_faces(int direction, - int **intersections, - int *plist, int *work, - struct processed_grid *out); - -static void -process_horizontal_faces(int **intersections, - int *plist, - struct processed_grid *out); - -static int -linearindex(const int dims[3], int i, int j, int k) -{ - assert (0 <= i); - assert (0 <= j); - assert (0 <= k); - - assert (i < dims[0]); - assert (j < dims[1]); - assert (k < dims[2]); - - return i + dims[0]*(j + dims[1]*k); -} - - -/*----------------------------------------------------------------- - Given a vector with k index running faster than i running - faster than j, and Cartesian dimensions , find pointers to the - (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. */ -static void -igetvectors(int dims[3], int i, int j, int *field, int *v[]) -{ - int im = MAX(1, i ) - 1; - int ip = MIN(dims[0], i+1) - 1; - int jm = MAX(1, j ) - 1; - int jp = MIN(dims[1], j+1) - 1; - - v[0] = field + dims[2]*(im + dims[0]* jm); - v[1] = field + dims[2]*(im + dims[0]* jp); - v[2] = field + dims[2]*(ip + dims[0]* jm); - v[3] = field + dims[2]*(ip + dims[0]* jp); -} - - - - - - -/*----------------------------------------------------------------- - Special purpose - - Convert from k-index to Cartesian index i+nx*(j+ny*k) for every - other element in neighbors. - -*/ -static void -compute_cell_index(const int dims[3], int i, int j, - int *neighbors, int len) -{ - int k; - - if (((i < 0) || (i >= dims[0])) || /* 'i' outside [0, dims[0]) */ - ((j < 0) || (j >= dims[1]))) { /* 'j' outside [0, dims[1]) */ - - for (k = 0; k < len; k += 2) { - neighbors[k] = -1; /* Neighbour is outside domain */ - } - } - else { - for (k = 0; k < len; k += 2) { - if (neighbors[k] != -1) { - neighbors[k] = linearindex(dims, i, j, neighbors[k]); - } - } - } -} - - -/*----------------------------------------------------------------- - Ensure there's sufficient memory */ -static int -checkmemory(int nz, struct processed_grid *out, int **intersections) -{ - int r, m, n, ok; - - /* Ensure there is enough space to manage the (pathological) case - * of every single cell on one side of a fault connecting to all - * cells on the other side of the fault (i.e., an all-to-all cell - * connectivity pairing). */ - r = (2*nz + 2) * (2*nz + 2); - m = out->m; - n = out->n; - - if (out->number_of_faces + r > m) { - m += MAX(m / 2, 2 * r); - } - if (out->face_ptr[out->number_of_faces] + 6*r > n) { - n += MAX(n / 2, 12 * r); - } - - ok = m == out->m; - if (! ok) { - void *p1, *p2, *p3, *p4; - - p1 = realloc(*intersections , 4*m * sizeof **intersections); - p2 = realloc(out->face_neighbors, 2*m * sizeof *out->face_neighbors); - p3 = realloc(out->face_ptr , (m+1) * sizeof *out->face_ptr); - p4 = realloc(out->face_tag , 1*m * sizeof *out->face_tag); - - if (p1 != NULL) { *intersections = p1; } - if (p2 != NULL) { out->face_neighbors = p2; } - if (p3 != NULL) { out->face_ptr = p3; } - if (p4 != NULL) { out->face_tag = p4; } - - ok = (p1 != NULL) && (p2 != NULL) && (p3 != NULL) && (p4 != NULL); - - if (ok) { out->m = m; } - } - - if (ok && (n != out->n)) { - void *p1; - - p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes); - - ok = p1 != NULL; - - if (ok) { - out->face_nodes = p1; - out->n = n; - } - } - - return ok; -} - -/*----------------------------------------------------------------- - For each vertical face (i.e. i or j constant), - -find point numbers for the corners and - -cell neighbors. - -new points on faults defined by two intgersecting lines. - - direction == 0 : constant-i faces. - direction == 1 : constant-j faces. -*/ -static void -process_vertical_faces(int direction, - int **intersections, - int *plist, int *work, - struct processed_grid *out) -{ - int i,j; - int *cornerpts[4]; - int d[3]; - int f; - enum face_tag tag[] = { LEFT, BACK }; - int *tmp; - int nx = out->dimensions[0]; - int ny = out->dimensions[1]; - int nz = out->dimensions[2]; - int startface; - int num_intersections; - int *ptr; - int len; - - assert ((direction == 0) || (direction == 1)); - - d[0] = 2 * (nx + 0); - d[1] = 2 * (ny + 0); - d[2] = 2 * (nz + 1); - - for (j = 0; j < ny + direction; ++j) { - for (i = 0; i < nx + (1 - direction); ++i) { - - if (! checkmemory(nz, out, intersections)) { - fprintf(stderr, - "Could not allocate enough space in " - "process_vertical_faces()\n"); - exit(1); - } - - /* Vectors of point numbers */ - igetvectors(d, 2*i + direction, 2*j + (1 - direction), - plist, cornerpts); - - if (direction == 1) { - /* 1 3 0 1 */ - /* ---> */ - /* 0 2 2 3 */ - /* rotate clockwise */ - tmp = cornerpts[1]; - cornerpts[1] = cornerpts[0]; - cornerpts[0] = cornerpts[2]; - cornerpts[2] = cornerpts[3]; - cornerpts[3] = tmp; - } - - /* int startface = ftab->position; */ - startface = out->number_of_faces; - /* int num_intersections = *npoints - npillarpoints; */ - num_intersections = out->number_of_nodes - - out->number_of_nodes_on_pillars; - - /* Establish new connections (faces) along pillar pair. */ - findconnections(2*nz + 2, cornerpts, - *intersections + 4*num_intersections, - work, out); - - /* Start of ->face_neighbors[] for this set of connections. */ - ptr = out->face_neighbors + 2*startface; - - /* Total number of cells (both sides) connected by this - * set of connections (faces). */ - len = 2*out->number_of_faces - 2*startface; - - /* Derive inter-cell connectivity (i.e. ->face_neighbors) - * of global (uncompressed) cells for this set of - * connections (faces). */ - compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr , len); - compute_cell_index(out->dimensions, i , j , ptr + 1, len); - - /* Tag the new faces */ - f = startface; - for (; f < out->number_of_faces; ++f) { - out->face_tag[f] = tag[direction]; - } - } - } -} - - -/*----------------------------------------------------------------- - For each horizontal face (i.e. k constant), - -find point numbers for the corners and - -cell neighbors. - - Also define map from logically Cartesian - cell index to local cell index 0, ..., #. Exclude - cells that are have collapsed coordinates. (This includes cells with - ACTNUM==0) - -*/ -static void -process_horizontal_faces(int **intersections, - int *plist, - struct processed_grid *out) -{ - int i,j,k; - - int nx = out->dimensions[0]; - int ny = out->dimensions[1]; - int nz = out->dimensions[2]; - - int *cell = out->local_cell_index; - int cellno = 0; - int *f, *n, *c[4]; - int prevcell, thiscell; - int idx; - - /* dimensions of plist */ - int d[3]; - d[0] = 2*nx; - d[1] = 2*ny; - d[2] = 2+2*nz; - - - for(j=0; jface_nodes + out->face_ptr[out->number_of_faces]; - n = out->face_neighbors + 2*out->number_of_faces; - - - /* Vectors of point numbers */ - igetvectors(d, 2*i+1, 2*j+1, plist, c); - - prevcell = -1; - - - for (k = 1; kdimensions, i,j,(k-1)/2); - cell[idx] = -1; - } - } - else{ - - if (k%2){ - /* Add face */ - *f++ = c[0][k]; - *f++ = c[2][k]; - *f++ = c[3][k]; - *f++ = c[1][k]; - - out->face_tag[ out->number_of_faces] = TOP; - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - - thiscell = linearindex(out->dimensions, i,j,(k-1)/2); - *n++ = prevcell; - *n++ = prevcell = thiscell; - - cell[thiscell] = cellno++; - - } - else{ - if (prevcell != -1){ - /* Add face */ - *f++ = c[0][k]; - *f++ = c[2][k]; - *f++ = c[3][k]; - *f++ = c[1][k]; - - out->face_tag[ out->number_of_faces] = TOP; - out->face_ptr[++out->number_of_faces] = f - out->face_nodes; - - *n++ = prevcell; - *n++ = prevcell = -1; - } - } - } - } - } - } - out->number_of_cells = cellno; -} - - -/*----------------------------------------------------------------- - On input, - L points to 4 ints that indirectly refers to points in c. - c points to array of coordinates [x0,y0,z0,x1,y1,z1,...,xn,yn,zn]. - pt points to array of 3 doubles. - - On output, - pt holds coordinates to intersection between lines given by point - numbers L[0]-L[1] and L[2]-L[3]. -*/ -static void approximate_intersection_pt(int *L, double *c, double *pt) -{ - double a; - double z0, z1, z2, z3; - double b1, b2; - double x1, y1; - double x2, y2; - double z; - - /* no intersection on pillars expected here! */ - assert (L[0] != L[2]); - assert (L[1] != L[3]); - - z0 = c[3*L[0] + 2]; - z1 = c[3*L[1] + 2]; - z2 = c[3*L[2] + 2]; - z3 = c[3*L[3] + 2]; - - /* find parameter a where lines L0L1 and L2L3 have same - * z-coordinate */ - if (fabs((z1 - z0) - (z3 - z2)) > 0.0) { - - a = (z2 - z0) / ((z1 - z0) - (z3 - z2)); - - } else { - - a = 0; - - } - - /* the corresponding z-coordinate is */ - z = z0*(1.0 - a) + z1*a; - - - /* find point (x1, y1, z) on pillar 1 */ - b1 = (z2 - z) / (z2 - z0); - b2 = (z - z0) / (z2 - z0); - x1 = c[3*L[0] + 0]*b1 + c[3*L[2] + 0]*b2; - y1 = c[3*L[0] + 1]*b1 + c[3*L[2] + 1]*b2; - - /* find point (x2, y2, z) on pillar 2 */ - b1 = (z - z3) / (z1 - z3); - b2 = (z1 - z) / (z1 - z3); - x2 = c[3*L[1] + 0]*b1 + c[3*L[3] + 0]*b2; - y2 = c[3*L[1] + 1]*b1 + c[3*L[3] + 1]*b2; - - /* horizontal lines are by definition ON the bilinear surface - spanned by L0, L1, L2 and L3. find point (x, y, z) on - horizontal line between point (x1, y1, z) and (x2, y2, z).*/ - pt[0] = x1*(1.0 - a) + x2*a; - pt[1] = y1*(1.0 - a) + y2*a; - pt[2] = z; -} - -/*----------------------------------------------------------------- - Compute x,y and z coordinates for points on each pillar. Then, - append x,y and z coordinates for extra points on faults. */ -static void -compute_intersection_coordinates(int *intersections, - struct processed_grid *out) -{ - int n = out->number_of_nodes; - int np = out->number_of_nodes_on_pillars; - int k; - double *pt; - int *itsct = intersections; - /* Make sure the space allocated for nodes match the number of - * node. */ - void *p = realloc (out->node_coordinates, 3*n*sizeof(double)); - if (p) { - out->node_coordinates = p; - } - else { - fprintf(stderr, "Could not allocate extra space for intersections\n"); - } - - - /* Append intersections */ - pt = out->node_coordinates + 3*np; - - for (k=np; knode_coordinates, pt); - pt += 3; - itsct += 4; - - } -} - - -/* ------------------------------------------------------------------ */ -static int* -copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out) -/* ------------------------------------------------------------------ */ -{ - int i,j,k; - int *ptr = out; - - /* Permute actnum such that values of each vertical stack of cells - * are adjacent in memory, i.e., - * - * out = [in(0,0,:), in(1,0,:),..., in(nx-1, ny-1,:)] - * - * in MATLAB pseudo-code. - */ - if (in != NULL) { - for (j = 0; j < ny; ++j) { - for (i = 0; i < nx; ++i) { - for (k = 0; k < nz; ++k) { - *ptr++ = in[i + nx*(j + ny*k)]; - } - } - } - } - else { - /* No explicit ACTNUM. Assume all cells active. */ - for (i = 0; i < nx * ny * nz; i++) { - out[ i ] = 1; - } - } - - return out; -} - -/* ------------------------------------------------------------------ */ -static double* -copy_and_permute_zcorn(int nx, int ny, int nz, const double *in, - double sign, double *out) -/* ------------------------------------------------------------------ */ -{ - int i,j,k; - double *ptr = out; - /* Permute zcorn such that values of each vertical stack of cells - * are adjacent in memory, i.e., - - out = [in(0,0,:), in(1,0,:),..., in(2*nx-1, 2*ny-1,:)] - - in Matlab pseudo-code. - */ - for (j=0; j<2*ny; ++j){ - for (i=0; i<2*nx; ++i){ - for (k=0; k<2*nz; ++k){ - *ptr++ = sign * in[i+2*nx*(j+2*ny*k)]; - } - } - } - return out; -} - -/* ------------------------------------------------------------------ */ -static int -get_zcorn_sign(int nx, int ny, int nz, const int *actnum, - const double *zcorn, int *error) -/* ------------------------------------------------------------------ */ -{ - /* Ensure that zcorn (i.e., depth) is strictly nondecreasing in - the k-direction. This is required by the processign algorithm. - - 1) if z(i,j,k) <= z(i,j,k+1) for all (i,j,k), return 1.0 - - 2) if -z(i,j,k) <=-z(i,j,k+1) for all (i,j,k), return -1.0 - - 3) if (1) and (2) fails, return -1.0, and set *error = 1. - - */ - int sign; - int i, j, k; - int c1, c2; - double z1, z2; - - for (sign = 1; sign>-2; sign = sign - 2) - { - *error = 0; - - for (j=0; j<2*ny; ++j){ - for (i=0; i<2*nx; ++i){ - for (k=0; k<2*nz-1; ++k){ - z1 = sign*zcorn[i+2*nx*(j+2*ny*(k))]; - z2 = sign*zcorn[i+2*nx*(j+2*ny*(k+1))]; - - c1 = i/2 + nx*(j/2 + ny*(k/2)); - c2 = i/2 + nx*(j/2 + ny*((k+1)/2)); - - assert (c1 < (nx * ny * nz)); - assert (c2 < (nx * ny * nz)); - - if (((actnum == NULL) || - (actnum[c1] && actnum[c2])) - && (z2 < z1)) { - - fprintf(stderr, "\nZCORN should be strictly " - "nondecreasing along pillars!\n"); - *error = 1; - goto end; - } - } - } - } - - end: - if (!*error){ - break; - } - } - - if (*error){ - fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n" - "Grid definition may be broken\n"); - } - - return sign; -} - - -static void -ind2sub(const size_t nx, - const size_t ny, - const size_t nz, - size_t c , - size_t *i, size_t *j, size_t *k) -{ - assert (c < (nx * ny * nz)); - -#if defined(NDEBUG) - (void) nz; -#endif - - *i = c % nx; c /= nx; - *j = c % ny; - *k = c / ny; -} - - -/* ---------------------------------------------------------------------- */ -static double -vert_size(const struct grdecl *in, - const size_t c , - const size_t off[8]) -/* ---------------------------------------------------------------------- */ -{ - size_t i, j, k, nx, ny, start; - double dz; - const double *zcorn; - - nx = in->dims[ 0 ]; - ny = in->dims[ 1 ]; - - ind2sub(nx, ny, in->dims[ 2 ], c, &i, &j, &k); - - zcorn = in->zcorn; - start = (2 * i) + (2 * nx)*((2 * j) + (2 * ny)*(2 * k)); - - for (k = 0, dz = 0.0; (! (fabs(dz) > 0)) && (k < 4); k++) { - dz = zcorn[start + off[k + 4]] - zcorn[start + off[k]]; - } - - return dz; -} - - -/* ---------------------------------------------------------------------- */ -static int -is_lefthanded(const struct grdecl *in) -/* ---------------------------------------------------------------------- */ -{ - int active, searching; - size_t nx, ny, nz, c; - size_t origin, imax, jmax; - size_t off[8]; - double dx[2], dy[2], dz, triple; - const double *pt_coord; - - nx = in->dims[0]; - ny = in->dims[1]; - nz = in->dims[2]; - - off[0] = 0; - off[1] = off[0] + 1; - off[2] = off[0] + (2 * nx); - off[3] = off[2] + 1; - off[4] = off[0] + ((2 * nx) * (2 * ny)); - off[5] = off[4] + 1; - off[6] = off[4] + (2 * nx); - off[7] = off[6] + 1; - - pt_coord = in->coord; - - origin = 0; - imax = (nx + 0) * 1 * (2 * 3); - jmax = (nx + 1) * (ny + 0) * (2 * 3); - - dx[0] = pt_coord[imax + 0] - pt_coord[origin + 0]; - dy[0] = pt_coord[imax + 1] - pt_coord[origin + 1]; - - dx[1] = pt_coord[jmax + 0] - pt_coord[origin + 0]; - dy[1] = pt_coord[jmax + 1] - pt_coord[origin + 1]; - - c = 0; dz = 0.0; - do { - active = (in->actnum == NULL) || (in->actnum[c] != 0); - - if (active) { - dz = vert_size(in, c, off); - } - - searching = ! (active && (fabs(dz) > 0.0)); - - c += 1; - } while (searching && (c < (nx * ny * nz))); - - assert (! searching); /* active && (fabs(dz) > 0) */ - - /* Compute vector triple product to distinguish left-handed (<0) - * from right-handed (>0) coordinate systems. */ - triple = dz * (dx[0]*dy[1] - dx[1]*dy[0]); - - assert (fabs(triple) > 0.0); - - return triple < 0.0; -} - - -/* ---------------------------------------------------------------------- */ -static void -reverse_face_nodes(struct processed_grid *out) -/* ---------------------------------------------------------------------- */ -{ - int f, t, *i, *j; - - for (f = 0; f < out->number_of_faces; f++) { - i = out->face_nodes + (out->face_ptr[f + 0] + 0); - j = out->face_nodes + (out->face_ptr[f + 1] - 1); - - assert (i <= j); - - while (i < j) { - t = *i; - *i = *j; - *j = t; - - i += 1; - j -= 1; - } - } -} - - -/*----------------------------------------------------------------- - Public interface -*/ -void process_grdecl(const struct grdecl *in, - double tolerance, - struct processed_grid *out) -{ - struct grdecl g; - - size_t i; - int sign, error, left_handed; - int cellnum; - - int *actnum, *iptr; - int *global_cell_index; - - double *zcorn; - - const size_t BIGNUM = 64; - const int nx = in->dims[0]; - const int ny = in->dims[1]; - const int nz = in->dims[2]; - const size_t nc = ((size_t) nx) * ((size_t) ny) * ((size_t) nz); - - /* internal work arrays */ - int *work; - int *plist; - int *intersections; - - - - - /* -----------------------------------------------------------------*/ - /* Initialize output structure: - 1) allocate space for grid topology (which may need to be - increased) - 2) set Cartesian imensions - */ - out->m = (int) (BIGNUM / 3); - out->n = (int) BIGNUM; - - out->face_neighbors = malloc( BIGNUM * sizeof *out->face_neighbors); - out->face_nodes = malloc( out->n * sizeof *out->face_nodes); - out->face_ptr = malloc((out->m + 1) * sizeof *out->face_ptr); - out->face_tag = malloc( out->m * sizeof *out->face_tag); - out->face_ptr[0] = 0; - - out->dimensions[0] = in->dims[0]; - out->dimensions[1] = in->dims[1]; - out->dimensions[2] = in->dims[2]; - out->number_of_faces = 0; - out->number_of_nodes = 0; - out->number_of_cells = 0; - - out->node_coordinates = NULL; - out->local_cell_index = malloc(nc * sizeof *out->local_cell_index); - - - - /* Do actual work here:*/ - - /* -----------------------------------------------------------------*/ - /* For each pillar, compare zcorn values for adjacent cells to - * find the unique node z-coordinates specified by the input. - * While here, enumerate unique points and assign point numbers - * (in plist) for each cornerpoint cell. In other words, plist has - * 8 node numbers for each cornerpoint cell.*/ - - /* initialize grdecl structure "g" that will be processd by - * "finduniquepoints" */ - g.dims[0] = in->dims[0]; - g.dims[1] = in->dims[1]; - g.dims[2] = in->dims[2]; - - actnum = malloc (nc * sizeof *actnum); - g.actnum = copy_and_permute_actnum(nx, ny, nz, in->actnum, actnum); - - zcorn = malloc (nc * 8 * sizeof *zcorn); - sign = get_zcorn_sign(nx, ny, nz, in->actnum, in->zcorn, &error); - g.zcorn = copy_and_permute_zcorn(nx, ny, nz, in->zcorn, sign, zcorn); - - g.coord = in->coord; - - - /* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) - * padding */ - plist = malloc(8 * (nc + ((size_t)nx)*((size_t)ny)) * sizeof *plist); - - finduniquepoints(&g, plist, tolerance, out); - - free (zcorn); - free (actnum); - - /* Determine if coordinate system is left handed or not. */ - left_handed = is_lefthanded(in); - if (left_handed) { - /* Reflect Y coordinates about XZ plane to create right-handed - * coordinate system whilst processing intersections. */ - for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) { - out->node_coordinates[i] = -out->node_coordinates[i]; - } - } - - /* -----------------------------------------------------------------*/ - /* Find face topology and face-to-cell connections */ - - /* internal */ - work = malloc(2 * ((size_t) (2*nz + 2)) * sizeof *work); - for(i = 0; i < ((size_t)4) * (nz + 1); ++i) { work[i] = -1; } - - /* internal array to store intersections */ - intersections = malloc(BIGNUM* sizeof(*intersections)); - - - - process_vertical_faces (0, &intersections, plist, work, out); - process_vertical_faces (1, &intersections, plist, work, out); - process_horizontal_faces ( &intersections, plist, out); - - free (plist); - free (work); - - /* -----------------------------------------------------------------*/ - /* (re)allocate space for and compute coordinates of nodes that - * arise from intersecting cells (faults) */ - compute_intersection_coordinates(intersections, out); - - free (intersections); - - /* -----------------------------------------------------------------*/ - /* Enumerate compressed cells: - -make array [0...#cells-1] of global cell numbers - -make [0...nx*ny*nz-1] array of local cell numbers, - lexicographically ordered, used to remap out->face_neighbors - */ - global_cell_index = malloc(nc * sizeof *global_cell_index); - cellnum = 0; - for (i = 0; i < nc; ++i) { - if (out->local_cell_index[i] != -1) { - global_cell_index[cellnum] = (int) i; - out->local_cell_index[i] = cellnum; - cellnum++; - } - } - - /* Remap out->face_neighbors */ - iptr = out->face_neighbors; - for (i = 0; i < ((size_t) 2) * out->number_of_faces; ++i, ++iptr) { - if (*iptr != -1){ - *iptr = out->local_cell_index[*iptr]; - } - } - - - free(out->local_cell_index); - out->local_cell_index = global_cell_index; - - /* Reflect Y coordinate back to original position if left-handed - * coordinate system was detected and handled earlier. */ - if (left_handed) { - for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) { - out->node_coordinates[i] = -out->node_coordinates[i]; - } - } - - /* if sign==-1 in ZCORN preprocessing, the sign of the - * z-coordinate need to change before we finish */ - if (sign == -1) - { - for (i = 2; i < ((size_t) 3) * out->number_of_nodes; i += 3) - out->node_coordinates[i] *= sign; - } - - /* If an odd number of coordinate reflections were applied, the - * processing routines--especially facetopology()--will produce - * node orderings that lead to normals pointing from 2 to 1. - * Reverse nodes to reestablish expected normal direction (and - * positive cell volumes). */ - if (left_handed ^ (sign == -1)) { - reverse_face_nodes(out); - } -} - -/*-------------------------------------------------------*/ -void free_processed_grid(struct processed_grid *g) -{ - if( g ){ - free ( g->face_nodes ); - free ( g->face_ptr ); - free ( g->face_tag ); - free ( g->face_neighbors ); - free ( g->node_coordinates ); - free ( g->local_cell_index ); - } -} - -/* Local Variables: */ -/* c-basic-offset:4 */ -/* End: */ diff --git a/opm/core/grid/cpgpreprocess/preprocess.h b/opm/core/grid/cpgpreprocess/preprocess.h deleted file mode 100644 index 499cf35e..00000000 --- a/opm/core/grid/cpgpreprocess/preprocess.h +++ /dev/null @@ -1,154 +0,0 @@ -/*=========================================================================== -// -// File: preprocess.h -// -// Created: Fri Jun 19 08:43:04 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//==========================================================================*/ - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_PREPROCESS_HEADER -#define OPM_PREPROCESS_HEADER - -/** - * \file - * Low-level corner-point processing routines and supporting data structures. - * - * User code should typically employ higher-level routines such as - * create_grid_cornerpoint() in order to construct fully formed UnstructuredGrid - * data structures from a corner-point specification. Incidentally, the routines - * provided by this module are used to implement function - * create_grid_cornerpoint(). - */ - -#ifdef __cplusplus -extern "C" { -#endif - - /** - * Raw corner-point specification of a particular geological model. - */ - struct grdecl { - int dims[3]; /**< Cartesian box dimensions. */ - const double *coord; /**< Pillar end-points. */ - const double *zcorn; /**< Corner-point depths. */ - const int *actnum; /**< Explicit "active" map. May be NULL.*/ - const double *mapaxes; /**< 6 Element rotation vector - can be NULL. */ - }; - - /** - * Connection taxonomy. - */ - enum face_tag { - LEFT, /**< Connection topologically parallel to J-K plane. */ - BACK, /**< Connection topologically parallel to I-K plane. */ - TOP /**< Connection topologically parallel to I-J plane. */ - }; - - - /** - * Result structure representing minimal derived topology and geometry of - * a geological model in corner-point format. - */ - struct processed_grid { - int m; /**< Upper bound on "number_of_faces". For internal use in - function process_grid()'s memory management. */ - int n; /**< Upper bound on "number_of_nodes". For internal use in - function process_grid()'s memory management. */ - - int dimensions[3]; /**< Cartesian box dimensions. */ - - int number_of_faces; /**< Total number of unique grid faces - (i.e., connections). */ - int *face_nodes; /**< Node (vertex) numbers of each face, - stored sequentially. */ - int *face_ptr; /**< Start position for each face's - `face_nodes'. */ - int *face_neighbors; /**< Global cell numbers. Two elements per - face, stored sequentially. */ - enum face_tag *face_tag; /**< Classification of grid's individual - connections (faces). */ - - int number_of_nodes; /**< Number of unique grid vertices. */ - int number_of_nodes_on_pillars; /**< Total number of unique cell - vertices that lie on pillars. */ - double *node_coordinates; /**< Vertex coordinates. Three doubles - (\f$x\f$, \f$y\f$, \f$z\f$) per vertex, - stored sequentially. */ - - int number_of_cells; /**< Number of active grid cells. */ - int *local_cell_index; /**< Deceptively named local-to-global cell - index mapping. */ - }; - - - /** - * Construct a prototypical grid representation from a corner-point - * specification. - * - * Pinched cells will be removed irrespective of any explicit "active" map - * in the geological model input specification. On input, the result - * structure "out" must point to a valid management structure. In other - * words, the result structure must point to a region of memory that is - * typically backed by automatic or allocated (dynamic) storage duration. - * - * @param[in] g Corner-point specification. If "actnum" is NULL, then - * the specification is interpreted as if all cells are - * initially active. - * @param[in] tol Absolute tolerance of node-coincidence. - * @param[in,out] out Minimal grid representation featuring face-to-cell - * neighbourship definition, vertex geometry, face's - * constituent vertices, and local-to-global cell - * mapping. - */ - void process_grdecl(const struct grdecl *g , - double tol, - struct processed_grid *out); - - /** - * Release memory resources acquired in previous grid processing using - * function process_grdecl(). - * - * Note: This function releases the resources associated to the individual - * fields of the processed_grid, but does not free() the structure itself. - * - * @param[in,out] g Prototypical grid representation obtained in an earlier - * call to function process_grdecl(). - */ - void free_processed_grid(struct processed_grid *g); -#ifdef __cplusplus -} -#endif - - -#endif /* OPM_PREPROCESS_HEADER */ - - -/* Local Variables: */ -/* c-basic-offset:4 */ -/* End: */ diff --git a/opm/core/grid/cpgpreprocess/uniquepoints.c b/opm/core/grid/cpgpreprocess/uniquepoints.c deleted file mode 100644 index 17090776..00000000 --- a/opm/core/grid/cpgpreprocess/uniquepoints.c +++ /dev/null @@ -1,383 +0,0 @@ -/*=========================================================================== -// -// File: uniquepoints.c -// -// Created: Fri Jun 19 08:46:05 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//==========================================================================*/ - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include -#include - - -#include "preprocess.h" -#include "uniquepoints.h" - -#define MIN(i,j) (((i) < (j)) ? (i) : (j)) -#define MAX(i,j) (((i) > (j)) ? (i) : (j)) - -/*----------------------------------------------------------------- - Compare function passed to qsortx */ -static int compare(const void *a, const void *b) -{ - const double a0 = *(const double*) a; - const double b0 = *(const double*) b; - - /* { -1, a < b - * compare(a,b) = { 0, a = b - * { 1, a > b */ - return (a0 > b0) - (a0 < b0); -} - -/*----------------------------------------------------------------- - Creat sorted list of z-values in zcorn with actnum==1x */ -static int createSortedList(double *list, int n, int m, - const double *z[], const int *a[]) -{ - int i,j; - double *ptr = list; - for (i=0; i apart in of increasing - doubles. */ -static int uniquify(int n, double *list, double tolerance) -{ - int i; - int pos; - double val; - - assert (!(tolerance < 0.0)); - - if (n<1) return 0; - pos = 0; - val = list[pos++];/* Keep first value */ - - for (i=1; i - away from the last point (list[n-1]), we remove this - second-to-last point as it cannot be distinguished from "final" - point. - */ - list[pos-1] = list[n-1]; - - return pos; -} - - -/*----------------------------------------------------------------- - Along single pillar: */ -static int assignPointNumbers(int begin, - int end, - const double *zlist, - int n, - const double *zcorn, - const int *actnum, - int *plist, - double tolerance) -{ - /* n - number of cells */ - /* zlist - list of len unique z-values */ - /* start - number of unique z-values processed before. */ - - int i, k; - /* All points should now be within tolerance of a listed point. */ - - - const double *z = zcorn; - const int *a = actnum; - int *p = plist; - - k = begin; - *p++ = INT_MIN; /* Padding to ease processing of faults */ - for (i=0; i with k index running faster than i running - faster than j, and Cartesian dimensions , find pointers to the - (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. */ -static void igetvectors(const int dims[3], int i, int j, - const int *field, const int *v[]) -{ - size_t p, start[4]; - - vector_positions(dims, i, j, start); - - for (p = 0; p < 4; p++) { - v[p] = field + start[p]; - } -} - - -/*----------------------------------------------------------------- - Given a vector with k index running faster than i running - faster than j, and Cartesian dimensions , find pointers to the - (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. */ -static void dgetvectors(const int dims[3], int i, int j, - const double *field, const double *v[]) -{ - size_t p, start[4]; - - vector_positions(dims, i, j, start); - - for (p = 0; p < 4; p++) { - v[p] = field + start[p]; - } -} - - -/*----------------------------------------------------------------- - Given a z coordinate, find x and y coordinates on line defined by - coord. Coord points to a vector of 6 doubles [x0,y0,z0,x1,y1,z1]. - */ -static void interpolate_pillar(const double *coord, double *pt) -{ - double a; - - if (fabs(coord[5] - coord[2]) > 0) { - - a = (pt[2] - coord[2]) / (coord[5] - coord[2]); - - } else { - - a = 0; - - } - -#if 0 - pt[0] = coord[0] + a*(coord[3]-coord[0]); - pt[1] = coord[1] + a*(coord[4]-coord[1]); -#else - pt[0] = (1.0 - a)*coord[0] + a*coord[3]; - pt[1] = (1.0 - a)*coord[1] + a*coord[4]; -#endif -} - -/*----------------------------------------------------------------- - Assign point numbers p such that "zlist(p)==zcorn". Assume that - coordinate number is arranged in a sequence such that the natural - index is (k,i,j) */ -int finduniquepoints(const struct grdecl *g, - /* return values: */ - int *plist, /* list of point numbers on - * each pillar*/ - double tolerance, - struct processed_grid *out) - -{ - - const int nx = out->dimensions[0]; - const int ny = out->dimensions[1]; - const int nz = out->dimensions[2]; - const int nc = g->dims[0]*g->dims[1]*g->dims[2]; - - - /* zlist may need extra space temporarily due to simple boundary - * treatement */ - int npillarpoints = 8*(nx+1)*(ny+1)*nz; - int npillars = (nx+1)*(ny+1); - - double *zlist = malloc(npillarpoints*sizeof *zlist); - int *zptr = malloc((npillars+1)*sizeof *zptr); - - - - - int i,j,k; - - int d1[3]; - int len = 0; - double *zout = zlist; - int pos = 0; - double *pt; - const double *z[4]; - const int *a[4]; - int *p; - int pix, cix; - int zix; - - const double *coord = g->coord; - - d1[0] = 2*g->dims[0]; - d1[1] = 2*g->dims[1]; - d1[2] = 2*g->dims[2]; - - out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates)); - - zptr[pos++] = zout - zlist; - - pt = out->node_coordinates; - - /* Loop over pillars, find unique points on each pillar */ - for (j=0; j < g->dims[1]+1; ++j){ - for (i=0; i < g->dims[0]+1; ++i){ - - /* Get positioned pointers for actnum and zcorn data */ - igetvectors(g->dims, i, j, g->actnum, a); - dgetvectors(d1, 2*i, 2*j, g->zcorn, z); - - len = createSortedList( zout, d1[2], 4, z, a); - len = uniquify (len, zout, tolerance); - - /* Assign unique points */ - for (k=0; knumber_of_nodes_on_pillars = zptr[pos-1]; - out->number_of_nodes = zptr[pos-1]; - - /* Loop over all vertical sets of zcorn values, assign point - * numbers */ - p = plist; - for (j=0; j < 2*g->dims[1]; ++j){ - for (i=0; i < 2*g->dims[0]; ++i){ - - /* pillar index */ - pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2); - - /* cell column position */ - cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]); - - /* zcorn column position */ - zix = 2*g->dims[2]*(i+2*g->dims[0]*j); - - if (!assignPointNumbers(zptr[pix], zptr[pix+1], zlist, - 2*g->dims[2], - g->zcorn + zix, g->actnum + cix, - p, tolerance)){ - fprintf(stderr, "Something went wrong in assignPointNumbers"); - return 0; - } - - p += 2 + 2*g->dims[2]; - } - } - - free(zptr); - free(zlist); - - return 1; -} - -/* Local Variables: */ -/* c-basic-offset:4 */ -/* End: */ diff --git a/opm/core/grid/cpgpreprocess/uniquepoints.h b/opm/core/grid/cpgpreprocess/uniquepoints.h deleted file mode 100644 index 08514a41..00000000 --- a/opm/core/grid/cpgpreprocess/uniquepoints.h +++ /dev/null @@ -1,47 +0,0 @@ -/*=========================================================================== -// -// File: uniquepoints.h -// -// Created: Fri Jun 19 08:46:30 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//==========================================================================*/ - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_UNIQUEPOINTS_HEADER -#define OPM_UNIQUEPOINTS_HEADER - -int finduniquepoints(const struct grdecl *g, /* input */ - int *p, /* for each z0 in zcorn, z0 = z[p0] */ - double t, /* tolerance*/ - struct processed_grid *out); - -#endif /* OPM_UNIQUEPOINTS_HEADER */ - -/* Local Variables: */ -/* c-basic-offset:4 */ -/* End: */ diff --git a/opm/core/grid/grid.c b/opm/core/grid/grid.c deleted file mode 100644 index 5d182fd3..00000000 --- a/opm/core/grid/grid.c +++ /dev/null @@ -1,562 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include "config.h" -#include - -#include -#include -#include -#include -#include -#include - - -void -destroy_grid(struct UnstructuredGrid *g) -{ - if (g!=NULL) - { - free(g->face_nodes); - free(g->face_nodepos); - free(g->face_cells); - free(g->cell_facepos); - free(g->cell_faces); - - free(g->node_coordinates); - free(g->face_centroids); - free(g->face_areas); - free(g->face_normals); - free(g->cell_centroids); - free(g->cell_volumes); - - free(g->global_cell); - free(g->cell_facetag); - - free(g->zcorn); - } - - free(g); -} - - -struct UnstructuredGrid * -create_grid_empty(void) -{ - struct UnstructuredGrid *G, g = { 0 }; - - G = malloc(1 * sizeof *G); - - if (G != NULL) { - *G = g; - } - - return G; -} - - -void -attach_zcorn_copy(struct UnstructuredGrid* G , const double * zcorn) -{ - size_t zcorn_elements = G->cartdims[0] * G->cartdims[1] * G->cartdims[2] * 8; - double * new_zcorn = realloc(G->zcorn , zcorn_elements * sizeof * G->zcorn ); - if (new_zcorn) { - G->zcorn = new_zcorn; - memcpy(G->zcorn , zcorn , zcorn_elements * sizeof * G->zcorn ); - } -} - - -struct UnstructuredGrid * -allocate_grid(size_t ndims , - size_t ncells , - size_t nfaces , - size_t nfacenodes, - size_t ncellfaces, - size_t nnodes ) -{ - size_t nel; - struct UnstructuredGrid *G; - - G = create_grid_empty(); - - if (G != NULL) { - /* zcorn cache - only for output. */ - G->zcorn = NULL; - - /* Grid fields ---------------------------------------- */ - G->dimensions = ndims; - G->number_of_cells = ncells; - G->number_of_faces = nfaces; - G->number_of_nodes = nnodes; - - /* Node fields ---------------------------------------- */ - nel = nnodes * ndims; - G->node_coordinates = malloc(nel * sizeof *G->node_coordinates); - - /* Face fields ---------------------------------------- */ - nel = nfacenodes; - G->face_nodes = malloc(nel * sizeof *G->face_nodes); - - nel = nfaces + 1; - G->face_nodepos = malloc(nel * sizeof *G->face_nodepos); - - nel = 2 * nfaces; - G->face_cells = malloc(nel * sizeof *G->face_cells); - - nel = nfaces * ndims; - G->face_centroids = malloc(nel * sizeof *G->face_centroids); - - nel = nfaces * ndims; - G->face_normals = malloc(nel * sizeof *G->face_normals); - - nel = nfaces * 1; - G->face_areas = malloc(nel * sizeof *G->face_areas); - - - /* Cell fields ---------------------------------------- */ - nel = ncellfaces; - G->cell_faces = malloc(nel * sizeof *G->cell_faces); - - G->cell_facetag = malloc(nel * sizeof *G->cell_facetag); - - nel = ncells + 1; - G->cell_facepos = malloc(nel * sizeof *G->cell_facepos); - - nel = ncells * ndims; - G->cell_centroids = malloc(nel * sizeof *G->cell_centroids); - - nel = ncells * 1; - G->cell_volumes = malloc(nel * sizeof *G->cell_volumes); - - if ((G->node_coordinates == NULL) || - (G->face_nodes == NULL) || - (G->face_nodepos == NULL) || - (G->face_cells == NULL) || - (G->face_centroids == NULL) || - (G->face_normals == NULL) || - (G->face_areas == NULL) || - (G->cell_faces == NULL) || - (G->cell_facetag == NULL) || - (G->cell_facepos == NULL) || - (G->cell_centroids == NULL) || - (G->cell_volumes == NULL) ) - { - destroy_grid(G); - G = NULL; - } - } - - return G; -} - - -#define GRID_NMETA 6 -#define GRID_NDIMS 0 -#define GRID_NCELLS 1 -#define GRID_NFACES 2 -#define GRID_NNODES 3 -#define GRID_NFACENODES 4 -#define GRID_NCELLFACES 5 - - -static void -input_error(FILE *fp, const char * const err) -{ - int save_errno = errno; - - if (ferror(fp)) { - fprintf(stderr, "%s: %s\n", err, strerror(save_errno)); - clearerr(fp); - } - else if (feof(fp)) { - fprintf(stderr, "%s: End-of-file\n", err); - } - - errno = save_errno; -} - - -static struct UnstructuredGrid * -allocate_grid_from_file(FILE *fp, int *has_tag, int *has_indexmap) -{ - struct UnstructuredGrid *G; - - int save_errno; - unsigned long tmp; - size_t dimens[GRID_NMETA], i; - - save_errno = errno; - - i = 0; - while ((i < GRID_NMETA) && (fscanf(fp, " %lu", &tmp) == 1)) { - dimens[i] = tmp; - - i += 1; - } - - if (i == GRID_NMETA) { - if (fscanf(fp, "%d %d", has_tag, has_indexmap) == 2) { - G = allocate_grid(dimens[GRID_NDIMS] , - dimens[GRID_NCELLS] , - dimens[GRID_NFACES] , - dimens[GRID_NFACENODES], - dimens[GRID_NCELLFACES], - dimens[GRID_NNODES] ); - - if (G != NULL) { - if (! *has_tag) { - free(G->cell_facetag); - G->cell_facetag = NULL; - } - - if (*has_indexmap) { - G->global_cell = - malloc(dimens[GRID_NCELLS] * sizeof *G->global_cell); - - /* Allocation failure checked elsewhere. */ - } - - G->number_of_cells = (int) dimens[GRID_NCELLS]; - G->number_of_faces = (int) dimens[GRID_NFACES]; - G->number_of_nodes = (int) dimens[GRID_NNODES]; - G->dimensions = (int) dimens[GRID_NDIMS]; - - i = 0; - while ((i < dimens[GRID_NDIMS]) && - (fscanf(fp, "%d", & G->cartdims[ i ]) == 1)) { - i += 1; - } - - if (i < dimens[GRID_NDIMS]) { - input_error(fp, "Unable to read Cartesian dimensions"); - - destroy_grid(G); - G = NULL; - } - else { - /* Account for dimens[GRID_DIMS] < 3 */ - size_t n = (sizeof G->cartdims) / (sizeof G->cartdims[0]); - for (; i < n; i++) { G->cartdims[ i ] = 1; } - } - } - } - else { - input_error(fp, "Unable to read grid predicates"); - - G = NULL; - } - } - else { - input_error(fp, "Unable to read grid dimensions"); - - G = NULL; - } - - errno = save_errno; - - return G; -} - - -static int -read_grid_nodes(FILE *fp, struct UnstructuredGrid *G) -{ - int save_errno; - size_t i, n; - - save_errno = errno; - - n = G->dimensions; - n *= G->number_of_nodes; - - i = 0; - while ((i < n) && - (fscanf(fp, " %lf", & G->node_coordinates[ i ]) == 1)) { - i += 1; - } - - if (i < n) { - input_error(fp, "Unable to read node coordinates"); - } - - errno = save_errno; - - return i == n; -} - - -static int -read_grid_faces(FILE *fp, struct UnstructuredGrid *G) -{ - int save_errno, ok; - size_t nf, nfn, i; - - save_errno = errno; - - nf = G->number_of_faces; - - /* G->face_nodepos */ - i = 0; - while ((i < nf + 1) && - (fscanf(fp, " %d", & G->face_nodepos[ i ]) == 1)) { - i += 1; - } - ok = i == nf + 1; - - if (! ok) { - input_error(fp, "Unable to read node indirection array"); - } - else { - /* G->face_nodes */ - nfn = G->face_nodepos[ nf ]; - - i = 0; - while ((i < nfn) && (fscanf(fp, " %d", & G->face_nodes[ i ]) == 1)) { - i += 1; - } - - ok = i == nfn; - if (! ok) { - input_error(fp, "Unable to read face-nodes"); - } - } - - if (ok) { - /* G->face_cells */ - i = 0; - while ((i < 2 * nf) && (fscanf(fp, " %d", & G->face_cells[ i ]) == 1)) { - i += 1; - } - - ok = i == 2 * nf; - if (! ok) { - input_error(fp, "Unable to read neighbourship"); - } - } - - if (ok) { - /* G->face_areas */ - i = 0; - while ((i < nf) && (fscanf(fp, " %lf", & G->face_areas[ i ]) == 1)) { - i += 1; - } - - ok = i == nf; - if (! ok) { - input_error(fp, "Unable to read face areas"); - } - } - - if (ok) { - /* G->face_centroids */ - size_t n; - - n = G->dimensions; - n *= nf; - - i = 0; - while ((i < n) && (fscanf(fp, " %lf", & G->face_centroids[ i ]) == 1)) { - i += 1; - } - - ok = i == n; - if (! ok) { - input_error(fp, "Unable to read face centroids"); - } - } - - if (ok) { - /* G->face_normals */ - size_t n; - - n = G->dimensions; - n *= nf; - - i = 0; - while ((i < n) && (fscanf(fp, " %lf", & G->face_normals[ i ]) == 1)) { - i += 1; - } - - ok = i == n; - if (! ok) { - input_error(fp, "Unable to read face normals"); - } - } - - errno = save_errno; - - return ok; -} - - -static int -read_grid_cells(FILE *fp, int has_tag, int has_indexmap, - struct UnstructuredGrid *G) -{ - int save_errno, ok; - size_t nc, ncf, i; - - save_errno = errno; - - nc = G->number_of_cells; - - /* G->cell_facepos */ - i = 0; - while ((i < nc + 1) && (fscanf(fp, " %d", & G->cell_facepos[ i ]) == 1)) { - i += 1; - } - ok = i == nc + 1; - - if (! ok) { - input_error(fp, "Unable to read face indirection array"); - } - else { - /* G->cell_faces (and G->cell_facetag if applicable) */ - ncf = G->cell_facepos[ nc ]; - i = 0; - - if (has_tag) { - assert (G->cell_facetag != NULL); - - while ((i < ncf) && - (fscanf(fp, " %d %d", - & G->cell_faces [ i ], - & G->cell_facetag[ i ]) == 2)) { - i += 1; - } - } - else { - while ((i < ncf) && - (fscanf(fp, " %d", & G->cell_faces[ i ]) == 1)) { - i += 1; - } - } - - ok = i == ncf; - if (! ok) { - input_error(fp, "Unable to read cell-faces"); - } - } - - if (ok) { - /* G->global_cell if applicable */ - if (has_indexmap) { - i = 0; - - if (G->global_cell != NULL) { - while ((i < nc) && - (fscanf(fp, " %d", & G->global_cell[ i ]) == 1)) { - i += 1; - } - } - else { - int discard; - - while ((i < nc) && (fscanf(fp, " %d", & discard) == 1)) { - i += 1; - } - } - } - else { - assert (G->global_cell == NULL); - i = nc; - } - - ok = i == nc; - if (! ok) { - input_error(fp, "Unable to read global cellmap"); - } - } - - if (ok) { - /* G->cell_volumes */ - i = 0; - while ((i < nc) && (fscanf(fp, " %lf", & G->cell_volumes[ i ]) == 1)) { - i += 1; - } - - ok = i == nc; - if (! ok) { - input_error(fp, "Unable to read cell volumes"); - } - } - - if (ok) { - /* G->cell_centroids */ - size_t n; - - n = G->dimensions; - n *= nc; - - i = 0; - while ((i < n) && (fscanf(fp, " %lf", & G->cell_centroids[ i ]) == 1)) { - i += 1; - } - - ok = i == n; - if (! ok) { - input_error(fp, "Unable to read cell centroids"); - } - } - - errno = save_errno; - - return ok; -} - - -struct UnstructuredGrid * -read_grid(const char *fname) -{ - struct UnstructuredGrid *G; - FILE *fp; - - int save_errno; - int has_tag, has_indexmap, ok; - - save_errno = errno; - - fp = fopen(fname, "rt"); - if (fp != NULL) { - G = allocate_grid_from_file(fp, & has_tag, & has_indexmap); - - ok = G != NULL; - - if (ok) { ok = read_grid_nodes(fp, G); } - if (ok) { ok = read_grid_faces(fp, G); } - if (ok) { ok = read_grid_cells(fp, has_tag, has_indexmap, G); } - - if (! ok) { - destroy_grid(G); - G = NULL; - } - - fclose(fp); - } - else { - G = NULL; - } - - errno = save_errno; - - return G; -} - - diff --git a/opm/core/grid/grid_equal.cpp b/opm/core/grid/grid_equal.cpp deleted file mode 100644 index f5394c77..00000000 --- a/opm/core/grid/grid_equal.cpp +++ /dev/null @@ -1,76 +0,0 @@ -#include // C string.h to get memcmp() - -#include - -#include - -/* - The grid_equal() function is separated out into a separate file to - ensure that it is compiled with a C++ compiler, so that the - namespace features used in the opm/common/util/numeric/cmp.cpp - implementation compiles. -*/ - - -bool -grid_equal(const struct UnstructuredGrid * grid1 , const struct UnstructuredGrid * grid2) { - if ((grid1->dimensions == grid2->dimensions) && - (grid1->number_of_cells == grid2->number_of_cells) && - (grid1->number_of_faces == grid2->number_of_faces) && - (grid1->number_of_nodes == grid2->number_of_nodes)) { - - // Exact integer comparisons - { - if (memcmp(grid1->face_nodepos , grid2->face_nodepos , (grid1->number_of_faces + 1) * sizeof * grid1->face_nodepos) != 0) - return false; - - if (memcmp(grid1->face_nodes , grid2->face_nodes , grid1->face_nodepos[grid1->number_of_faces] * sizeof * grid1->face_nodes) != 0) - return false; - - if (memcmp(grid1->face_cells , grid2->face_cells , 2 * grid1->number_of_faces * sizeof * grid1->face_cells) != 0) - return false; - - if (memcmp(grid1->cell_faces , grid2->cell_faces , grid1->cell_facepos[grid1->number_of_cells] * sizeof * grid1->cell_faces) != 0) - return false; - - if (memcmp(grid1->cell_facepos , grid2->cell_facepos , (grid1->number_of_cells + 1) * sizeof * grid1->cell_facepos) != 0) - return false; - - if (grid1->global_cell && grid2->global_cell) { - if (memcmp(grid1->global_cell , grid2->global_cell , grid1->number_of_cells * sizeof * grid1->global_cell) != 0) - return false; - } else if (grid1->global_cell != grid2->global_cell) - return false; - - if (grid1->cell_facetag && grid2->cell_facetag) { - if (memcmp(grid1->cell_facetag , grid2->cell_facetag , grid1->cell_facepos[grid1->number_of_cells] * sizeof * grid1->cell_facetag) != 0) - return false; - } else if (grid1->cell_facetag != grid2->cell_facetag) - return false; - } - - - // Floating point comparisons. - { - if (!Opm::cmp::array_equal( grid1->node_coordinates , grid2->node_coordinates , static_cast(grid1->dimensions * grid1->number_of_nodes))) - return false; - - if (!Opm::cmp::array_equal( grid1->face_centroids , grid2->face_centroids , static_cast(grid1->dimensions * grid1->number_of_faces))) - return false; - - if (!Opm::cmp::array_equal( grid1->face_areas , grid2->face_areas , static_cast(grid1->number_of_faces))) - return false; - - if (!Opm::cmp::array_equal( grid1->face_normals , grid2->face_normals , static_cast(grid1->dimensions * grid1->number_of_faces))) - return false; - - if (!Opm::cmp::array_equal( grid1->cell_centroids , grid2->cell_centroids , static_cast(grid1->dimensions * grid1->number_of_cells))) - return false; - - if (!Opm::cmp::array_equal( grid1->cell_volumes , grid2->cell_volumes , static_cast(grid1->number_of_cells))) - return false; - } - return true; - } else - return false; -} diff --git a/opm/core/utility/SparseTable.hpp b/opm/core/utility/SparseTable.hpp deleted file mode 100644 index 60c14ccf..00000000 --- a/opm/core/utility/SparseTable.hpp +++ /dev/null @@ -1,244 +0,0 @@ -//=========================================================================== -// -// File: SparseTable.hpp -// -// Created: Fri Apr 24 09:50:27 2009 -// -// Author(s): Atgeirr F Rasmussen -// -// $Date$ -// -// $Revision$ -// -//=========================================================================== - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_SPARSETABLE_HEADER -#define OPM_SPARSETABLE_HEADER - -#include -#include -#include -#include -#include - -#include - -namespace Opm -{ - - /// A SparseTable stores a table with rows of varying size - /// as efficiently as possible. - /// It is supposed to behave similarly to a vector of vectors. - /// Its behaviour is similar to compressed row sparse matrices. - template - class SparseTable - { - public: - /// Default constructor. Yields an empty SparseTable. - SparseTable() - { - } - - /// A constructor taking all the data for the table and row sizes. - /// \param data_beg The start of the table data. - /// \param data_end One-beyond-end of the table data. - /// \param rowsize_beg The start of the row length data. - /// \param rowsize_end One beyond the end of the row length data. - template - SparseTable(DataIter data_beg, DataIter data_end, - IntegerIter rowsize_beg, IntegerIter rowsize_end) - : data_(data_beg, data_end) - { - setRowStartsFromSizes(rowsize_beg, rowsize_end); - } - - - /// Sets the table to contain the given data, organized into - /// rows as indicated by the given row sizes. - /// \param data_beg The start of the table data. - /// \param data_end One-beyond-end of the table data. - /// \param rowsize_beg The start of the row length data. - /// \param rowsize_end One beyond the end of the row length data. - template - void assign(DataIter data_beg, DataIter data_end, - IntegerIter rowsize_beg, IntegerIter rowsize_end) - { - data_.assign(data_beg, data_end); - setRowStartsFromSizes(rowsize_beg, rowsize_end); - } - - - /// Request storage for table of given size. - /// \param rowsize_beg Start of row size data. - /// \param rowsize_end One beyond end of row size data. - template - void allocate(IntegerIter rowsize_beg, IntegerIter rowsize_end) - { - typedef typename std::vector::size_type sz_t; - - sz_t ndata = std::accumulate(rowsize_beg, rowsize_end, sz_t(0)); - data_.resize(ndata); - setRowStartsFromSizes(rowsize_beg, rowsize_end); - } - - - /// Appends a row to the table. - template - void appendRow(DataIter row_beg, DataIter row_end) - { - data_.insert(data_.end(), row_beg, row_end); - if (row_start_.empty()) { - row_start_.reserve(2); - row_start_.push_back(0); - } - row_start_.push_back(data_.size()); - } - - /// True if the table contains no rows. - bool empty() const - { - return row_start_.empty(); - } - - /// Returns the number of rows in the table. - int size() const - { - return empty() ? 0 : row_start_.size() - 1; - } - - /// Allocate storage for table of expected size - void reserve(int exptd_nrows, int exptd_ndata) - { - row_start_.reserve(exptd_nrows + 1); - data_.reserve(exptd_ndata); - } - - /// Swap contents for other SparseTable - void swap(SparseTable& other) - { - row_start_.swap(other.row_start_); - data_.swap(other.data_); - } - - /// Returns the number of data elements. - int dataSize() const - { - return data_.size(); - } - - /// Returns the size of a table row. - int rowSize(int row) const - { -#ifndef NDEBUG - OPM_ERROR_IF(row < 0 || row >= size(), "Row index " << row << " is out of range"); -#endif - return row_start_[row + 1] - row_start_[row]; - } - - /// Makes the table empty(). - void clear() - { - data_.clear(); - row_start_.clear(); - } - - /// Defining the row type, returned by operator[]. - typedef boost::iterator_range row_type; - typedef boost::iterator_range mutable_row_type; - - /// Returns a row of the table. - row_type operator[](int row) const - { - assert(row >= 0 && row < size()); - const T* start_ptr = data_.empty() ? 0 : &data_[0]; - return row_type(start_ptr + row_start_[row], start_ptr + row_start_[row + 1]); - } - - /// Returns a mutable row of the table. - mutable_row_type operator[](int row) - { - assert(row >= 0 && row < size()); - T* start_ptr = data_.empty() ? 0 : &data_[0]; - return mutable_row_type(start_ptr + row_start_[row], start_ptr + row_start_[row + 1]); - } - - /// Equality. - bool operator==(const SparseTable& other) const - { - return data_ == other.data_ && row_start_ == other.row_start_; - } - - template - void print(std::basic_ostream& os) const - { - os << "Number of rows: " << size() << '\n'; - - os << "Row starts = ["; - std::copy(row_start_.begin(), row_start_.end(), - std::ostream_iterator(os, " ")); - os << "\b]\n"; - - os << "Data values = ["; - std::copy(data_.begin(), data_.end(), - std::ostream_iterator(os, " ")); - os << "\b]\n"; - } - const T data(int i)const { - return data_[i]; - } - - private: - std::vector data_; - // Like in the compressed row sparse matrix format, - // row_start_.size() is equal to the number of rows + 1. - std::vector row_start_; - - template - void setRowStartsFromSizes(IntegerIter rowsize_beg, IntegerIter rowsize_end) - { - // Since we do not store the row sizes, but cumulative row sizes, - // we have to create the cumulative ones. - int num_rows = rowsize_end - rowsize_beg; - if (num_rows < 1) { - OPM_THROW(std::runtime_error, "Must have at least one row. Got " << num_rows << " rows."); - } -#ifndef NDEBUG - if (*std::min_element(rowsize_beg, rowsize_end) < 0) { - OPM_THROW(std::runtime_error, "All row sizes must be at least 0."); - } -#endif - row_start_.resize(num_rows + 1); - row_start_[0] = 0; - std::partial_sum(rowsize_beg, rowsize_end, row_start_.begin() + 1); - // Check that data_ and row_start_ match. - if (int(data_.size()) != row_start_.back()) { - OPM_THROW(std::runtime_error, "End of row start indices different from data size."); - } - - } - }; - -} // namespace Opm - - -#endif // OPM_SPARSETABLE_HEADER diff --git a/opm/core/utility/StopWatch.cpp b/opm/core/utility/StopWatch.cpp deleted file mode 100644 index 3ed60c7f..00000000 --- a/opm/core/utility/StopWatch.cpp +++ /dev/null @@ -1,106 +0,0 @@ -//=========================================================================== -// -// File: StopWatch.cpp -// -// Created: Thu Jul 2 23:04:51 2009 -// -// Author(s): Atgeirr F Rasmussen -// -// $Date$ -// -// $Revision$ -// -//=========================================================================== - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#if HAVE_CONFIG_H -#include "config.h" -#endif -#include -#include -#include - -namespace Opm -{ - - namespace time - { - - StopWatch::StopWatch() - : state_(UnStarted) - { - } - - - void StopWatch::start() - { - start_time_ = boost::posix_time::microsec_clock::local_time(); - last_time_ = start_time_; - state_ = Running; - } - - void StopWatch::stop() - { - if (state_ != Running) { - OPM_THROW(std::runtime_error, "Called stop() on a StopWatch that was not running."); - } - stop_time_ = boost::posix_time::microsec_clock::local_time(); - state_ = Stopped; - } - - double StopWatch::secsSinceLast() - { - boost::posix_time::ptime run_time; - if (state_ == Running) { - run_time = boost::posix_time::microsec_clock::local_time(); - } else if (state_ == Stopped) { - run_time = stop_time_; - } else { - assert(state_ == UnStarted); - OPM_THROW(std::runtime_error, "Called secsSinceLast() on a StopWatch that had not been started."); - } - boost::posix_time::time_duration dur = run_time - last_time_; - last_time_ = run_time; - return double(dur.total_microseconds())/1000000.0; - } - - double StopWatch::secsSinceStart() - { - boost::posix_time::ptime run_time; - if (state_ == Running) { - run_time = boost::posix_time::microsec_clock::local_time(); - } else if (state_ == Stopped) { - run_time = stop_time_; - } else { - assert(state_ == UnStarted); - OPM_THROW(std::runtime_error, "Called secsSinceStart() on a StopWatch that had not been started."); - } - boost::posix_time::time_duration dur = run_time - start_time_; - return double(dur.total_microseconds())/1000000.0; - } - - } // namespace time - -} // namespace Opm - - - diff --git a/opm/core/utility/StopWatch.hpp b/opm/core/utility/StopWatch.hpp deleted file mode 100644 index ec3099e9..00000000 --- a/opm/core/utility/StopWatch.hpp +++ /dev/null @@ -1,81 +0,0 @@ -//=========================================================================== -// -// File: StopWatch.hpp -// -// Created: Thu Jul 2 23:04:17 2009 -// -// Author(s): Atgeirr F Rasmussen -// -// $Date$ -// -// $Revision$ -// -//=========================================================================== - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_STOPWATCH_HEADER -#define OPM_STOPWATCH_HEADER - -#include - -namespace Opm -{ - - namespace time - { - - class StopWatch - { - public: - /// Default constructor. Before the StopWatch is start()-ed, - /// it is an error to call anything other than start(). - StopWatch(); - - /// Starts the StopWatch. It is always legal to call - /// start(), even if not stop()-ped. - void start(); - /// Stops the StopWatch. The watch no longer runs, until - /// restarted by a call to start(). - void stop(); - - /// \return the number of running seconds that have passed - /// since last call to start(), secsSinceLast() or - /// secsSinceStart() - double secsSinceLast(); - /// \return the number of running seconds that have passed - /// since last call to start(). - double secsSinceStart(); - - private: - enum StopWatchState { UnStarted, Running, Stopped }; - - StopWatchState state_; - boost::posix_time::ptime start_time_; - boost::posix_time::ptime last_time_; - boost::posix_time::ptime stop_time_; - }; - - } // namespace time - -} // namespace Opm - -#endif // OPM_STOPWATCH_HEADER diff --git a/tests/test_sparsetable.cpp b/tests/test_sparsetable.cpp deleted file mode 100644 index 6aa461f3..00000000 --- a/tests/test_sparsetable.cpp +++ /dev/null @@ -1,132 +0,0 @@ -//=========================================================================== -// -// File: test_sparsetable.cpp -// -// Created: Thu May 28 10:01:46 2009 -// -// Author(s): Atgeirr F Rasmussen -// Bård Skaflestad -// -// $Date$ -// -// $Revision$ -// -//=========================================================================== - -/* - Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. - Copyright 2009, 2010 Statoil ASA. - - This file is part of The Open Reservoir Simulator Project (OpenRS). - - OpenRS is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenRS is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OpenRS. If not, see . -*/ - -#include - -#if HAVE_DYNAMIC_BOOST_TEST -#define BOOST_TEST_DYN_LINK -#endif -#define NVERBOSE // to suppress our messages when throwing - -#define BOOST_TEST_MODULE SparseTableTest -#include - -#include - -using namespace Opm; - -BOOST_AUTO_TEST_CASE(construction_and_queries) -{ - const SparseTable st1; - BOOST_CHECK(st1.empty()); - BOOST_CHECK_EQUAL(st1.size(), 0); - BOOST_CHECK_EQUAL(st1.dataSize(), 0); - - // This should be getting us a table like this: - // ---------------- - // 0 - // - // 1 2 - // 3 4 5 6 - // 7 8 9 - // ---------------- - const int num_elem = 10; - const int elem[num_elem] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }; - const int num_rows = 5; - const int rowsizes[num_rows] = { 1, 0, 2, 4, 3 }; - const SparseTable st2(elem, elem + num_elem, rowsizes, rowsizes + num_rows); - BOOST_CHECK(!st2.empty()); - BOOST_CHECK_EQUAL(st2.size(), num_rows); - BOOST_CHECK_EQUAL(st2.dataSize(), num_elem); - BOOST_CHECK_EQUAL(st2[0][0], 0); - BOOST_CHECK_EQUAL(st2.rowSize(0), 1); - BOOST_CHECK(st2[1].empty()); - BOOST_CHECK_EQUAL(st2.rowSize(1), 0); - BOOST_CHECK_EQUAL(st2[3][1], 4); - BOOST_CHECK_EQUAL(st2[4][2], 9); - BOOST_CHECK(int(st2[4].size()) == rowsizes[4]); - const SparseTable st2_again(elem, elem + num_elem, rowsizes, rowsizes + num_rows); - BOOST_CHECK(st2 == st2_again); - SparseTable st2_byassign; - st2_byassign.assign(elem, elem + num_elem, rowsizes, rowsizes + num_rows); - BOOST_CHECK(st2 == st2_byassign); - const int last_row_size = rowsizes[num_rows - 1]; - SparseTable st2_append(elem, elem + num_elem - last_row_size, rowsizes, rowsizes + num_rows - 1); - BOOST_CHECK_EQUAL(st2_append.dataSize(), num_elem - last_row_size); - st2_append.appendRow(elem + num_elem - last_row_size, elem + num_elem); - BOOST_CHECK(st2 == st2_append); - SparseTable st2_append2; - st2_append2.appendRow(elem, elem + 1); - st2_append2.appendRow(elem + 1, elem + 1); - st2_append2.appendRow(elem + 1, elem + 3); - st2_append2.appendRow(elem + 3, elem + 7); - st2_append2.appendRow(elem + 7, elem + 10); - BOOST_CHECK(st2 == st2_append2); - st2_append2.clear(); - SparseTable st_empty; - BOOST_CHECK(st2_append2 == st_empty); - - SparseTable st2_allocate; - st2_allocate.allocate(rowsizes, rowsizes + num_rows); - BOOST_CHECK_EQUAL(st2_allocate.size(), num_rows); - BOOST_CHECK_EQUAL(st2_allocate.dataSize(), num_elem); - int s = 0; - for (int i = 0; i < num_rows; ++i) { - SparseTable::mutable_row_type row = st2_allocate[i]; - for (int j = 0; j < rowsizes[i]; ++j, ++s) - row[j] = elem[s]; - } - BOOST_CHECK(st2 == st2_allocate); - - // One element too few. - BOOST_CHECK_THROW(const SparseTable st3(elem, elem + num_elem - 1, rowsizes, rowsizes + num_rows), std::exception); - - // A few elements too many. - BOOST_CHECK_THROW(const SparseTable st4(elem, elem + num_elem, rowsizes, rowsizes + num_rows - 1), std::exception); - - // Need at least one row. - BOOST_CHECK_THROW(const SparseTable st5(elem, elem + num_elem, rowsizes, rowsizes), std::exception); - - // Tests that only run in debug mode. -#ifndef NDEBUG - // Do not ask for wrong row numbers. - BOOST_CHECK_THROW(st1.rowSize(0), std::exception); - BOOST_CHECK_THROW(st2.rowSize(-1), std::exception); - BOOST_CHECK_THROW(st2.rowSize(st2.size()), std::exception); - // No negative row sizes. - const int err_rs[num_rows] = { 1, 0, -1, 7, 3 }; - BOOST_CHECK_THROW(const SparseTable st6(elem, elem + num_elem, err_rs, err_rs + num_rows), std::exception); -#endif -}