OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Remove files moved to opm-grid.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2016-11-29 13:25:22 +01:00
parent 647e005297
commit c7cc24385a
31 changed files with 0 additions and 7516 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
CMakeLists_files.cmake
View File

@ -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

332
opm/core/grid.h
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GRID_HEADER_INCLUDED
#define OPM_GRID_HEADER_INCLUDED
#include <stddef.h>
#include <stdbool.h>
/**
* \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
<code>NULL</code> as appropriate. <code>NULL</code> 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
<code>global_cell</code> allocated, but not filled with meaningful
values. <code>NULL</code> 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 */

260
opm/core/grid/CellQuadrature.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_CELLQUADRATURE_HEADER_INCLUDED
#define OPM_CELLQUADRATURE_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/common/ErrorMacros.hpp>
#include <algorithm>
#include <cmath>
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

122
opm/core/grid/ColumnExtract.hpp
View File

@ -1,122 +0,0 @@
#include <opm/core/grid.h>
#include <vector>
#include <map>
#include <algorithm>
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<std::vector<int> >& columns )
{
const int* dims = grid.cartdims;
// Keeps track of column_index ---> index of vector
std::map<int, int> 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<int, int>::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<int>());
}
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<int> > 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<int>());
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

189
opm/core/grid/FaceQuadrature.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FACEQUADRATURE_HEADER_INCLUDED
#define OPM_FACEQUADRATURE_HEADER_INCLUDED
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/grid.h>
#include <cmath>
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

177
opm/core/grid/GridHelpers.cpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/grid/GridHelpers.hpp>
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<grid.cell_facepos[cell_index+1]; ++i) {
// Traverse the vertices associated with each face
assert(grid.face_nodepos[grid.cell_faces[i]+1] - grid.face_nodepos[grid.cell_faces[i]] == nv/2);
for (int j=grid.face_nodepos[grid.cell_faces[i]]; j<grid.face_nodepos[grid.cell_faces[i]+1]; ++j) {
zz += (grid.node_coordinates+nd*(grid.face_nodes[j]))[nd-1];
}
}
return zz/nv;
}
double cellCentroidCoordinate(const UnstructuredGrid& grid, int cell_index,
int coordinate)
{
return grid.cell_centroids[grid.dimensions*cell_index+coordinate];
}
const double*
cellCentroid(const UnstructuredGrid& grid, int cell_index)
{
return grid.cell_centroids+(cell_index*grid.dimensions);
}
const double* beginCellVolumes(const UnstructuredGrid& grid)
{
return grid.cell_volumes;
}
const double* endCellVolumes(const UnstructuredGrid& grid)
{
return grid.cell_volumes+numCells(grid);
}
const double* beginFaceCentroids(const UnstructuredGrid& grid)
{
return grid.face_centroids;
}
const double* faceCentroid(const UnstructuredGrid& grid, int face_index)
{
return grid.face_centroids+face_index*grid.dimensions;
}
const double* faceNormal(const UnstructuredGrid& grid, int face_index)
{
return grid.face_normals+face_index*grid.dimensions;
}
double faceArea(const UnstructuredGrid& grid, int face_index)
{
return grid.face_areas[face_index];
}
int faceTag(const UnstructuredGrid& grid,
boost::iterator_range<const int*>::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<UnstructuredGrid>::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<size_t>(dims[0])) &&
(inputGrid.getNY( ) == static_cast<size_t>(dims[1])) &&
(inputGrid.getNZ( ) == static_cast<size_t>(dims[2]))) {
std::vector<int> 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");
}
}
}
}

381
opm/core/grid/GridHelpers.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
#define OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/range/iterator_range.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
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<const int*>
{
public:
typedef boost::iterator_range<const int*> 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<const int*> 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<class G>
struct CellCentroidTraits
{
};
template<>
struct CellCentroidTraits<UnstructuredGrid>
{
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<UnstructuredGrid>::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<class T>
struct CellVolumeIteratorTraits
{
};
template<>
struct CellVolumeIteratorTraits<UnstructuredGrid>
{
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<class G>
struct FaceCentroidTraits
{
};
template<>
struct FaceCentroidTraits<UnstructuredGrid>
{
typedef const double* IteratorType;
typedef const double* ValueType;
};
/// \brief Get an iterator over the face centroids positioned at the first cell.
FaceCentroidTraits<UnstructuredGrid>::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<UnstructuredGrid>::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 int*>::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<class T>
struct Cell2FacesTraits
{
};
template<>
struct Cell2FacesTraits<UnstructuredGrid>
{
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<class T>
struct Face2VerticesTraits
{
};
template<>
struct Face2VerticesTraits<UnstructuredGrid>
{
typedef SparseTableView Type;
};
/// \brief Get the cell to faces mapping of a grid.
Cell2FacesTraits<UnstructuredGrid>::Type
cell2Faces(const UnstructuredGrid& grid);
/// \brief Get the face to vertices mapping of a grid.
Face2VerticesTraits<UnstructuredGrid>::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<class T>
struct FaceCellTraits
{};
template<>
struct FaceCellTraits<UnstructuredGrid>
{
typedef FaceCellsProxy Type;
};
/// \brief Get the face to cell mapping of a grid.
FaceCellTraits<UnstructuredGrid>::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<class T>
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<class T>
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<class T>
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<class T>
double getCoordinate(T t, int i)
{
return (*t)[i];
}
} // end namespace UGGridHelpers
} // end namespace OPM
#endif

233
opm/core/grid/GridManager.cpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/grid/MinpvProcessor.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/DeckItem.hpp>
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <array>
#include <algorithm>
#include <numeric>
namespace Opm
{
/// Construct a 3d corner-point grid from a deck.
GridManager::GridManager(const Opm::EclipseGrid& inputGrid)
: ug_(0)
{
initFromEclipseGrid(inputGrid, std::vector<double>());
}
GridManager::GridManager(const Opm::EclipseGrid& inputGrid,
const std::vector<double>& 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<double>& poreVolumes)
{
struct grdecl g;
int cells_modified = 0;
std::vector<int> actnum;
std::vector<double> coord;
std::vector<double> zcorn;
std::vector<double> 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<double>& zcorn = deck.getKeyword("ZCORN").getSIDoubleData();
const std::vector<double>& coord = deck.getKeyword("COORD").getSIDoubleData();
const int* actnum = NULL;
if (deck.hasKeyword("ACTNUM")) {
actnum = &(deck.getKeyword("ACTNUM").getIntData()[0]);
}
std::array<int, 3> 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<double*>(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

98
opm/core/grid/GridManager.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GRIDMANAGER_HEADER_INCLUDED
#define OPM_GRIDMANAGER_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <string>
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<double>& 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<double>& poreVolumes);
// The managed UnstructuredGrid.
UnstructuredGrid* ug_;
};
} // namespace Opm
#endif // OPM_GRIDMANAGER_HEADER_INCLUDED

133
opm/core/grid/GridUtilities.cpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <opm/core/grid/GridUtilities.hpp>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/math/constants/constants.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <set>
#include <vector>
#include <cmath>
#include <algorithm>
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<int> 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<std::set<int>> 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<int> 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<int> 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<int>& 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<double, int> AngleAndPos;
std::vector<AngleAndPos> angle_and_pos;
std::vector<int> 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<double>() - 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

42
opm/core/grid/GridUtilities.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GRIDUTILITIES_HEADER_INCLUDED
#define OPM_GRIDUTILITIES_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/utility/SparseTable.hpp>
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<int> 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<int>& nb);
} // namespace Opm
#endif // OPM_GRIDUTILITIES_HEADER_INCLUDED

176
opm/core/grid/MinpvProcessor.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MINPVPROCESSOR_HEADER_INCLUDED
#define OPM_MINPVPROCESSOR_HEADER_INCLUDED
#include <opm/common/ErrorMacros.hpp>
#include <array>
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<double>& pv,
const double minpv,
const std::vector<int>& actnum,
const bool mergeMinPVCells,
double* zcorn) const;
private:
std::array<int,8> cornerIndices(const int i, const int j, const int k) const;
std::array<double, 8> 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<double, 8>& cellz, double* z) const;
std::array<int, 3> dims_;
std::array<int, 3> 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<double>& pv,
const double minpv,
const std::vector<int>& 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<double, 8> 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<double, 8> 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<int,8> 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<int, 8> 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<double, 8> MinpvProcessor::getCellZcorn(const int i, const int j, const int k, const double* z) const
{
const std::array<int, 8> ixs = cornerIndices(i, j, k);
std::array<double, 8> 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<double, 8>& cellz, double* z) const
{
const std::array<int, 8> ixs = cornerIndices(i, j, k);
for (int count = 0; count < 8; ++count) {
z[ixs[count]] = cellz[count];
}
}
} // namespace Opm
#endif // OPM_MINPVPROCESSOR_HEADER_INCLUDED

488
opm/core/grid/PinchProcessor.hpp
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PINCHPROCESSOR_HEADER_INCLUDED
#define OPM_PINCHPROCESSOR_HEADER_INCLUDED
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/PinchMode.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <array>
#include <iostream>
#include <algorithm>
#include <unordered_map>
#include <limits>
namespace Opm
{
template <class Grid>
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<double>& htrans,
const std::vector<int>& actnum,
const std::vector<double>& multz,
const std::vector<double>& pv,
NNC& nnc);
private:
double minpvValue_;
double thickness_;
PinchMode::ModeEnum transMode_;
PinchMode::ModeEnum multzMode_;
/// Mark minpved cells.
std::vector<int> getMinpvCells_(const std::vector<int>& actnum,
const std::vector<double>& 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<std::vector<int> >
getPinchoutsColumn_(const Grid& grid,
const std::vector<int>& actnum,
const std::vector<double>& pv);
/// Get global cell index.
int getGlobalIndex_(const int i, const int j, const int k, const int* dims);
/// Get cartesian index.
std::array<int, 3> getCartIndex_(const int idx,
const int* dims);
/// Compute transmissibility for nnc.
std::vector<double> transCompute_(const Grid& grid,
const std::vector<double>& htrans,
const std::vector<int>& pinCells,
const std::vector<int>& pinFaces);
/// Get map between half-trans index and the pair of face index and cell index.
std::vector<int> 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<double>& htrans,
const std::vector<int>& actnum,
const std::vector<double>& multz,
const std::vector<double>& pv,
NNC& nnc);
/// Item 5 in PINCH keyword.
std::unordered_multimap<int, double> multzOptions_(const Grid& grid,
const std::vector<int>& pinCells,
const std::vector<int>& pinFaces,
const std::vector<double>& multz,
const std::vector<std::vector<int> >& seg);
/// Apply multz vector to face transmissibility.
void applyMultz_(std::vector<double>& trans,
const std::unordered_multimap<int, double>& multzmap);
};
template <class Grid>
inline PinchProcessor<Grid>::PinchProcessor(const double minpv,
const double thickness,
const PinchMode::ModeEnum transMode,
const PinchMode::ModeEnum multzMode)
{
minpvValue_ = minpv;
thickness_ = thickness;
transMode_ = transMode;
multzMode_ = multzMode;
}
template <class Grid>
inline int PinchProcessor<Grid>::getGlobalIndex_(const int i, const int j, const int k, const int* dims)
{
return i + dims[0] * (j + dims[1] * k);
}
template <class Grid>
inline std::array<int, 3> PinchProcessor<Grid>::getCartIndex_(const int idx,
const int* dims)
{
std::array<int, 3> ijk;
ijk[0] = (idx % dims[0]);
ijk[1] = ((idx / dims[0]) % dims[1]);
ijk[2] = ((idx / dims[0]) / dims[1]);
return ijk;
}
template<class Grid>
inline int PinchProcessor<Grid>::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<< "("<<ijk1[0]<<","<<ijk1[1]<<","<<ijk1[2]<<")"
<< " and " << cellIdx2<<"("<<ijk2[0]<<","<<ijk2[1]<<","<<ijk2[2]<<")");
}
return commonFace;
}
template<class Grid>
inline int PinchProcessor<Grid>::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<class Grid>
inline std::vector<int> PinchProcessor<Grid>::getMinpvCells_(const std::vector<int>& actnum,
const std::vector<double>& pv)
{
std::vector<int> minpvCells(pv.size(), 0);
for (int idx = 0; idx < static_cast<int>(pv.size()); ++idx) {
if (actnum[idx]) {
if (pv[idx] < minpvValue_) {
minpvCells[idx] = 1;
}
}
}
return minpvCells;
}
template<class Grid>
inline std::vector<int> PinchProcessor<Grid>::getHfIdxMap_(const Grid& grid)
{
std::vector<int> 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<class Grid>
inline int PinchProcessor<Grid>::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<class Grid>
inline std::vector<double> PinchProcessor<Grid>::transCompute_(const Grid& grid,
const std::vector<double>& htrans,
const std::vector<int>& pinCells,
const std::vector<int>& pinFaces)
{
const int nc = Opm::UgGridHelpers::numCells(grid);
const int nf = Opm::UgGridHelpers::numFaces(grid);
std::vector<double> 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<class Grid>
inline std::vector<std::vector<int>> PinchProcessor<Grid>::getPinchoutsColumn_(const Grid& grid,
const std::vector<int>& actnum,
const std::vector<double>& pv)
{
const int* dims = Opm::UgGridHelpers::cartDims(grid);
std::vector<int> minpvCells = getMinpvCells_(actnum, pv);
std::vector<std::vector<int>> 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<int> 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<class Grid>
inline void PinchProcessor<Grid>::transTopbot_(const Grid& grid,
const std::vector<double>& htrans,
const std::vector<int>& actnum,
const std::vector<double>& multz,
const std::vector<double>& pv,
NNC& nnc)
{
const int* dims = Opm::UgGridHelpers::cartDims(grid);
std::vector<int> pinFaces;
std::vector<int> pinCells;
std::vector<std::vector<int> > newSeg;
auto minpvSeg = getPinchoutsColumn_(grid, actnum, pv);
for (auto& seg : minpvSeg) {
std::array<int, 3> ijk1 = getCartIndex_(seg.front(), dims);
std::array<int, 3> 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<int>(pinCells.size())/2; ++i) {
nnc.addNNC(static_cast<int>(pinCells[2*i]), static_cast<int>(pinCells[2*i+1]), faceTrans[pinFaces[2*i]]);
}
}
template<class Grid>
inline std::unordered_multimap<int, double> PinchProcessor<Grid>::multzOptions_(const Grid& grid,
const std::vector<int>& pinCells,
const std::vector<int>& pinFaces,
const std::vector<double>& multz,
const std::vector<std::vector<int> >& segs)
{
std::unordered_multimap<int, double> multzmap;
if (multzMode_ == PinchMode::ModeEnum::TOP) {
for (int i = 0; i < static_cast<int>(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<double>::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<class Grid>
inline void PinchProcessor<Grid>::applyMultz_(std::vector<double>& trans,
const std::unordered_multimap<int, double>& multzmap)
{
for (auto& x : multzmap) {
trans[x.first] *= x.second;
}
}
template<class Grid>
inline void PinchProcessor<Grid>::process(const Grid& grid,
const std::vector<double>& htrans,
const std::vector<int>& actnum,
const std::vector<double>& multz,
const std::vector<double>& pv,
NNC& nnc)
{
transTopbot_(grid, htrans, actnum, multz, pv, nnc);
}
} // namespace Opm
#endif // OPM_PINCHPROCESSOR_HEADER_INCLUDED

708
opm/core/grid/cart_grid.c
View File

@ -1,708 +0,0 @@
/*===========================================================================
//
// File: cart_grid.c
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <stdlib.h>
#include <stdio.h>
#include <opm/core/grid.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/grid/cart_grid.h>
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; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*cfaces++ = i+ Nx*(j+ ny* k );
*cfaces++ = i+1+Nx*(j+ ny* k );
*cfaces++ = i+ nx*(j+ Ny* k ) +nxf;
*cfaces++ = i+ nx*(j+1+Ny* k ) +nxf;
*cfaces++ = i+ nx*(j+ ny* k ) +nxf+nyf;
*cfaces++ = i+ nx*(j+ ny*(k+1)) +nxf+nyf;
cfacepos[1] = cfacepos[0]+6;
++cfacepos;
}
}
}
for (k = 0; k < nx * ny * nz; ++k) {
for (i = 0; i < 6; ++i) {
G->cell_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; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
*fnodes++ = i+Nx*(j + Ny * k );
*fnodes++ = i+Nx*(j+1 + Ny * k );
*fnodes++ = i+Nx*(j+1 + Ny *(k+1));
*fnodes++ = i+Nx*(j + Ny *(k+1));
fnodepos[1] = fnodepos[0] + 4;
++fnodepos;
if (i==0) {
*fcells++ = -1;
*fcells++ = i+nx*(j+ny*k);
}
else if (i == nx) {
*fcells++ = i-1+nx*(j+ny*k);
*fcells++ = -1;
}
else {
*fcells++ = i-1 + nx*(j+ny*k);
*fcells++ = i + nx*(j+ny*k);
}
}
}
}
/* Faces with y-normal */
for (k=0; k<nz; ++k) {
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
*fnodes++ = i+ Nx*(j + Ny * k );
*fnodes++ = i + Nx*(j + Ny *(k+1));
*fnodes++ = i+1 + Nx*(j + Ny *(k+1));
*fnodes++ = i+1 + Nx*(j + Ny * k );
fnodepos[1] = fnodepos[0] + 4;
++fnodepos;
if (j==0) {
*fcells++ = -1;
*fcells++ = i+nx*(j+ny*k);
}
else if (j == ny) {
*fcells++ = i+nx*(j-1+ny*k);
*fcells++ = -1;
}
else {
*fcells++ = i+nx*(j-1+ny*k);
*fcells++ = i+nx*(j+ny*k);
}
}
}
}
/* Faces with z-normal */
for (k=0; k<nz+1; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*fnodes++ = i+ Nx*(j + Ny * k);
*fnodes++ = i+1 + Nx*(j + Ny * k);
*fnodes++ = i+1 + Nx*(j+1 + Ny * k);
*fnodes++ = i+ Nx*(j+1 + Ny * k);
fnodepos[1] = fnodepos[0] + 4;
++fnodepos;
if (k==0) {
*fcells++ = -1;
*fcells++ = i+nx*(j+ny*k);
}
else if (k == nz) {
*fcells++ = i+nx*(j+ny*(k-1));
*fcells++ = -1;
}
else {
*fcells++ = i+nx*(j+ny*(k-1));
*fcells++ = i+nx*(j+ny*k);
}
}
}
}
}
/* --------------------------------------------------------------------- */
static void
fill_cart_geometry_3d(struct UnstructuredGrid *G,
const double *x,
const double *y,
const double *z)
{
int nx, ny, nz;
int i,j,k;
double dx, dy, dz;
double *coord, *ccentroids, *cvolumes;
double *fnormals, *fcentroids, *fareas;
nx = G->cartdims[0];
ny = G->cartdims[1];
nz = G->cartdims[2];
ccentroids = G->cell_centroids;
cvolumes = G->cell_volumes;
for (k=0; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*ccentroids++ = (x[i] + x[i + 1]) / 2.0;
*ccentroids++ = (y[j] + y[j + 1]) / 2.0;
*ccentroids++ = (z[k] + z[k + 1]) / 2.0;
dx = x[i + 1] - x[i];
dy = y[j + 1] - y[j];
dz = z[k + 1] - z[k];
*cvolumes++ = dx * dy * dz;
}
}
}
fnormals = G->face_normals;
fcentroids = G->face_centroids;
fareas = G->face_areas;
/* Faces with x-normal */
for (k=0; k<nz; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
dy = y[j + 1] - y[j];
dz = z[k + 1] - z[k];
*fnormals++ = dy * dz;
*fnormals++ = 0;
*fnormals++ = 0;
*fcentroids++ = x[i];
*fcentroids++ = (y[j] + y[j + 1]) / 2.0;
*fcentroids++ = (z[k] + z[k + 1]) / 2.0;
*fareas++ = dy * dz;
}
}
}
/* Faces with y-normal */
for (k=0; k<nz; ++k) {
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
dx = x[i + 1] - x[i];
dz = z[k + 1] - z[k];
*fnormals++ = 0;
*fnormals++ = dx * dz;
*fnormals++ = 0;
*fcentroids++ = (x[i] + x[i + 1]) / 2.0;
*fcentroids++ = y[j];
*fcentroids++ = (z[k] + z[k + 1]) / 2.0;
*fareas++ = dx * dz;
}
}
}
/* Faces with z-normal */
for (k=0; k<nz+1; ++k) {
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
dx = x[i + 1] - x[i];
dy = y[j + 1] - y[j];
*fnormals++ = 0;
*fnormals++ = 0;
*fnormals++ = dx * dy;
*fcentroids++ = (x[i] + x[i + 1]) / 2.0;
*fcentroids++ = (y[j] + y[j + 1]) / 2.0;
*fcentroids++ = z[k];
*fareas++ = dx * dy;
}
}
}
coord = G->node_coordinates;
for (k=0; k<nz+1; ++k) {
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx+1; ++i) {
*coord++ = x[i];
*coord++ = y[j];
*coord++ = z[k];
}
}
}
}
/* ---------------------------------------------------------------------- */
static void
fill_layered_geometry_3d(struct UnstructuredGrid *G,
const double *x,
const double *y,
const double *z,
const double *depthz)
{
int i , j , k ;
int nx, ny, nz;
const double *depth;
double *coord;
nx = G->cartdims[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; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*cfaces++ = i+ Nx*j;
*cfaces++ = i+ nx*j +nxf;
*cfaces++ = i+1+Nx*j;
*cfaces++ = i+ nx*(j+1)+nxf;
cfacepos[1] = cfacepos[0]+4;
++cfacepos;
}
}
for (j = 0; j < nx * ny; ++j) {
G->cell_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; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
*fnodes++ = i+Nx*j;
*fnodes++ = i+Nx*(j+1);
fnodepos[1] = fnodepos[0] + 2;
++fnodepos;
if (i==0) {
*fcells++ = -1;
*fcells++ = i+nx*j;
}
else if (i == nx) {
*fcells++ = i-1+nx*j;
*fcells++ = -1;
}
else {
*fcells++ = i-1 + nx*j;
*fcells++ = i + nx*j;
}
}
}
/* Faces with y-normal */
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
*fnodes++ = i+1 + Nx*j;
*fnodes++ = i+ Nx*j;
fnodepos[1] = fnodepos[0] + 2;
++fnodepos;
if (j==0) {
*fcells++ = -1;
*fcells++ = i+nx*j;
}
else if (j == ny) {
*fcells++ = i+nx*(j-1);
*fcells++ = -1;
}
else {
*fcells++ = i+nx*(j-1);
*fcells++ = i+nx*j;
}
}
}
}
/* --------------------------------------------------------------------- */
static void
fill_cart_geometry_2d(struct UnstructuredGrid *G,
const double *x,
const double *y)
{
int i,j;
int nx, ny;
double dx, dy;
double *coord, *ccentroids, *cvolumes;
double *fnormals, *fcentroids, *fareas;
nx = G->cartdims[0];
ny = G->cartdims[1];
ccentroids = G->cell_centroids;
cvolumes = G->cell_volumes;
for (j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
*ccentroids++ = (x[i] + x[i + 1]) / 2.0;
*ccentroids++ = (y[j] + y[j + 1]) / 2.0;
dx = x[i + 1] - x[i];
dy = y[j + 1] - y[j];
*cvolumes++ = dx * dy;
}
}
fnormals = G->face_normals;
fcentroids = G->face_centroids;
fareas = G->face_areas;
/* Faces with x-normal */
for (j=0; j<ny; ++j) {
for (i=0; i<nx+1; ++i) {
dy = y[j + 1] - y[j];
*fnormals++ = dy;
*fnormals++ = 0;
*fcentroids++ = x[i];
*fcentroids++ = (y[j] + y[j + 1]) / 2.0;
*fareas++ = dy;
}
}
/* Faces with y-normal */
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx; ++i) {
dx = x[i + 1] - x[i];
*fnormals++ = 0;
*fnormals++ = dx;
*fcentroids++ = (x[i] + x[i + 1]) / 2.0;
*fcentroids++ = y[j];
*fareas++ = dx;
}
}
coord = G->node_coordinates;
for (j=0; j<ny+1; ++j) {
for (i=0; i<nx+1; ++i) {
*coord++ = x[i];
*coord++ = y[j];
}
}
}

172
opm/core/grid/cart_grid.h
View File

@ -1,172 +0,0 @@
/*===========================================================================
//
// File: cart_grid.h
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#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 <CODE>nx + 1</CODE>.
* @param[in] y Position along @c y axis of each grid line with constant @c y
* coordinate. Array of size <CODE>ny + 1</CODE>.
*
* @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 <CODE>nx + 1</CODE>.
* @param[in] y Position along @c y axis of each grid line with constant @c y
* coordinate. Array of size <CODE>ny + 1</CODE>.
* @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
* <CODE>nz + 1</CODE>.
*
* @param[in] depthz
* Top-layer topography specification. If @c NULL, interpreted as
* horizontal top-layer at <CODE>z=0</CODE>. Otherwise, must be
* an array of size <CODE>(nx + 1) * (ny + 1)</CODE>, 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 */

235
opm/core/grid/cornerpoint_grid.c
View File

@ -1,235 +0,0 @@
/*===========================================================================
//
// File: cgridinterface.c
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
*/
#include "config.h"
#include <assert.h>
#include <stdlib.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/grid/cpgpreprocess/geometry.h>
#include <opm/core/grid/cpgpreprocess/preprocess.h>
#include <opm/core/grid.h>
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;
}

103
opm/core/grid/cornerpoint_grid.h
View File

@ -1,103 +0,0 @@
/*===========================================================================
//
// File: preprocess.h
//
// Created: Fri Jun 19 08:43:04 2009
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#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 <opm/core/grid.h>
#include <opm/core/grid/cpgpreprocess/preprocess.h>
#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), <CODE>g->dimensions</CODE> scalars per face
* stored sequentially in <CODE>g->face_centroids</CODE>.
* -# Areas, one scalar per face stored sequentially in
* <CODE>g->face_areas</CODE>.
* -# Normals, <CODE>g->dimensions</CODE> scalars per face stored
* sequentially in <CODE>g->face_normals</CODE>. The Euclidian norm of
* each normal is equal to the corresponding face's area.
*
* - Quantities pertaining to cells (volumes)
* -# Barycenters (centroids), <CODE>g->dimensions</CODE> scalars per cell
* stored sequentially in <CODE>g->cell_centroids</CODE>.
* -# Volumes, one scalar per cell stored sequentially in
* <CODE>g->cell_volumes</CODE>.
*
* 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 */

369
opm/core/grid/cpgpreprocess/facetopology.c
View File

@ -1,369 +0,0 @@
/*===========================================================================
//
// File: facetopology.c
//
// Created: Fri Jun 19 08:46:53 2009
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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<n; ++j) { itop[j]=-1; }
/* Set j to appropriate start position for next i */
j = MIN(k1, k2);
}
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

48
opm/core/grid/cpgpreprocess/facetopology.h
View File

@ -1,48 +0,0 @@
/*===========================================================================
//
// File: facetopology.h
//
// Created: Fri Jun 19 08:47:10 2009
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#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: */

460
opm/core/grid/cpgpreprocess/geometry.c
View File

@ -1,460 +0,0 @@
/*
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
*/
#include "config.h"
#include <math.h>
#include <stdio.h>
#include "geometry.h"
#include <assert.h>
/* ------------------------------------------------------------------ */
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; f<nfaces; ++f)
{
for(i=0; i<ndims; ++i) x[i] = 0.0;
for(i=0; i<ndims; ++i) n[i] = 0.0;
for(i=0; i<ndims; ++i) cface[i] = 0.0;
/* average node */
for(k=nodepos[f]; k<nodepos[f+1]; ++k)
{
node = facenodes[k];
for (i=0; i<ndims; ++i) x[i] += coords[3*node+i];
}
num_face_nodes = nodepos[f+1] - nodepos[f];
for(i=0; i<ndims; ++i) x[i] /= num_face_nodes;
/* compute first vector u (to the last node in the face) */
node = facenodes[nodepos[f+1]-1];
for(i=0; i<ndims; ++i) u[i] = coords[3*node+i] - x[i];
area=0.0;
/* Compute triangular contrib. to face normal and face centroid*/
for(k=nodepos[f]; k<nodepos[f+1]; ++k)
{
node = facenodes[k];
for (i=0; i<ndims; ++i) v[i] = coords[3*node+i] - x[i];
cross(u,v,w);
a = 0.5*norm(w);
area += a;
/* if(!(a>0))
{
fprintf(stderr, "Internal error in compute_face_geometry.");
}
*/
/* face normal */
for (i=0; i<ndims; ++i) n[i] += w[i];
/* face centroid */
for (i=0; i<ndims; ++i)
cface[i] += a*(x[i]+twothirds*0.5*(u[i]+v[i]));
/* Store v in u for next iteration */
for (i=0; i<ndims; ++i) u[i] = v[i];
}
/* Store face normal and face centroid */
for (i=0; i<ndims; ++i)
{
/* normal is scaled with face area */
fnormals [3*f+i] = 0.5*n[i];
fcentroids[3*f+i] = cface[i]/area;
}
fareas[f] = area;
}
}
/* ------------------------------------------------------------------ */
static void
compute_edge_geometry_2d(
/* in */ double *node_coords,
/* in */ int num_edges,
/* in */ int *edge_node_pos,
/* in */ int *edge_nodes,
/* out */ double *edge_normals,
/* out */ double *edge_midpoints,
/* out */ double *edge_lengths)
{
const int num_dims = 2;
/* offsets to each of the nodes in a compacted edge */
const int a_ofs = 0;
const int b_ofs = 1;
/* offsets to each dimension is a compacted point */
const int x_ofs = 0;
const int y_ofs = 1;
int edge; /* edge index */
int a_nod, b_nod; /* node indices */
double a_x, a_y, b_x, b_y; /* node coordinates */
double v_x, v_y; /* vector elements */
/* decompose each edge into a tuple (a,b) between two points and
* compute properties for that face. hopefully the host has enough
* cache pages to keep both input and output at the same time, and
* registers for all the local variables */
for (edge = 0; edge < num_edges; ++edge)
{
/* an edge in 2D can only have starting and ending point
* check that there are exactly two nodes till the next edge */
assert (edge_node_pos[edge + 1] - edge_node_pos[edge] == num_dims);
/* get the first and last point on the edge */
a_nod = edge_nodes[edge_node_pos[edge] + a_ofs];
b_nod = edge_nodes[edge_node_pos[edge] + b_ofs];
/* extract individual coordinates for the points */
a_x = node_coords[a_nod * num_dims + x_ofs];
a_y = node_coords[a_nod * num_dims + y_ofs];
b_x = node_coords[b_nod * num_dims + x_ofs];
b_y = node_coords[b_nod * num_dims + y_ofs];
/* compute edge center -- average of node coordinates */
edge_midpoints[edge * num_dims + x_ofs] = (a_x + b_x) * 0.5;
edge_midpoints[edge * num_dims + y_ofs] = (a_y + b_y) * 0.5;
/* vector from first to last point */
v_x = b_x - a_x;
v_y = b_y - a_y;
/* two-dimensional (unary) cross product analog that makes the
* "triple" (dot-cross) product zero, i.e. it's a normal; the
* direction of this vector is such that it will be pointing
* inwards when enumerating nodes clock-wise */
edge_normals[edge * num_dims + x_ofs] = +v_y;
edge_normals[edge * num_dims + y_ofs] = -v_x;
/* Euclidian norm in two dimensions is magnitude of edge */
edge_lengths[edge] = sqrt(v_x*v_x + v_y*v_y);
}
}
/* ------------------------------------------------------------------ */
void
compute_face_geometry(int ndims, double *coords, int nfaces,
int *nodepos, int *facenodes, double *fnormals,
double *fcentroids, double *fareas)
/* ------------------------------------------------------------------ */
{
if (ndims == 3)
{
compute_face_geometry_3d(coords, nfaces, nodepos, facenodes,
fnormals, fcentroids, fareas);
}
else if (ndims == 2)
{
/* two-dimensional interfaces are called 'edges' */
compute_edge_geometry_2d(coords, nfaces, nodepos, facenodes,
fnormals, fcentroids, fareas);
}
else
{
assert(0);
}
}
/* ------------------------------------------------------------------ */
static void
compute_cell_geometry_3d(double *coords,
int *nodepos, int *facenodes, int *neighbors,
double *fnormals,
double *fcentroids,
int ncells, int *facepos, int *cellfaces,
double *ccentroids, double *cvolumes)
/* ------------------------------------------------------------------ */
{
const int ndims = 3;
int i,k, f,c;
int face,node;
double x[3];
double u[3];
double v[3];
double w[3];
double xcell[3];
double ccell[3];
double cface[3] = {0};
int num_faces;
double volume;
double tet_volume, subnormal_sign;
double twothirds = 0.666666666666666666666666666667;
#pragma omp parallel for default(none) \
private(i,k,f,c,face,node,x,u,v,w,xcell \
,ccell ,cface,num_faces,volume, tet_volume, subnormal_sign) \
shared(coords,nodepos,facenodes,neighbors,twothirds, \
fnormals,fcentroids,facepos,cellfaces,ccentroids,cvolumes) \
firstprivate(ncells)
for (c=0; c<ncells; ++c)
{
for(i=0; i<ndims; ++i) xcell[i] = 0.0;
for(i=0; i<ndims; ++i) ccell[i] = 0.0;
/*
* Approximate cell center as average of face centroids
*/
for(f=facepos[c]; f<facepos[c+1]; ++f)
{
face = cellfaces[f];
for (i=0; i<ndims; ++i) xcell[i] += fcentroids[3*face+i];
}
num_faces = facepos[c+1] - facepos[c];
for(i=0; i<ndims; ++i) xcell[i] /= num_faces;
/*
* For all faces, add tetrahedron's volume and centroid to
* 'cvolume' and 'ccentroid'.
*/
volume=0.0;
for(f=facepos[c]; f<facepos[c+1]; ++f)
{
int num_face_nodes;
for(i=0; i<ndims; ++i) x[i] = 0.0;
for(i=0; i<ndims; ++i) cface[i] = 0.0;
face = cellfaces[f];
/* average face node x */
for(k=nodepos[face]; k<nodepos[face+1]; ++k)
{
node = facenodes[k];
for (i=0; i<ndims; ++i) x[i] += coords[3*node+i];
}
num_face_nodes = nodepos[face+1] - nodepos[face];
for(i=0; i<ndims; ++i) x[i] /= num_face_nodes;
/* compute first vector u (to the last node in the face) */
node = facenodes[nodepos[face+1]-1];
for(i=0; i<ndims; ++i) u[i] = coords[3*node+i] - x[i];
/* Compute triangular contributions to face normal and face centroid */
for(k=nodepos[face]; k<nodepos[face+1]; ++k)
{
node = facenodes[k];
for (i=0; i<ndims; ++i) v[i] = coords[3*node+i] - x[i];
cross(u,v,w);
tet_volume = 0.0;
for(i=0; i<ndims; ++i){
tet_volume += w[i]*(x[i]-xcell[i]);
}
tet_volume *= 0.5 / 3;
subnormal_sign=0.0;
for(i=0; i<ndims; ++i){
subnormal_sign += w[i]*fnormals[3*face+i];
}
if(subnormal_sign < 0.0){
tet_volume = -tet_volume;
}
if(!(neighbors[2*face+0]==c)){
tet_volume = -tet_volume;
}
volume += tet_volume;
/* face centroid of triangle */
for (i=0; i<ndims; ++i) cface[i] = (x[i]+(twothirds)*0.5*(u[i]+v[i]));
/* Cell centroid */
for (i=0; i<ndims; ++i) ccell[i] += tet_volume * 3/4.0*(cface[i] - xcell[i]);
/* Store v in u for next iteration */
for (i=0; i<ndims; ++i) u[i] = v[i];
}
}
for (i=0; i<ndims; ++i) ccentroids[3*c+i] = xcell[i] + ccell[i]/volume;
cvolumes[c] = volume;
}
}
/* ------------------------------------------------------------------ */
static void
compute_cell_geometry_2d(
/* in */ double *node_coords,
/* in */ int *edge_node_pos,
/* in */ int *edge_nodes,
/* in */ double *edge_midpoints,
/* in */ int num_cells,
/* in */ int *cell_edge_pos,
/* in */ int *cell_edges,
/* out */ double *cell_centers,
/* out */ double *cell_areas)
{
const int num_dims = 2;
/* offsets to each of the nodes in a compacted edge */
const int a_ofs = 0;
const int b_ofs = 1;
/* offsets to each dimension is a compacted point */
const int x_ofs = 0;
const int y_ofs = 1;
int cell; /* cell index */
int num_nodes; /* number of vertices in current cell */
int edge_ndx; /* relative edge index within cell */
int edge; /* absolute cell index */
double center_x; /* x-coordinate for cell barycenter */
double center_y; /* y-coordinate for cell barycenter */
double area; /* (accumulated) cell area */
int a_nod, b_nod; /* node indices for edge start and end points */
double a_x, a_y,
b_x, b_y; /* vectors from center to edge points */
for (cell = 0; cell < num_cells; ++cell)
{
/* since the cell is a closed polygon, each point serves as the starting
* point of one edge and the ending point of another; thus there is as
* many vertices as there are edges */
num_nodes = cell_edge_pos[cell + 1] - cell_edge_pos[cell];
/* to enumerate all vertices of a cell, we would have to expand the
* edges and then remove duplicates. however, the centroid of each
* edge contains half of the two vertices that are incident on it. if
* we instead sum all the face centroids, we get the sum of all the
* vertices */
center_x = 0.;
center_y = 0.;
for (edge_ndx = cell_edge_pos[cell];
edge_ndx < cell_edge_pos[cell + 1]; ++edge_ndx)
{
edge = cell_edges[edge_ndx];
center_x += edge_midpoints[edge * num_dims + x_ofs];
center_y += edge_midpoints[edge * num_dims + y_ofs];
}
center_x /= (double) num_nodes;
center_y /= (double) num_nodes;
cell_centers[cell * num_dims + x_ofs] = center_x;
cell_centers[cell * num_dims + y_ofs] = center_y;
/* triangulate the polygon by introducing the cell center and then new
* internal edges from this center to the vertices. the total area of
* the cell is the sum of area of these sub-triangles */
area = 0.;
for (edge_ndx = cell_edge_pos[cell];
edge_ndx < cell_edge_pos[cell + 1]; ++edge_ndx)
{
/* indirect lookup of edge index (from array that contains all the
* edge indices for a certain cell) */
edge = cell_edges[edge_ndx];
/* get the first and last point on the edge */
a_nod = edge_nodes[edge_node_pos[edge] + a_ofs];
b_nod = edge_nodes[edge_node_pos[edge] + b_ofs];
/* vector from center to each of the nodes */
a_x = node_coords[a_nod * num_dims + x_ofs] - center_x;
a_y = node_coords[a_nod * num_dims + y_ofs] - center_y;
b_x = node_coords[b_nod * num_dims + x_ofs] - center_x;
b_y = node_coords[b_nod * num_dims + y_ofs] - center_y;
/* two-dimensional (binary) cross product analog that has length
* equal to the parallelogram spanned by the two vectors (but which
* is a scalar). the sign tells us the orientation between the nodes
* a and b, but we are not interested in that, just the area */
area += fabs(a_x * b_y - a_y * b_x);
}
/* we summed parallelograms which are twice the size of the triangles
* that make up the cell; divide out the half for all terms here */
area *= 0.5;
cell_areas[cell] = area;
}
}
/* ------------------------------------------------------------------ */
void
compute_cell_geometry(int ndims, double *coords,
int *nodepos, int *facenodes, int *neighbors,
double *fnormals,
double *fcentroids,
int ncells, int *facepos, int *cellfaces,
double *ccentroids, double *cvolumes)
/* ------------------------------------------------------------------ */
{
if (ndims == 3)
{
compute_cell_geometry_3d(coords, nodepos, facenodes,
neighbors, fnormals, fcentroids, ncells,
facepos, cellfaces, ccentroids, cvolumes);
}
else if (ndims == 2)
{
compute_cell_geometry_2d(coords, nodepos, facenodes, fcentroids,
ncells, facepos, cellfaces, ccentroids,
cvolumes);
}
else
{
assert(0);
}
}

19
opm/core/grid/cpgpreprocess/geometry.h
View File

@ -1,19 +0,0 @@
/*
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
*/
#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 */

956
opm/core/grid/cpgpreprocess/preprocess.c
View File

@ -1,956 +0,0 @@
/*===========================================================================
//
// File: preprocess.c
//
// Created: Fri Jun 19 08:42:39 2009
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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 <field> with k index running faster than i running
faster than j, and Cartesian dimensions <dims>, 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, ..., #<active cells>. 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; j<ny; ++j) {
for (i=0; i<nx; ++i) {
if (! checkmemory(nz, out, intersections)) {
fprintf(stderr,
"Could not allocate enough space in "
"process_horizontal_faces()\n");
exit(1);
}
f = out->face_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; k<nz*2+1; ++k){
/* Skip if space between face k and face k+1 is collapsed. */
/* Note that inactive cells (with ACTNUM==0) have all been */
/* collapsed in finduniquepoints. */
if (c[0][k] == c[0][k+1] && c[1][k] == c[1][k+1] &&
c[2][k] == c[2][k+1] && c[3][k] == c[3][k+1]){
/* If the pinch is a cell: */
if (k%2){
idx = linearindex(out->dimensions, 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; k<n; ++k){
approximate_intersection_pt(itsct, out->node_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: */

154
opm/core/grid/cpgpreprocess/preprocess.h
View File

@ -1,154 +0,0 @@
/*===========================================================================
//
// File: preprocess.h
//
// Created: Fri Jun 19 08:43:04 2009
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#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: */

383
opm/core/grid/cpgpreprocess/uniquepoints.c
View File

@ -1,383 +0,0 @@
/*===========================================================================
//
// File: uniquepoints.c
//
// Created: Fri Jun 19 08:46:05 2009
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <assert.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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<n; ++i){
for (j=0; j<m; ++j){
if (a[j][i/2]) *ptr++ = z[j][i];
/* else fprintf(stderr, "skipping point in inactive cell\n"); */
}
}
qsort(list, ptr-list, sizeof(double), compare);
return ptr-list;
}
/*-----------------------------------------------------------------
Remove points less than <tolerance> apart in <list> 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<n; ++i){
if (val + tolerance < list [i]){
val = list[i];
list[pos++] = val;
}
}
/*
Preserve outer z-boundary.
This operation is a no-op in the case
list[pos-2] + tolerance < list[n-1].
If, however, the second to last point is less than <tolerance>
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<n; ++i){
/* Skip inactive cells */
if (!a[i/2]) {
p[0] = p[-1]; /* Inactive cells are collapsed leaving
* void space.*/
++p;
continue;
}
/* Find next k such that zlist[k] < z[i] < zlist[k+1] */
while ((k < end) && (zlist[k] + tolerance < z[i])){
k++;
}
/* assert (k < len && z[i] - zlist[k] <= tolerance) */
if ((k == end) || ( zlist[k] + tolerance < z[i])){
fprintf(stderr, "Cannot associate zcorn values with given list\n");
fprintf(stderr, "of z-coordinates to given tolerance\n");
return 0;
}
*p++ = k;
}
*p++ = INT_MAX;/* Padding to ease processing of faults */
return 1;
}
/* ---------------------------------------------------------------------- */
static void
vector_positions(const int dims[3] ,
const int i ,
const int j ,
size_t start[4])
/* ---------------------------------------------------------------------- */
{
size_t im, ip, jm, jp;
im = MAX(1, i ) - 1;
jm = MAX(1, j ) - 1;
ip = MIN(dims[0], i+1) - 1;
jp = MIN(dims[1], j+1) - 1;
start[ 0 ] = dims[2] * (im + dims[0]*jm);
start[ 1 ] = dims[2] * (im + dims[0]*jp);
start[ 2 ] = dims[2] * (ip + dims[0]*jm);
start[ 3 ] = dims[2] * (ip + dims[0]*jp);
}
/*-----------------------------------------------------------------
Given a vector <field> with k index running faster than i running
faster than j, and Cartesian dimensions <dims>, 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 <field> with k index running faster than i running
faster than j, and Cartesian dimensions <dims>, 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; k<len; ++k){
pt[2] = zout[k];
interpolate_pillar(coord, pt);
pt += 3;
}
/* Increment pointer to sparse table of unique zcorn
* values */
zout = zout + len;
zptr[pos++] = zout - zlist;
coord += 6;
}
}
out->number_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: */

47
opm/core/grid/cpgpreprocess/uniquepoints.h
View File

@ -1,47 +0,0 @@
/*===========================================================================
//
// File: uniquepoints.h
//
// Created: Fri Jun 19 08:46:30 2009
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#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: */

562
opm/core/grid/grid.c
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/grid.h>
#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
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;
}

76
opm/core/grid/grid_equal.cpp
View File

@ -1,76 +0,0 @@
#include <string.h> // C string.h to get memcmp()
#include <opm/common/util/numeric/cmp.hpp>
#include <opm/core/grid.h>
/*
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<double>( grid1->node_coordinates , grid2->node_coordinates , static_cast<size_t>(grid1->dimensions * grid1->number_of_nodes)))
return false;
if (!Opm::cmp::array_equal<double>( grid1->face_centroids , grid2->face_centroids , static_cast<size_t>(grid1->dimensions * grid1->number_of_faces)))
return false;
if (!Opm::cmp::array_equal<double>( grid1->face_areas , grid2->face_areas , static_cast<size_t>(grid1->number_of_faces)))
return false;
if (!Opm::cmp::array_equal<double>( grid1->face_normals , grid2->face_normals , static_cast<size_t>(grid1->dimensions * grid1->number_of_faces)))
return false;
if (!Opm::cmp::array_equal<double>( grid1->cell_centroids , grid2->cell_centroids , static_cast<size_t>(grid1->dimensions * grid1->number_of_cells)))
return false;
if (!Opm::cmp::array_equal<double>( grid1->cell_volumes , grid2->cell_volumes , static_cast<size_t>(grid1->number_of_cells)))
return false;
}
return true;
} else
return false;
}

244
opm/core/utility/SparseTable.hpp
View File

@ -1,244 +0,0 @@
//===========================================================================
//
// File: SparseTable.hpp
//
// Created: Fri Apr 24 09:50:27 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SPARSETABLE_HEADER
#define OPM_SPARSETABLE_HEADER
#include <vector>
#include <numeric>
#include <algorithm>
#include <boost/range/iterator_range.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <ostream>
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 <typename T>
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 <typename DataIter, typename IntegerIter>
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 <typename DataIter, typename IntegerIter>
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 <typename IntegerIter>
void allocate(IntegerIter rowsize_beg, IntegerIter rowsize_end)
{
typedef typename std::vector<T>::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 <typename DataIter>
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<T>
void swap(SparseTable<T>& 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<const T*> row_type;
typedef boost::iterator_range<T*> 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<class charT, class traits>
void print(std::basic_ostream<charT, traits>& os) const
{
os << "Number of rows: " << size() << '\n';
os << "Row starts = [";
std::copy(row_start_.begin(), row_start_.end(),
std::ostream_iterator<int>(os, " "));
os << "\b]\n";
os << "Data values = [";
std::copy(data_.begin(), data_.end(),
std::ostream_iterator<T>(os, " "));
os << "\b]\n";
}
const T data(int i)const {
return data_[i];
}
private:
std::vector<T> data_;
// Like in the compressed row sparse matrix format,
// row_start_.size() is equal to the number of rows + 1.
std::vector<int> row_start_;
template <class IntegerIter>
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

106
opm/core/utility/StopWatch.cpp
View File

@ -1,106 +0,0 @@
//===========================================================================
//
// File: StopWatch.cpp
//
// Created: Thu Jul 2 23:04:51 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include <boost/date_time/posix_time/posix_time.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/common/ErrorMacros.hpp>
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

81
opm/core/utility/StopWatch.hpp
View File

@ -1,81 +0,0 @@
//===========================================================================
//
// File: StopWatch.hpp
//
// Created: Thu Jul 2 23:04:17 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_STOPWATCH_HEADER
#define OPM_STOPWATCH_HEADER
#include <boost/date_time/posix_time/posix_time.hpp>
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

132
tests/test_sparsetable.cpp
View File

@ -1,132 +0,0 @@
//===========================================================================
//
// File: test_sparsetable.cpp
//
// Created: Thu May 28 10:01:46 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <bard.skaflestad@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#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 <boost/test/unit_test.hpp>
#include <opm/core/utility/SparseTable.hpp>
using namespace Opm;
BOOST_AUTO_TEST_CASE(construction_and_queries)
{
const SparseTable<int> 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
// <empty row>
// 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<int> 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<int> st2_again(elem, elem + num_elem, rowsizes, rowsizes + num_rows);
BOOST_CHECK(st2 == st2_again);
SparseTable<int> 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<int> 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<int> 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<int> st_empty;
BOOST_CHECK(st2_append2 == st_empty);
SparseTable<int> 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<int>::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<int> st3(elem, elem + num_elem - 1, rowsizes, rowsizes + num_rows), std::exception);
// A few elements too many.
BOOST_CHECK_THROW(const SparseTable<int> st4(elem, elem + num_elem, rowsizes, rowsizes + num_rows - 1), std::exception);
// Need at least one row.
BOOST_CHECK_THROW(const SparseTable<int> 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<int> st6(elem, elem + num_elem, err_rs, err_rs + num_rows), std::exception);
#endif
}