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

Merged master into ert branch

Browse Source
This commit is contained in:
Joakim Hove 2012-06-27 08:22:26 +02:00
commit d06028cb49
20 changed files with 2144 additions and 1819 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

222
Makefile.am
View File

@ -42,6 +42,7 @@ $(ERT_LDFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
#
# Please try to keep the list sorted.
<<<<<<< HEAD
libopmcore_la_SOURCES = \
opm/core/GridManager.cpp \
opm/core/eclipse/EclipseGridInspector.cpp \
@ -265,6 +266,227 @@ opm/core/wells/InjectionSpecification.hpp \
opm/core/wells/ProductionSpecification.hpp \
opm/core/wells/WellCollection.hpp \
opm/core/wells/WellsGroup.hpp \
=======
libopmcore_la_SOURCES = \
opm/core/GridManager.cpp \
opm/core/eclipse/EclipseGridInspector.cpp \
opm/core/eclipse/EclipseGridParser.cpp \
opm/core/fluid/BlackoilPropertiesBasic.cpp \
opm/core/fluid/BlackoilPropertiesFromDeck.cpp \
opm/core/fluid/IncompPropertiesBasic.cpp \
opm/core/fluid/IncompPropertiesFromDeck.cpp \
opm/core/fluid/PvtPropertiesBasic.cpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp \
opm/core/fluid/RockBasic.cpp \
opm/core/fluid/RockCompressibility.cpp \
opm/core/fluid/RockFromDeck.cpp \
opm/core/fluid/SaturationPropsBasic.cpp \
opm/core/fluid/SaturationPropsFromDeck.cpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \
opm/core/fluid/blackoil/SinglePvtDead.cpp \
opm/core/fluid/blackoil/SinglePvtInterface.cpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.cpp \
opm/core/grid.c \
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/linalg/LinearSolverFactory.cpp \
opm/core/linalg/LinearSolverInterface.cpp \
opm/core/linalg/sparse_sys.c \
opm/core/newwells.c \
opm/core/pressure/CompressibleTpfa.cpp \
opm/core/pressure/FlowBCManager.cpp \
opm/core/pressure/IncompTpfa.cpp \
opm/core/pressure/cfsh.c \
opm/core/pressure/flow_bc.c \
opm/core/pressure/fsh.c \
opm/core/pressure/fsh_common_impl.c \
opm/core/pressure/ifsh.c \
opm/core/pressure/mimetic/hybsys.c \
opm/core/pressure/mimetic/hybsys_global.c \
opm/core/pressure/mimetic/mimetic.c \
opm/core/pressure/msmfem/coarse_conn.c \
opm/core/pressure/msmfem/coarse_sys.c \
opm/core/pressure/msmfem/dfs.c \
opm/core/pressure/msmfem/hash_set.c \
opm/core/pressure/msmfem/ifsh_ms.c \
opm/core/pressure/msmfem/partition.c \
opm/core/pressure/tpfa/cfs_tpfa.c \
opm/core/pressure/tpfa/cfs_tpfa_residual.c \
opm/core/pressure/tpfa/compr_bc.c \
opm/core/pressure/tpfa/compr_quant.c \
opm/core/pressure/tpfa/compr_quant_general.c \
opm/core/pressure/tpfa/compr_source.c \
opm/core/pressure/tpfa/ifs_tpfa.c \
opm/core/pressure/tpfa/trans_tpfa.c \
opm/core/pressure/well.c \
opm/core/simulator/SimulatorReport.cpp \
opm/core/simulator/SimulatorTimer.cpp \
opm/core/simulator/SimulatorTwophase.cpp \
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \
opm/core/transport/reorder/TransportModelInterface.cpp \
opm/core/transport/reorder/TransportModelTwophase.cpp \
opm/core/transport/reorder/nlsolvers.c \
opm/core/transport/reorder/reordersequence.cpp \
opm/core/transport/reorder/tarjan.c \
opm/core/transport/spu_explicit.c \
opm/core/transport/spu_implicit.c \
opm/core/transport/transport_source.c \
opm/core/utility/MonotCubicInterpolator.cpp \
opm/core/utility/StopWatch.cpp \
opm/core/utility/miscUtilities.cpp \
opm/core/utility/miscUtilitiesBlackoil.cpp \
opm/core/utility/parameters/Parameter.cpp \
opm/core/utility/parameters/ParameterGroup.cpp \
opm/core/utility/parameters/ParameterTools.cpp \
opm/core/utility/parameters/ParameterXML.cpp \
opm/core/utility/parameters/tinyxml/tinystr.cpp \
opm/core/utility/parameters/tinyxml/tinyxml.cpp \
opm/core/utility/parameters/tinyxml/tinyxmlerror.cpp \
opm/core/utility/parameters/tinyxml/tinyxmlparser.cpp \
opm/core/utility/writeVtkData.cpp \
opm/core/vag_format/vag.cpp \
opm/core/wells/InjectionSpecification.cpp \
opm/core/wells/ProductionSpecification.cpp \
opm/core/wells/WellCollection.cpp \
opm/core/wells/WellsGroup.cpp \
opm/core/wells/WellsManager.cpp
nobase_include_HEADERS = \
opm/core/GridAdapter.hpp \
opm/core/GridManager.hpp \
opm/core/eclipse/CornerpointChopper.hpp \
opm/core/eclipse/EclipseGridInspector.hpp \
opm/core/eclipse/EclipseGridParser.hpp \
opm/core/eclipse/EclipseGridParserHelpers.hpp \
opm/core/eclipse/EclipseUnits.hpp \
opm/core/eclipse/SpecialEclipseFields.hpp \
opm/core/fluid/BlackoilPropertiesBasic.hpp \
opm/core/fluid/BlackoilPropertiesFromDeck.hpp \
opm/core/fluid/BlackoilPropertiesInterface.hpp \
opm/core/fluid/IncompPropertiesBasic.hpp \
opm/core/fluid/IncompPropertiesFromDeck.hpp \
opm/core/fluid/IncompPropertiesInterface.hpp \
opm/core/fluid/PvtPropertiesBasic.hpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \
opm/core/fluid/RockBasic.hpp \
opm/core/fluid/RockCompressibility.hpp \
opm/core/fluid/RockFromDeck.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \
opm/core/fluid/SimpleFluid2p.hpp \
opm/core/fluid/blackoil/BlackoilPhases.hpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \
opm/core/fluid/blackoil/SinglePvtDead.hpp \
opm/core/fluid/blackoil/SinglePvtInterface.hpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.hpp \
opm/core/fluid/blackoil/phaseUsageFromDeck.hpp \
opm/core/grid.h \
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/grdecl.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/blas_lapack.h \
opm/core/linalg/sparse_sys.h \
opm/core/newwells.h \
opm/core/pressure/CompressibleTpfa.hpp \
opm/core/pressure/FlowBCManager.hpp \
opm/core/pressure/HybridPressureSolver.hpp \
opm/core/pressure/IncompTpfa.hpp \
opm/core/pressure/TPFACompressiblePressureSolver.hpp \
opm/core/pressure/TPFAPressureSolver.hpp \
opm/core/pressure/flow_bc.h \
opm/core/pressure/fsh.h \
opm/core/pressure/fsh_common_impl.h \
opm/core/pressure/mimetic/hybsys.h \
opm/core/pressure/mimetic/hybsys_global.h \
opm/core/pressure/mimetic/mimetic.h \
opm/core/pressure/msmfem/coarse_conn.h \
opm/core/pressure/msmfem/coarse_sys.h \
opm/core/pressure/msmfem/dfs.h \
opm/core/pressure/msmfem/hash_set.h \
opm/core/pressure/msmfem/ifsh_ms.h \
opm/core/pressure/msmfem/partition.h \
opm/core/pressure/tpfa/cfs_tpfa.h \
opm/core/pressure/tpfa/cfs_tpfa_residual.h \
opm/core/pressure/tpfa/compr_bc.h \
opm/core/pressure/tpfa/compr_quant.h \
opm/core/pressure/tpfa/compr_quant_general.h \
opm/core/pressure/tpfa/compr_source.h \
opm/core/pressure/tpfa/ifs_tpfa.h \
opm/core/pressure/tpfa/trans_tpfa.h \
opm/core/simulator/BlackoilState.hpp \
opm/core/simulator/SimulatorReport.hpp \
opm/core/simulator/SimulatorTimer.hpp \
opm/core/simulator/SimulatorTwophase.hpp \
opm/core/simulator/TwophaseState.hpp \
opm/core/simulator/WellState.hpp \
opm/core/transport/CSRMatrixBlockAssembler.hpp \
opm/core/transport/CSRMatrixUmfpackSolver.hpp \
opm/core/transport/GravityColumnSolver.hpp \
opm/core/transport/GravityColumnSolver_impl.hpp \
opm/core/transport/ImplicitAssembly.hpp \
opm/core/transport/ImplicitTransport.hpp \
opm/core/transport/JacobianSystem.hpp \
opm/core/transport/NormSupport.hpp \
opm/core/transport/SimpleFluid2pWrapper.hpp \
opm/core/transport/SinglePointUpwindTwoPhase.hpp \
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp \
opm/core/transport/reorder/TransportModelInterface.hpp \
opm/core/transport/reorder/TransportModelTwophase.hpp \
opm/core/transport/reorder/nlsolvers.h \
opm/core/transport/reorder/reordersequence.h \
opm/core/transport/reorder/tarjan.h \
opm/core/transport/spu_explicit.h \
opm/core/transport/spu_implicit.h \
opm/core/transport/transport_source.h \
opm/core/utility/Average.hpp \
opm/core/utility/ColumnExtract.hpp \
opm/core/utility/ErrorMacros.hpp \
opm/core/utility/Factory.hpp \
opm/core/utility/MonotCubicInterpolator.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/Units.hpp \
opm/core/utility/buildUniformMonotoneTable.hpp \
opm/core/utility/initState.hpp \
opm/core/utility/initState_impl.hpp \
opm/core/utility/linInt.hpp \
opm/core/utility/linearInterpolation.hpp \
opm/core/utility/miscUtilities.hpp \
opm/core/utility/miscUtilitiesBlackoil.hpp \
opm/core/utility/parameters/Parameter.hpp \
opm/core/utility/parameters/ParameterGroup.hpp \
opm/core/utility/parameters/ParameterGroup_impl.hpp \
opm/core/utility/parameters/ParameterMapItem.hpp \
opm/core/utility/parameters/ParameterRequirement.hpp \
opm/core/utility/parameters/ParameterStrings.hpp \
opm/core/utility/parameters/ParameterTools.hpp \
opm/core/utility/parameters/ParameterXML.hpp \
opm/core/utility/parameters/tinyxml/tinystr.h \
opm/core/utility/parameters/tinyxml/tinyxml.h \
opm/core/utility/writeVtkData.hpp \
opm/core/vag_format/vag.hpp \
opm/core/well.h \
opm/core/wells/InjectionSpecification.hpp \
opm/core/wells/ProductionSpecification.hpp \
opm/core/wells/WellCollection.hpp \
opm/core/wells/WellsGroup.hpp \
>>>>>>> master
opm/core/wells/WellsManager.hpp
# ----------------------------------------------------------------------

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

@ -32,20 +32,19 @@ You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <stdio.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "preprocess.h"
#include "sparsetable.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 MIN(i,j) ((i)<(j) ? (i) : (j))
#define MAX(i,j) ((i)>(j) ? (i) : (j))
#define DEBUG 1
@ -59,17 +58,18 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
/*------------------------------------------------------*/
/* Determine face geometry first, then compute intersections. */
/* Determine face topology first, then compute intersection. */
/* All intersections that occur are present in the final face geometry.*/
static int *computeFaceTopology(int *a1,
int *a2,
int *b1,
int *b2,
int intersect[4],
int *faces)
static int *
computeFaceTopology(const int *a1, const int *a2,
const int *b1, const int *b2,
int intersect[4], int *faces)
{
int mask[8] = {-1};
int k, *f;
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]; }
@ -146,7 +146,6 @@ static int *computeFaceTopology(int *a1,
#if DEBUG>1
/* Check for repeated nodes:*/
{
int i;
fprintf(stderr, "face: ");
for (i=0; i<8; ++i){
@ -158,7 +157,6 @@ static int *computeFaceTopology(int *a1,
}
}
fprintf(stderr, "\n");
}
#endif
@ -171,27 +169,37 @@ static int *computeFaceTopology(int *a1,
/* a) If we assume that the index increase when z increase for
each pillar (but only separately), we can use only the point indices.
/* 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 no intersections occur on the first and last lines.
This is convenient in the identification of (unique) intersections.
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 lineintersection(a1,a2,b1,b2)(((a1>b1)&&(a2<b2))||((a1<b1)&&(a2>b2)))
static int faceintersection(int *a1, int *a2, int *b1, int *b2)
#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]) ||
lineintersection(a1[0], a2[0], b1[0], b2[0]) ||
lineintersection(a1[1], a2[1], b1[1], b2[1]);
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 !((a1[i] ==INT_MIN) && (b1[j] ==INT_MIN)) && \
!((a1[i+1]==INT_MAX) && (b1[j+1]==INT_MAX))
#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],
@ -210,24 +218,34 @@ void findconnections(int n, int *pts[4],
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 *tmp;
int k1 = 0;
int k2 = 0;
int i,j=0;
int intersect[4]= {-1};
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) {
if (a1[i] == a1[i+1] && a2[i] == a2[i+1]) continue;
/* 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])) {
while(j<n-1 && (b1[j] < a1[i+1] || b2[j] < a2[i+1])){
if (b1[j] == b1[j+1] && b2[j] == b2[j+1]){
itop[j+1] = itop[j];
++j;
continue;
@ -241,14 +259,13 @@ void findconnections(int n, int *pts[4],
/* 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 (((a1[i] == b1[j]) && (a1[i+1] == b1[j+1])) &&
((a2[i] == b2[j]) && (a2[i+1] == b2[j+1]))) {
/* Add face to list of faces if not any first or last points are involved. */
if (meaningful_face){
if (MEANINGFUL_FACE(i, j)) {
int cell_a = i%2 ? (i-1)/2 : -1;
int cell_b = j%2 ? (j-1)/2 : -1;
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;
@ -257,21 +274,17 @@ void findconnections(int n, int *pts[4],
/* 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];
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{
;
/*
fprintf(stderr,
"Warning. For some reason we get face connecting void spaces\n");
*/
}
}
}
@ -279,7 +292,8 @@ void findconnections(int n, int *pts[4],
else{
/* Find new intersection */
if (lineintersection(a1[i+1],a2[i+1],b1[j+1],b2[j+1])) {
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 */
@ -300,15 +314,16 @@ void findconnections(int n, int *pts[4],
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){
/* 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 ? (i-1)/2 : -1;
int cell_b = j%2 ? (j-1)/2 : -1;
int cell_a = i%2 != 0 ? (i-1)/2 : -1;
int cell_b = j%2 != 0 ? (j-1)/2 : -1;
@ -322,32 +337,29 @@ void findconnections(int n, int *pts[4],
}
else{
;
/*
fprintf(stderr,
"Warning. For some reason we get face connecting void spaces\n");
*/
}
}
}
}
/* 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;
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] */
/* 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;
for(j=0;j<n; ++j) { itop[j]=-1; }
/* Set j to appropriate start position for next i */
j = min(k1, k2);
j = MIN(k1, k2);
}
}

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

@ -32,8 +32,8 @@ You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_FACETOPOLOGY_HEADER
#define OPENRS_FACETOPOLOGY_HEADER
#ifndef OPM_FACETOPOLOGY_HEADER
#define OPM_FACETOPOLOGY_HEADER
void findconnections(int n, int *pts[4],
@ -41,7 +41,7 @@ void findconnections(int n, int *pts[4],
int *work,
struct processed_grid *out);
#endif /* OPENRS_FACETOPOLOGY_HEADER */
#endif /* OPM_FACETOPOLOGY_HEADER */
/* Local Variables: */
/* c-basic-offset:4 */

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

@ -131,6 +131,12 @@ compute_cell_geometry(int ndims, double *coords,
double ccell[3];
double cface[3] = {0};
const double twothirds = 0.666666666666666666666666666667;
int ndigits;
ndigits = ((int) (log(ncells) / log(10.0))) + 1;
for (c=0; c<ncells; ++c)
{
int num_faces;
@ -192,31 +198,25 @@ compute_cell_geometry(int ndims, double *coords,
cross(u,v,w);
tet_volume = 0;
for(i=0; i<ndims; ++i){
tet_volume += 0.5/3 * w[i]*(x[i]-xcell[i]);
/*tet_volume += 0.5/3 * w[i]*(x[i]-xcell[i]);*/
tet_volume += w[i] * (x[i] - xcell[i]);
}
/*tet_volume = fabs(tet_volume);*/
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){
tet_volume*=-1.0;
if (subnormal_sign < 0.0) {
tet_volume = - tet_volume;
}
if(!(neighbors[2*face+0]==c)){
tet_volume*=-1.0;
if (neighbors[2*face + 0] != c) {
tet_volume = - tet_volume;
}
volume += tet_volume;
if(!(volume>0)){
fprintf(stderr, "Internal error in compute_cell_geometry: negative volume\n");
}
/* face centroid of triangle */
for (i=0; i<ndims; ++i) cface[i] = (x[i]+twothirds*0.5*(u[i]+v[i]));
@ -229,6 +229,13 @@ compute_cell_geometry(int ndims, double *coords,
for (i=0; i<ndims; ++i) u[i] = v[i];
}
}
if (! (volume > 0.0)) {
fprintf(stderr,
"Internal error in mex_compute_geometry(%*d): "
"negative volume\n", ndigits, c);
}
for (i=0; i<ndims; ++i) ccentroids[3*c+i] = xcell[i] + ccell[i]/volume;
cvolumes[c] = volume;
}

10
opm/core/grid/cpgpreprocess/grdecl.h
View File

@ -1,5 +1,5 @@
#ifndef GRDECL_H
#define GRDECL_H
#ifndef GRDECL_H_INCLUDED
#define GRDECL_H_INCLUDED
struct grdecl{
int dims[3];
@ -9,4 +9,8 @@ struct grdecl{
};
#endif
#endif /* GRDECL_H_INCLUDED */
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

62
opm/core/grid/cpgpreprocess/mxgrdecl.c
View File

@ -1,4 +1,4 @@
//===========================================================================
/*=========================================================================
//
// File: mxgrdecl.c
//
@ -10,7 +10,7 @@
//
// $Revision$
//
//===========================================================================
//=======================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
@ -32,36 +32,33 @@ You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <mex.h>
#include <grdecl.h>
#include "grdecl.h"
#include "mxgrdecl.h"
/* Get COORD, ZCORN, ACTNUM and DIMS from mxArray. */
/*-------------------------------------------------------*/
void mx_init_grdecl(struct grdecl *g, const mxArray *s)
{
int i,n;
mxArray *field;
int numel;
size_t numel;
mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL;
if (!mxIsStruct(s)
|| !(cartdims = mxGetField(s, 0, "cartDims"))
|| !(actnum = mxGetField(s, 0, "ACTNUM"))
|| !(coord = mxGetField(s, 0, "COORD"))
|| !(zcorn = mxGetField(s, 0, "ZCORN"))
)
{
char str[]="Input must be a single Matlab struct with fields\n"
"cartDims, ACTNUM, COORD and ZCORN\n";
char str[]="Input must be a single MATLAB struct with fields\n"
"'cartDims', 'COORD' and 'ZCORN'. ACTNUM may be included.\n";
mexErrMsgTxt(str);
}
@ -71,36 +68,49 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s)
mexErrMsgTxt("cartDims field must be 3 numbers");
}
if (mxIsDouble(cartdims)) {
double *tmp = mxGetPr(cartdims);
n = 1;
for (i = 0; i < 3; ++i) {
g->dims[i] = tmp[i];
n *= tmp[i];
g->dims[i] = (int) tmp[i];
}
}
else if (mxIsInt32(cartdims)) {
int *tmp = mxGetData(cartdims);
memcpy(g->dims, tmp, 3 * sizeof *g->dims);
}
n = g->dims[0];
for (i = 1; i < 3; i++) { n *= g->dims[ i ]; }
if ((actnum = mxGetField(s, 0, "ACTNUM")) != NULL) {
numel = mxGetNumberOfElements(actnum);
if (mxGetClassID(actnum) != mxINT32_CLASS ||
numel != g->dims[0]*g->dims[1]*g->dims[2] ){
if ((! mxIsInt32(actnum)) || (numel != (size_t)(n))) {
mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32");
}
g->actnum = mxGetData(actnum);
}
else {
g->actnum = NULL;
}
field = mxGetField(s, 0, "COORD");
numel = mxGetNumberOfElements(coord);
if (mxGetClassID(coord) != mxDOUBLE_CLASS ||
numel != 6*(g->dims[0]+1)*(g->dims[1]+1)){
if ((! mxIsDouble(coord)) ||
numel != (size_t)(6*(g->dims[0]+1)*(g->dims[1]+1))) {
mexErrMsgTxt("COORD field must have 6*(nx+1)*(ny+1) doubles.");
}
g->coord = mxGetPr(coord);
numel = mxGetNumberOfElements(zcorn);
if (mxGetClassID(zcorn) != mxDOUBLE_CLASS ||
numel != 8*g->dims[0]*g->dims[1]*g->dims[2]){
if ((! mxIsDouble(zcorn)) ||
numel != (size_t)(8*g->dims[0]*g->dims[1]*g->dims[2])) {
mexErrMsgTxt("ZCORN field must have 8*nx*ny*nz doubles.");
}
g->zcorn = mxGetPr(zcorn);
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

14
opm/core/grid/cpgpreprocess/mxgrdecl.h
View File

@ -1,4 +1,4 @@
//===========================================================================
/*=========================================================================
//
// File: mxgrdecl.h
//
@ -10,7 +10,7 @@
//
// $Revision$
//
//===========================================================================
//=======================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
@ -32,11 +32,13 @@ You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_MXGRDECL_HEADER
#define OPENRS_MXGRDECL_HEADER
#ifndef OPM_MXGRDECL_HEADER
#define OPM_MXGRDECL_HEADER
void mx_init_grdecl (struct grdecl *g, const mxArray *s);
#endif // OPENRS_MXGRDECL_HEADER
#endif /* OPM_MXGRDECL_HEADER */
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

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

@ -32,35 +32,65 @@ You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <stdio.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))
#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[])
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;
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);
@ -76,23 +106,27 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
/*-----------------------------------------------------------------
Special purpose
Convert from k-index to Cartesian index i+nx*(j+ny*k) for every other
element in neighbors.
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)
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] || j<0 || j >= dims[1]){
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;
neighbors[k] = -1; /* Neighbour is outside domain */
}
}else{
}
else {
for (k = 0; k < len; k += 2) {
if (neighbors[k] != -1) {
int tmp = i + dims[0]*(j + dims[1]*neighbors[k]);
neighbors[k] = tmp;
neighbors[k] = linearindex(dims, i, j, neighbors[k]);
}
}
}
@ -100,52 +134,60 @@ static void compute_cell_index(const int dims[3], int i, int j, int *neighbors,
/*-----------------------------------------------------------------
Ensure there's sufficient memory
*/
static int checkmemeory(int nz, struct processed_grid *out, int **intersections)
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 */
int r = (2*nz+2)*(2*nz+2);
int m = out->m;
int n = out->n;
/* 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*0.5, 2*r);
m += MAX(m / 2, 2 * r);
}
if (out->face_ptr[out->number_of_faces] + 6*r > n) {
n += max(n*0.5, 12*r);
n += MAX(n / 2, 12 * r);
}
if (m != out->m){
void *p1 = realloc(out->face_neighbors, 2*m * sizeof(*out->face_neighbors));
void *p2 = realloc(*intersections, 4*m * sizeof(**intersections));
if (p1 && p2) {
out->face_neighbors = p1;
*intersections = p2;
} else {
return 0;
}
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 (m != out->m || n != out->n) {
void *p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes);
void *p2 = realloc(out->face_ptr, (m+1) * sizeof *out->face_ptr);
void *p3 = realloc(out->face_tag, m * sizeof *out->face_tag);
if (ok && (n != out->n)) {
void *p1;
if (p1 && p2 && p3) {
p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes);
ok = p1 != NULL;
if (ok) {
out->face_nodes = p1;
out->face_ptr = p2;
out->face_tag = p3;
} else {
return 0;
}
out->m = m;
out->n = n;
}
}
return 1;
return ok;
}
/*-----------------------------------------------------------------
@ -157,41 +199,52 @@ static int checkmemeory(int nz, struct processed_grid *out, int **intersections)
direction == 0 : constant-i faces.
direction == 1 : constant-j faces.
*/
static void process_vertical_faces(int direction,
static void
process_vertical_faces(int direction,
int **intersections,
int *plist, int *work,
struct processed_grid *out)
{
int i,j;
int d[3];
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;
int startface, num_intersections, f, len, *ptr;
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){
for (i = 0; i < nx + (1 - direction); ++i) {
if (!checkmemeory(nz, out, intersections)){
fprintf(stderr, "Could not allocat enough space in process_vertical_faces\n");
if (! checkmemory(nz, out, intersections)) {
fprintf(stderr,
"Could not allocate enough space in "
"process_vertical_faces()\n");
exit(1);
}
/* Vectors of point numbers */
d[0] = 2*nx; d[1] = 2*ny; d[2] = 2+2*nz;
igetvectors(d, 2*i+direction, 2*j+1-direction, plist, cornerpts);
igetvectors(d, 2*i + direction, 2*j + (1 - direction),
plist, cornerpts);
if (direction == 1) {
/* 1 3 0 1 */
/* ---> */
/* 0 2 2 3 */
/* rotate clockwise */
int *tmp = cornerpts[1];
tmp = cornerpts[1];
cornerpts[1] = cornerpts[0];
cornerpts[0] = cornerpts[2];
cornerpts[2] = cornerpts[3];
@ -204,13 +257,21 @@ static void process_vertical_faces(int direction,
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);
@ -223,11 +284,6 @@ static void process_vertical_faces(int direction,
}
}
static int linearindex(const int dims[3], int i, int j, int k)
{
return i+dims[0]*(j+dims[1]*k);
}
/*-----------------------------------------------------------------
For each horizontal face (i.e. k constant),
@ -240,7 +296,8 @@ static int linearindex(const int dims[3], int i, int j, int k)
ACTNUM==0)
*/
static void process_horizontal_faces(int **intersections,
static void
process_horizontal_faces(int **intersections,
int *plist,
struct processed_grid *out)
{
@ -252,20 +309,25 @@ static void process_horizontal_faces(int **intersections,
int *cell = out->local_cell_index;
int cellno = 0;
int *f, *n, *c[4], prevcell, thiscell;
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;
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 (!checkmemeory(nz, out, intersections)){
fprintf(stderr, "Could not allocat enough space in process_horizontal_faces\n");
if (! checkmemory(nz, out, intersections)) {
fprintf(stderr,
"Could not allocate enough space in "
"process_horizontal_faces()\n");
exit(1);
}
@ -290,7 +352,7 @@ static void process_horizontal_faces(int **intersections,
/* If the pinch is a cell: */
if (k%2){
int idx = linearindex(out->dimensions, i,j,(k-1)/2);
idx = linearindex(out->dimensions, i,j,(k-1)/2);
cell[idx] = -1;
}
}
@ -348,39 +410,72 @@ static void process_horizontal_faces(int **intersections,
*/
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;
double z0 = c[3*L[0]+2];
double z1 = c[3*L[1]+2];
double z2 = c[3*L[2]+2];
double z3 = c[3*L[3]+2];
/* 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 {
double a = (z2-z0)/(z1-z0 - (z3-z2));
if (isinf(a) || isnan(a)){
a = 0;
}
/* the corresponding z-coordinate is */
z = z0*(1.0 - a) + z1*a;
pt[0] = c[3*L[0]+0]* (1.0-a) + c[3*L[1]+0]* a;
pt[1] = c[3*L[0]+1]* (1.0-a) + c[3*L[1]+1]* a;
pt[2] = c[3*L[0]+2]* (1.0-a) + c[3*L[1]+2]* 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.
*/
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, *itsct;
int k;
double *pt;
/* Make sure the space allocated for nodes match the number of node. */
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;
@ -392,7 +487,6 @@ compute_intersection_coordinates(int *intersections,
/* Append intersections */
pt = out->node_coordinates + 3*np;
itsct = intersections;
for (k=np; k<n; ++k){
approximate_intersection_pt(itsct, out->node_coordinates, pt);
@ -403,6 +497,257 @@ compute_intersection_coordinates(int *intersections,
}
/* ------------------------------------------------------------------ */
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);
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, start;
double dz;
const double *zcorn;
ind2sub(in->dims[0], in->dims[1], in->dims[2], c, &i, &j, &k);
zcorn = in->zcorn;
start = (k * off[4]) + (j * off[2]) + (i * 2);
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
*/
@ -410,105 +755,39 @@ void process_grdecl(const struct grdecl *in,
double tolerance,
struct processed_grid *out)
{
int i,j,k;
int nx = in->dims[0];
int ny = in->dims[1];
int nz = in->dims[2];
int nc = nx*ny*nz;
struct grdecl g;
int *actnum, *iptr, sign, error;
double *zcorn , *dptr;
size_t i;
int sign, error, left_handed;
int cellnum;
int *intersections, *work, *plist, *ptr, *global_cell_index;
int *actnum, *iptr;
int *global_cell_index;
const int BIGNUM = 64;
double *zcorn;
g.dims[0] = nx;
g.dims[1] = ny;
g.dims[2] = nz;
actnum = malloc (nc * sizeof *actnum);
zcorn = malloc (nc * 8 * sizeof *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);
/* Permute actnum */
iptr = actnum;
if (in->actnum != NULL) {
for (j=0; j<ny; ++j){
for (i=0; i<nx; ++i){
for (k=0; k<nz; ++k){
*iptr++ = in->actnum[i+nx*(j+ny*k)];
}
}
}
} else {
for (i=0; i<nx*ny*nz; ++i){
*iptr++ = 1;
}
}
g.actnum = actnum;
/* HACK */
/* Check that ZCORN is strictly nodecreaing along pillars. If */
/* not, check if -ZCORN is strictly nondecreasing. */
for (sign = 1; sign>-2; sign = sign - 2){
error = 0;
/* Ensure that zcorn is strictly nondecreasing in k-direction */
for (j=0; j<2*ny; ++j){
for (i=0; i<2*nx; ++i){
for (k=0; k<2*nz-1; ++k){
double z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))];
double z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))];
int c1 = i/2 + nx*(j/2 + ny*(k/2));
int c2 = i/2 + nx*(j/2 + ny*((k+1)/2));
if (((in->actnum == NULL) || (in->actnum[c1] && in->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");
}
/* Permute zcorn */
dptr = zcorn;
for (j=0; j<2*ny; ++j){
for (i=0; i<2*nx; ++i){
for (k=0; k<2*nz; ++k){
*dptr++ = sign*in->zcorn[i+2*nx*(j+2*ny*k)];
}
}
}
/* internal work arrays */
int *work;
int *plist;
int *intersections;
g.zcorn = zcorn;
g.coord = in->coord;
/* Code below assumes k index runs fastests, ie. that dimensions of
table is permuted to (dims[2], dims[0], dims[1]) */
out->m = BIGNUM/3;
out->n = BIGNUM;
/* -----------------------------------------------------------------*/
/* 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);
@ -516,35 +795,74 @@ void process_grdecl(const struct grdecl *in,
out->face_tag = malloc( out->m * sizeof *out->face_tag);
out->face_ptr[0] = 0;
out->dimensions[0] = g.dims[0];
out->dimensions[1] = g.dims[1];
out->dimensions[2] = g.dims[2];
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(nx*ny*nz * sizeof *out->local_cell_index);
out->local_cell_index = malloc(nc * sizeof *out->local_cell_index);
/* internal */
intersections = malloc(BIGNUM* sizeof(*intersections));
/* internal */
work = malloc(2* (2*nz+2)* sizeof(*work));
for(i=0; i<4*(nz+1); ++i) work[i] = -1;
/* Allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */
plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist);
/* 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);
@ -552,34 +870,65 @@ void process_grdecl(const struct grdecl *in,
free (plist);
free (work);
/* -----------------------------------------------------------------*/
/* (re)allocate space for and compute coordinates of nodes that
* arise from intesecting cells (faults) */
compute_intersection_coordinates(intersections, out);
free (intersections);
/* Convert to local cell numbers in face_neighbors */
ptr=out->face_neighbors;;
for (i=0; i<out->number_of_faces*2; ++i, ++ptr){
if (*ptr != -1){
*ptr = out->local_cell_index[*ptr];
}
}
/* Invert global-to-local map */
global_cell_index = malloc(out->number_of_cells *
sizeof (*global_cell_index));
for (i=0; i<nx*ny*nz; ++i){
/* -----------------------------------------------------------------*/
/* 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[out->local_cell_index[i]] = i;
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;
/* HACK: undo sign reversal in ZCORN prepprocessing */
for (i=2; i<3*out->number_of_nodes; i = i+3)
out->node_coordinates[i]=sign*out->node_coordinates[i];
/* 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);
}
}
/*-------------------------------------------------------*/

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

@ -53,8 +53,8 @@ extern "C" {
/* Output structure holding grid topology */
struct processed_grid{
int m;
int n;
int m; /** Upper bound on "number_of_faces" */
int n; /** Upper bound on "number_of_nodes" */
int dimensions[3]; /* Cartesian dimension */
int number_of_faces;
int *face_nodes; /* Nodes numbers of each face sequentially. */
@ -63,7 +63,7 @@ extern "C" {
enum face_tag *face_tag;
int number_of_nodes;
int number_of_nodes_on_pillars;
int number_of_nodes_on_pillars; /** Total number of unique cell vertices that lie on pillars. */
double *node_coordinates; /* 3 doubles per node, sequentially */
int number_of_cells; /* number of active cells */

594
opm/core/grid/cpgpreprocess/processgrid.c
View File

@ -13,8 +13,8 @@
//===========================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2011, 2012 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
@ -36,232 +36,430 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
/* Mex gateway by Jostein R. Natvig, SINTEF ICT. */
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <mex.h>
#include "preprocess.h"
#include "mxgrdecl.h"
void fill_grid(mxArray **out, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
static mxArray *
allocate_nodes(size_t nnodes)
/* ---------------------------------------------------------------------- */
{
const char *names[] = {"nodes", "faces", "cells", "cartDims", "type"};
mxArray *G = mxCreateStructMatrix(1,1,sizeof names / sizeof names[0],names);
size_t nflds;
const char *fields[] = { "num", "coords" };
int i,j;
double *ptr;
mxArray *nodes, *num, *coords;
nflds = sizeof(fields) / sizeof(fields[0]);
nodes = mxCreateStructMatrix(1, 1, nflds, fields);
num = mxCreateDoubleScalar(nnodes);
coords = mxCreateDoubleMatrix(nnodes, 3, mxREAL);
if ((nodes != NULL) && (num != NULL) && (coords != NULL)) {
mxSetField(nodes, 0, "num" , num);
mxSetField(nodes, 0, "coords", coords);
} else {
if (coords != NULL) { mxDestroyArray(coords); }
if (num != NULL) { mxDestroyArray(num); }
if (nodes != NULL) { mxDestroyArray(nodes); }
nodes = NULL;
}
return nodes;
}
/* nodes */
const char *n2[] = {"num", "coords"};
mxArray *nodes = mxCreateStructMatrix(1,1,sizeof n2 / sizeof n2[0],n2);
mxSetField(nodes, 0, "num", mxCreateDoubleScalar(grid->number_of_nodes));
mxArray *nodecoords = mxCreateDoubleMatrix(grid->number_of_nodes, 3, mxREAL);
ptr = mxGetPr(nodecoords);
for (j=0;j<3;++j)
/* ---------------------------------------------------------------------- */
static void
fill_nodes(mxArray *nodes, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
for(i=0; i<grid->number_of_nodes; ++i)
{
ptr[i+grid->number_of_nodes*j] = grid->node_coordinates[3*i+j];
size_t i, j, nnodes;
double *coords;
nnodes = grid->number_of_nodes;
coords = mxGetPr(mxGetField(nodes, 0, "coords"));
for (j = 0; j < 3; ++j) {
for (i = 0; i < nnodes; ++i) {
*coords++ = grid->node_coordinates[3*i + j];
}
}
mxSetField(nodes, 0, "coords", nodecoords);
}
/* ---------------------------------------------------------------------- */
static mxArray *
allocate_faces(size_t nf, size_t nfacenodes)
/* ---------------------------------------------------------------------- */
{
size_t nflds;
const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" };
mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag;
nflds = sizeof(fields) / sizeof(fields[0]);
faces = mxCreateStructMatrix(1, 1, nflds, fields);
num = mxCreateDoubleScalar (nf);
neighbors = mxCreateNumericMatrix(nf , 2, mxINT32_CLASS, mxREAL);
nodePos = mxCreateNumericMatrix(nf + 1 , 1, mxINT32_CLASS, mxREAL);
nodes = mxCreateNumericMatrix(nfacenodes, 1, mxINT32_CLASS, mxREAL);
tag = mxCreateNumericMatrix(nf , 1, mxINT32_CLASS, mxREAL);
if ((faces != NULL) && (num != NULL) && (neighbors != NULL) &&
(nodePos != NULL) && (nodes != NULL) && (tag != NULL)) {
mxSetField(faces, 0, "num" , num);
mxSetField(faces, 0, "neighbors", neighbors);
mxSetField(faces, 0, "nodePos" , nodePos);
mxSetField(faces, 0, "nodes" , nodes);
mxSetField(faces, 0, "tag" , tag);
} else {
if (tag != NULL) { mxDestroyArray(tag); }
if (nodes != NULL) { mxDestroyArray(nodes); }
if (nodePos != NULL) { mxDestroyArray(nodePos); }
if (neighbors != NULL) { mxDestroyArray(neighbors); }
if (num != NULL) { mxDestroyArray(num); }
if (faces != NULL) { mxDestroyArray(faces); }
faces = NULL;
}
return faces;
}
/* ---------------------------------------------------------------------- */
static void
fill_faces(mxArray *faces, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
size_t i, f, nf, nfn;
int *pi;
nf = grid->number_of_faces;
/* Fill faces.neighbors */
pi = mxGetData(mxGetField(faces, 0, "neighbors"));
for (i = 0; i < 2; i++) {
for (f = 0; f < nf; f++) {
/* Add one for one-based indexing in M */
*pi++ = grid->face_neighbors[2*f + i] + 1;
}
}
/* Fill faces.nodePos */
pi = mxGetData(mxGetField(faces, 0, "nodePos"));
for (i = 0; i <= nf; i++) { pi[i] = grid->face_ptr[i] + 1; }
/* Fill faces.nodes */
pi = mxGetData(mxGetField(faces, 0, "nodes"));
nfn = grid->face_ptr[nf]; /* Total number of face nodes */
for (i = 0; i < nfn; i++) { pi[i] = grid->face_nodes[i] + 1; }
/* Fill faces.tag */
pi = mxGetData(mxGetField(faces, 0, "tag"));
for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[f] + 1; }
}
/* ---------------------------------------------------------------------- */
static size_t
count_halffaces(size_t nf, const int *neighbors)
/* ---------------------------------------------------------------------- */
{
int c1, c2;
size_t nhf, f;
for (f = nhf = 0; f < nf; f++) {
c1 = neighbors[2*f + 0];
c2 = neighbors[2*f + 1];
nhf += c1 >= 0;
nhf += c2 >= 0;
}
return nhf;
}
/* ---------------------------------------------------------------------- */
static mxArray *
allocate_cells(size_t nc, size_t ncf)
/* ---------------------------------------------------------------------- */
{
size_t nflds;
const char *fields[] = { "num", "facePos", "faces", "indexMap" };
mxArray *cells, *num, *facePos, *faces, *indexMap;
nflds = sizeof(fields) / sizeof(fields[0]);
cells = mxCreateStructMatrix(1, 1, nflds, fields);
num = mxCreateDoubleScalar (nc);
facePos = mxCreateNumericMatrix(nc + 1, 1, mxINT32_CLASS, mxREAL);
faces = mxCreateNumericMatrix(ncf , 2, mxINT32_CLASS, mxREAL);
indexMap = mxCreateNumericMatrix(nc , 1, mxINT32_CLASS, mxREAL);
if ((cells != NULL) && (num != NULL) && (facePos != NULL) &&
(faces != NULL) && (indexMap != NULL)) {
mxSetField(cells, 0, "num" , num );
mxSetField(cells, 0, "facePos" , facePos );
mxSetField(cells, 0, "faces" , faces );
mxSetField(cells, 0, "indexMap", indexMap);
} else {
if (indexMap != NULL) { mxDestroyArray(indexMap); }
if (faces != NULL) { mxDestroyArray(faces); }
if (facePos != NULL) { mxDestroyArray(facePos); }
if (num != NULL) { mxDestroyArray(num); }
if (cells != NULL) { mxDestroyArray(cells); }
cells = NULL;
}
return cells;
}
/* ---------------------------------------------------------------------- */
static void
fill_cells(mxArray *cells, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
size_t c, nc, f, nf, i;
int c1, c2, cf_tag, nhf;
int *pi1, *pi2;
nc = grid->number_of_cells;
nf = grid->number_of_faces;
/* Simultaneously fill cells.facePos and cells.faces by transposing the
* neighbours mapping. */
pi1 = mxGetData(mxGetField(cells, 0, "facePos"));
pi2 = mxGetData(mxGetField(cells, 0, "faces" ));
for (i = 0; i < nc + 1; i++) { pi1[i] = 0; }
/* 1) Count connections (i.e., faces per cell). */
for (f = 0; f < nf; f++) {
c1 = grid->face_neighbors[2*f + 0];
c2 = grid->face_neighbors[2*f + 1];
if (c1 >= 0) { pi1[c1 + 1] += 1; }
if (c2 >= 0) { pi1[c2 + 1] += 1; }
}
/* 2) Define start pointers (really, position *end* pointers at start). */
for (c = 1; c <= nc; c++) {
pi1[0] += pi1[c];
pi1[c] = pi1[0] - pi1[c];
}
/* 3) Fill connection structure whilst advancing end pointers. */
nhf = pi1[0];
pi1[0] = 0;
mxAssert (((size_t) nhf) == mxGetM(mxGetField(cells, 0, "faces")),
"Number of half faces (SIZE(cells.faces,1)) incorrectly "
"determined earlier.");
for (f = 0; f < nf; f++) {
cf_tag = 2*grid->face_tag[f] + 1; /* [1, 3, 5] */
c1 = grid->face_neighbors[2*f + 0];
c2 = grid->face_neighbors[2*f + 1];
if (c1 >= 0) {
pi2[ pi1[ c1 + 1 ] + 0*nhf ] = f + 1;
pi2[ pi1[ c1 + 1 ] + 1*nhf ] = cf_tag + 1; /* out */
pi1[ c1 + 1 ] += 1;
}
if (c2 >= 0) {
pi2[ pi1[ c2 + 1 ] + 0*nhf ] = f + 1;
pi2[ pi1[ c2 + 1 ] + 1*nhf ] = cf_tag + 0; /* in */
pi1[ c2 + 1 ] += 1;
}
}
/* Finally, adjust pointer array for one-based indexing in M. */
for (i = 0; i < nc + 1; i++) { pi1[i] += 1; }
/* Fill cells.indexMap. Note that despite the name, 'local_cell_index'
* really *is* the (zero-based) indexMap of the 'processed_grid'. */
pi1 = mxGetData(mxGetField(cells, 0, "indexMap"));
for (c = 0; c < nc; c++) { pi1[c] = grid->local_cell_index[c] + 1; }
}
/* ---------------------------------------------------------------------- */
static mxArray *
allocate_grid(struct processed_grid *grid, const char *func)
/* ---------------------------------------------------------------------- */
{
size_t nflds, nhf;
const char *fields[] = { "nodes", "faces", "cells",
"type", "cartDims", "griddim" };
mxArray *G, *nodes, *faces, *cells;
mxArray *type, *typestr, *cartDims, *griddim;
nflds = sizeof(fields) / sizeof(fields[0]);
nhf = count_halffaces(grid->number_of_faces, grid->face_neighbors);
G = mxCreateStructMatrix(1, 1, nflds, fields);
nodes = allocate_nodes(grid->number_of_nodes);
faces = allocate_faces(grid->number_of_faces,
grid->face_ptr[ grid->number_of_faces ]);
cells = allocate_cells(grid->number_of_cells, nhf);
type = mxCreateCellMatrix(1, 1);
typestr = mxCreateString(func);
cartDims = mxCreateDoubleMatrix(1, 3, mxREAL);
griddim = mxCreateDoubleScalar(3);
if ((G != NULL) && (nodes != NULL) && (faces != NULL) &&
(cells != NULL) && (type != NULL) && (typestr != NULL) &&
(cartDims != NULL) && (griddim != NULL)) {
mxSetCell(type, 0, typestr);
mxSetField(G, 0, "nodes" , nodes );
/* faces */
const char *n3[] = {"num", "neighbors", "nodes", "numNodes", "nodePos", "tag"};
mxArray *faces = mxCreateStructMatrix(1,1,sizeof n3 / sizeof n3[0], n3);
mxSetField(faces, 0, "num", mxCreateDoubleScalar(grid->number_of_faces));
mxArray *faceneighbors = mxCreateNumericMatrix(grid->number_of_faces, 2,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(faceneighbors);
for(j=0; j<2; ++j)
{
for (i=0; i<grid->number_of_faces; ++i)
{
int ix = grid->face_neighbors[2*i+j];
if (ix == -1)
{
iptr[i+grid->number_of_faces*j] = 0;
}
else
{
iptr[i+grid->number_of_faces*j] = ix+1;
}
}
}
}
mxSetField(faces, 0, "neighbors", faceneighbors);
mxArray *numnodes = mxCreateNumericMatrix(grid->number_of_faces, 1,
mxINT32_CLASS, mxREAL);
mxArray *nodepos = mxCreateNumericMatrix(grid->number_of_faces+1, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr1 = mxGetData(numnodes);
int *iptr2 = mxGetData(nodepos);
iptr2[0] = 1;
for (i=0; i<grid->number_of_faces; ++i)
{
iptr1[i] = grid->face_ptr[i+1] - grid->face_ptr[i];
iptr2[i+1] = iptr2[i] + iptr1[i];
}
}
mxSetField(faces, 0, "numNodes", numnodes);
mxSetField(faces, 0, "nodePos", nodepos);
mxArray *tags = mxCreateNumericMatrix(grid->number_of_faces, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(tags);
for (i = 0; i < grid->number_of_faces; ++i)
{
iptr[i] = grid->face_tag[i] + 1;
}
}
mxSetField(faces, 0, "tag", tags);
const char *n4[] = {"num", "faces", "facePos", "indexMap"};
mxArray *cells = mxCreateStructMatrix(1,1,sizeof n4 / sizeof n4[0], n4);
mxSetField(cells, 0, "num", mxCreateDoubleScalar(grid->number_of_cells));
mxArray *map = mxCreateNumericMatrix(grid->number_of_cells, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(map);
for(i=0; i<grid->number_of_cells; ++i)
{
iptr[i] = grid->local_cell_index[i]+1;
}
}
mxSetField(cells, 0, "indexMap", map);
mxArray *facepos = mxCreateNumericMatrix(grid->number_of_cells+1, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(facepos);
for(i=0; i<grid->number_of_cells+1; ++i)
{
iptr[i] = 0;
}
for (i=0; i<2*grid->number_of_faces; ++i)
{
int c=grid->face_neighbors[i];
if(c != -1)
{
iptr[c+1]++;
}
}
iptr[0] = 1;
for(i=0; i<grid->number_of_cells; ++i)
{
iptr[i+1] += iptr[i];
}
}
mxSetField(cells, 0, "facePos", facepos);
int *counter = calloc(grid->number_of_cells, sizeof(*counter));
int num_half_faces = 0;
{
int *iptr = mxGetData(facepos);
for(i=0; i<grid->number_of_cells; ++i)
{
counter[i] = num_half_faces;
num_half_faces += iptr[i+1]-iptr[i];
}
}
mxArray *cellfaces = mxCreateNumericMatrix(num_half_faces, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(cellfaces);
for (i=0; i<grid->number_of_faces; ++i)
{
int c1 = grid->face_neighbors[2*i];
int c2 = grid->face_neighbors[2*i+1];
if(c1 != -1) iptr[counter[c1]++] = i+1;
if(c2 != -1) iptr[counter[c2]++] = i+1;
}
}
mxSetField(cells, 0, "faces", cellfaces);
mxSetField(G, 0, "cells", cells);
int n = grid->face_ptr[grid->number_of_faces];
mxArray *facenodes = mxCreateNumericMatrix(n, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(facenodes);
for (i=0; i<n; ++i)
{
iptr[i] = grid->face_nodes[i]+1;
}
}
mxSetField(faces, 0, "nodes", facenodes);
mxSetField(G, 0, "faces" , faces );
free(counter);
/* cartDims */
mxArray *cartDims = mxCreateDoubleMatrix(1, 3, mxREAL);
ptr = mxGetPr(cartDims);
ptr[0] = grid->dimensions[0];
ptr[1] = grid->dimensions[1];
ptr[2] = grid->dimensions[2];
mxSetField(G, 0, "cartDims", cartDims);
/* type */
mxArray *type = mxCreateCellMatrix(1, 1);
mxSetCell(type, 0, mxCreateString("processgrid"));
mxSetField(G, 0, "cells" , cells );
mxSetField(G, 0, "type" , type );
mxSetField(G, 0, "cartDims", cartDims);
mxSetField(G, 0, "griddim" , griddim );
} else {
if (griddim != NULL) { mxDestroyArray(griddim); }
if (cartDims != NULL) { mxDestroyArray(cartDims); }
if (typestr != NULL) { mxDestroyArray(typestr); }
if (type != NULL) { mxDestroyArray(type); }
if (cells != NULL) { mxDestroyArray(cells); }
if (faces != NULL) { mxDestroyArray(faces); }
if (nodes != NULL) { mxDestroyArray(nodes); }
if (G != NULL) { mxDestroyArray(G); }
out[0] = G;
G = NULL;
}
return G;
}
/* Gateway routine for Matlab mex function. */
/*-------------------------------------------------------*/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
static void
fill_grid(mxArray *G, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
/* Set up data passed from Matlab */
struct grdecl g;
struct processed_grid out;
double tolerance = 0.0;
double *pr;
mx_init_grdecl(&g, prhs[0]);
if (nrhs == 2)
pr = mxGetPr(mxGetField(G, 0, "cartDims"));
pr[0] = grid->dimensions[0];
pr[1] = grid->dimensions[1];
pr[2] = grid->dimensions[2];
fill_nodes(mxGetField(G, 0, "nodes"), grid);
fill_faces(mxGetField(G, 0, "faces"), grid);
fill_cells(mxGetField(G, 0, "cells"), grid);
}
/* ---------------------------------------------------------------------- */
static int
args_ok(int nlhs, int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
{
tolerance = mxGetScalar (prhs[1]);
int ok;
ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2));
ok = ok && !mxIsEmpty(prhs[0]);
ok = ok && mxIsStruct(prhs[0]);
ok = ok && (mxGetFieldNumber(prhs[0], "cartDims") >= 0);
ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0);
ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0);
if (ok && (nrhs == 2)) {
ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1);
}
process_grdecl(&g, tolerance, &out);
return ok;
}
if (plhs >0)
/* ---------------------------------------------------------------------- */
static double
define_tolerance(int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
{
/* write to matlab struct */
fill_grid(plhs, &out);
double tol;
tol = 0.0;
if (nrhs == 2) {
tol = mxGetScalar(prhs[1]);
}
return tol;
}
/* Free whatever was allocated in initGrdecl. */
free_processed_grid(&out);
/* G = processgrid(grdecl)
G = processgrid(grdecl, tolerance)
*/
/* ---------------------------------------------------------------------- */
void
mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
{
double tolerance;
char errmsg[1023 + 1];
struct grdecl grdecl;
struct processed_grid g;
if (args_ok(nlhs, nrhs, prhs)) {
mx_init_grdecl(&grdecl, prhs[0]);
tolerance = define_tolerance(nrhs, prhs);
process_grdecl(&grdecl, tolerance, &g);
plhs[0] = allocate_grid(&g, mexFunctionName());
if (plhs[0] != NULL) {
fill_grid(plhs[0], &g);
} else {
/* Failed to create grid structure. Return empty. */
plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL);
}
free_processed_grid(&g);
} else {
sprintf(errmsg,
"Calling sequence is\n\t"
"G = %s(grdecl)\t%%or\n\t"
"G = %s(grdecl, tolerance)\n"
"The 'grdecl' must be a valid structure with fields\n"
"\t'cartDims', 'COORD', 'ZCORN'",
mexFunctionName(), mexFunctionName());
mexErrMsgTxt(errmsg);
}
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

41
opm/core/grid/cpgpreprocess/processgrid.m
View File

@ -1,12 +1,14 @@
function varargout = processgrid(varargin)
function G = processgrid(varargin)
%Compute grid topology and geometry from pillar grid description.
%
% SYNOPSIS:
% G = processGRDECL(grdecl)
% G = processgrid(grdecl)
% G = processgrid(grdecl,ztol)
%
% PARAMETERS:
% grdecl - Raw pillar grid structure, as defined by function
% 'readGRDECL', with fields COORDS, ZCORN and, possibly, ACTNUM.
% ztol - tolerance for unique points
%
% RETURNS:
% G - Valid grid definition containing connectivity, cell
@ -32,10 +34,33 @@ function varargout = processgrid(varargin)
% $Date$
% $Revision$
% Build MEX edition of same.
%
buildmex CFLAGS='$CFLAGS -Wall -fPIC' processgrid.c preprocess.c ...
uniquepoints.c facetopology.c sparsetable.c mxgrdecl.c
G = processgrid_mex(varargin{:});
G.griddim = 3;
G = splitDisconnectedGrid(G, false);
end
% Call MEX edition.
[varargout{1:nargout}] = processgrid(varargin{:});
function G = splitDisconnectedGrid(G, verbose)
% Check if grid is connected
ix = all(G.faces.neighbors~=0, 2);
I = [G.faces.neighbors(ix,1);G.faces.neighbors(ix,2)];
J = [G.faces.neighbors(ix,2);G.faces.neighbors(ix,1)];
N = double(max(G.faces.neighbors(:)));
A = sparse(double(I),double(J),1,N,N)+speye(N);
clear ix I J
[a,b,c,d]=dmperm(A); %#ok
clear A b d
if numel(c) > 2,
dispif(verbose, '\nGrid has %d disconnected components\n', ...
numel(c)- 1);
% Partition grid into connected subgrids
for i = 1:numel(c) - 1,
g(i) = extractSubgrid(G, a(c(i):c(i+1)-1)); %#ok
sz(i) = g(i).cells.num; %#ok
g(i).cartDims = G.cartDims; %#ok
end
% Return largest (in number of cells) grid first
[i,i] = sort(-sz); %#ok
G = g(i);
end
end

282
opm/core/grid/cpgpreprocess/readvector.cpp
View File

@ -1,282 +0,0 @@
/*===========================================================================
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <assert.h>
#include <vector>
#include "readvector.hpp"
struct ParserState {
int error;
};
template <typename T>
static void parse_next_token (FILE *fp, struct ParserState *s,
T *v, int *count);
static char *get_token_string (FILE *fp, struct ParserState *s, char *buf);
static void skip_line (FILE *fp, struct ParserState *s);
static char *convert_string (int *value, char *str);
static char *convert_string (double *value, char *str);
static void split_token_string(char *str, char *rstr, char *vstr);
static int get_repeat_count (char *rstr, struct ParserState *s);
template <typename T>
static void get_value (char *vstr, struct ParserState *s, T *v);
/* ------------------------------------------------------------------ */
static const int bufsz = 1024;
/* ------------------------------------------------------------------ */
template <typename T>
static void readvector(FILE *fp, struct ParserState *s, std::vector<T>& v)
{
T value;
int count, i;
count = 1;
s->error = 0;
while (count > 0){
parse_next_token(fp, s, &value, &count);
for (i=0; i<count; ++i)
{
v.push_back(value);
}
}
if (s->error) {
assert(0);
}
}
/* ------------------------------------------------------------------ */
template <typename T>
static void parse_next_token(FILE *fp, struct ParserState *s, T *v, int *r)
{
char str[bufsz];
char rstr[bufsz];
char vstr[bufsz];
str[0]='\0';
while (strlen(str)==0)
{
get_token_string(fp, s, str);
}
if (str[0] == '/')
{
/* end-of-vector marker */
*r = 0;
}
else
{
split_token_string(str, rstr, vstr);
*r = get_repeat_count(rstr, s);
get_value(vstr, s, v);
if (s->error)
{
/* signal abort to caller */
*r = 0;
}
}
}
/* ------------------------------------------------------------------ */
/* split string 'rrrr*vvvv' in strings 'rrrr' and 'vvvv' */
static void split_token_string(char *str, char *rstr, char *vstr)
{
char *ptr;
ptr=strchr(str, '*');
if ((ptr != NULL) && (ptr!=str)) {
while(str!=ptr)
{
*rstr++=*str++;
}
*rstr='\0';
str++;
}
else
{
rstr[0]='\0';
}
while((*vstr++=*str++) )
{
continue ;
}
*vstr='\0';
}
/* ------------------------------------------------------------------ */
static int get_repeat_count(char *rstr, struct ParserState *s)
{
int r;
char *ptr;
if (strlen(rstr)>0)
{
ptr=convert_string(&r, rstr);
if ((ptr==rstr) || (strlen(ptr)>0) || (r<1))
{
s->error = 1;
}
}
else
{
r = 1;
}
return r;
}
/* ------------------------------------------------------------------ */
template <typename T>
static void get_value(char *vstr, struct ParserState *s, T *v)
{
char *ptr;
if (strlen(vstr)>0)
{
/* Convert FORTRAN style floats 3.13D-6 to IEEE */
for(ptr=vstr; *ptr; ++ptr)
{
if (*ptr=='D')
{
*ptr='E';
}
}
ptr = convert_string(v, vstr);
if ((ptr==vstr) || (strlen(ptr)>0))
{
s->error = 1;
}
}
else
{
s->error = 1;
}
}
/* ------------------------------------------------------------------ */
static char *get_token_string(FILE *fp, struct ParserState *s, char *str)
{
char *ptr = str;
int c;
/* Skip leading blanks */
while((c = fgetc(fp)) != EOF && isspace(c))
{
;
}
*ptr++ = c;
/* Catch end marker */
if (c == '/') {
skip_line(fp, s);
*ptr++ = '\0';
return str;
}
while((c = fgetc(fp)) != EOF && !isspace(c))
{
*ptr++ = c;
/* Break and skip rest of line if '--' if encountered */
if (*(ptr-1) == '-' && ptr-str>1 && *(ptr-2) == '-'){
ptr = ptr - 2;
skip_line(fp, s);
break;
}
/* If end marker is encontered, push character back onto stream. */
if (c=='/') {
ungetc(c, fp);
ptr--;
break;
}
assert(ptr-str < bufsz);
}
*ptr='\0';
return str;
}
/* ------------------------------------------------------------------ */
static void skip_line(FILE *fp, struct ParserState *s)
{
static_cast<void>(s);
int c;
while((c = fgetc(fp))!=EOF && c != '\n') {
;
}
}
/* ------------------------------------------------------------------ */
// int version
static char *convert_string(int *value, char *str)
{
char *q;
*value = strtol(str, &q, 10);
return q;
}
/* ------------------------------------------------------------------ */
// double version
static char *convert_string(double *value, char *str)
{
char *q;
*value = strtod(str, &q);
return q;
}
/* ------------------------------------------------------------------ */
template <typename T>
void read_vector_from_file(const std::string& fn, std::vector<T>& v)
{
FILE *fp = fopen(fn.c_str(), "r");
struct ParserState s = { 0 };
readvector(fp, &s, v);
fclose(fp);
}
void read_vector_from_file(const std::string& fn, std::vector<int>& v)
{
read_vector_from_file<int>(fn, v);
}
void read_vector_from_file(const std::string& fn, std::vector<double>& v)
{
read_vector_from_file<double>(fn, v);
}

24
opm/core/grid/cpgpreprocess/readvector.hpp
View File

@ -1,24 +0,0 @@
/*===========================================================================
//
// File: readvector.hpp
//
// Created: 2011-11-30 09:35:14+0100
//
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
*/
#ifndef READVECTOR_HPP_HEADER
#define READVECTOR_HPP_HEADER
#include <string>
void read_vector_from_file(const std::string&, std::vector<int>& v);
void read_vector_from_file(const std::string&, std::vector<double>& v);
#endif /* READVECTOR_HPP_HEADER */

100
opm/core/grid/cpgpreprocess/sparsetable.c
View File

@ -1,100 +0,0 @@
/*===========================================================================
//
// File: sparsetable.c
//
// Created: Fri Jun 19 08:48: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 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 <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <stdio.h>
#include "sparsetable.h"
void free_sparse_table (sparse_table_t *tab)
{
if(tab->ptr) free(tab->ptr);
if(tab->data) free(tab->data);
free(tab);
}
sparse_table_t *malloc_sparse_table(int m, int n, int datasz)
{
size_t alloc_sz;
sparse_table_t *tab = malloc(sizeof *tab);
tab->m = m;
tab->n = n;
tab->position = 0;
if (!(tab->ptr = malloc((m+1) * sizeof (*tab->ptr)))){
fprintf(stderr, "Could not allocate space for sparse ptr\n");
free_sparse_table(tab);
return NULL;
}
alloc_sz = datasz;
alloc_sz *= n;
if(!(tab->data = malloc(alloc_sz))){
fprintf(stderr, "Could not allocate space for sparse data(%d)\n", n);
free_sparse_table(tab);
return NULL;
}
return tab;
}
sparse_table_t *realloc_sparse_table(sparse_table_t *tab, int m, int n, int datasz)
{
void *p = realloc(tab->ptr, (m+1) * sizeof (*tab->ptr));
if (p){
tab->ptr = p;
}else{
fprintf(stderr, "Could not reallocate space for sparse ptr\n");
free_sparse_table(tab);
return NULL;
}
p = realloc(tab->data, n * datasz);
if(p){
tab->data = p;
}else{
fprintf(stderr, "Could not reallocate space for sparse data(%d)\n", n);
free_sparse_table(tab);
return NULL;
}
tab->m = m;
tab->n = n;
return tab;
}

53
opm/core/grid/cpgpreprocess/sparsetable.h
View File

@ -1,53 +0,0 @@
/*===========================================================================
//
// File: sparsetable.h
//
// Created: Fri Jun 19 08:47:45 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 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/>.
*/
#ifndef OPENRS_SPARSETABLE_HEADER
#define OPENRS_SPARSETABLE_HEADER
typedef struct{
int m; /* number of rows */
int *ptr; /* row pointer of size m+1 */
int position; /* first position in ptr that is not filled. */
int n; /* size of data */
void *data; /* sparse table data */
} sparse_table_t;
void free_sparse_table (sparse_table_t *tab);
sparse_table_t *malloc_sparse_table (int m, int n, int datasz);
sparse_table_t *realloc_sparse_table (sparse_table_t *tab, int m, int n, int datasz);
#endif /* OPENRS_SPARSETABLE_HEADER */

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

@ -33,35 +33,35 @@
*/
#include <assert.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <limits.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "sparsetable.h"
#include "preprocess.h"
#include "uniquepoints.h"
#define min(i,j) ((i)<(j) ? (i) : (j))
#define max(i,j) ((i)>(j) ? (i) : (j))
#define overlap(a1,a2,b1,b2) max(a1,b1) < min(a2,b2)
#define MIN(i,j) (((i) < (j)) ? (i) : (j))
#define MAX(i,j) (((i) > (j)) ? (i) : (j))
/*-----------------------------------------------------------------
Compare function passed to qsort
*/
Compare function passed to qsortx */
static int compare(const void *a, const void *b)
{
if (*(const double*)a < *(const double*) b) return -1;
else return 1;
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==1
*/
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[])
{
@ -81,17 +81,17 @@ static int createSortedList(double *list, int n, int m,
/*-----------------------------------------------------------------
Remove points less than <tolerance> apart in <list> of increasing
doubles.
*/
doubles. */
static int uniquify(int n, double *list, double tolerance)
{
int i, pos = 0;
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){
@ -120,8 +120,7 @@ static int uniquify(int n, double *list, double tolerance)
/*-----------------------------------------------------------------
Along single pillar:
*/
Along single pillar: */
static int assignPointNumbers(int begin,
int end,
const double *zlist,
@ -149,7 +148,8 @@ static int assignPointNumbers(int begin,
/* Skip inactive cells */
if (!a[i/2]) {
p[0] = p[-1]; /* Inactive cells are collapsed leaving void space.*/
p[0] = p[-1]; /* Inactive cells are collapsed leaving
* void space.*/
++p;
continue;
}
@ -175,26 +175,25 @@ static int assignPointNumbers(int begin,
}
/*-----------------------------------------------------------------
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[])
/* ---------------------------------------------------------------------- */
static void
vector_positions(const int dims[3] ,
const int i ,
const int j ,
size_t start[4])
/* ---------------------------------------------------------------------- */
{
size_t im, ip, jm, jp;
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;
im = MAX(1, i ) - 1;
jm = MAX(1, j ) - 1;
ip = MIN(dims[0], i+1) - 1;
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);
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);
}
@ -202,21 +201,35 @@ static void igetvectors(const int dims[3], int i, int j,
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.
*/
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];
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;
vector_positions(dims, i, j, start);
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);
for (p = 0; p < 4; p++) {
v[p] = field + start[p];
}
}
@ -226,10 +239,18 @@ static void dgetvectors(const int dims[3], int i, int j,
*/
static void interpolate_pillar(const double *coord, double *pt)
{
double a = (pt[2]-coord[2])/(coord[5]-coord[2]);
if (isinf(a) || isnan(a)){
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]);
@ -240,34 +261,34 @@ static void interpolate_pillar(const double *coord, double *pt)
}
/*-----------------------------------------------------------------
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)
*/
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*/
int *plist, /* list of point numbers on
* each pillar*/
double tolerance,
struct processed_grid *out)
{
int nx = out->dimensions[0];
int ny = out->dimensions[1];
int nz = out->dimensions[2];
/* ztab->data may need extra space temporarily due to simple boundary treatement */
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);
sparse_table_t *ztab = malloc_sparse_table(npillars,
npillarpoints,
sizeof(double));
double *zlist = malloc(npillarpoints*sizeof *zlist);
int *zptr = malloc((npillars+1)*sizeof *zptr);
int nc = g->dims[0]*g->dims[1]*g->dims[2];
double *zlist = ztab->data; /* casting void* to double* */
int *zptr = ztab->ptr;
int i,j,k;
@ -275,14 +296,21 @@ int finduniquepoints(const struct grdecl *g,
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;
double *pt;
int *p;
d1[0] = 2*g->dims[0]; d1[1] = 2*g->dims[1]; d1[2] = 2*g->dims[2];
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;
@ -291,9 +319,6 @@ int finduniquepoints(const struct grdecl *g,
for (j=0; j < g->dims[1]+1; ++j){
for (i=0; i < g->dims[0]+1; ++i){
const int *a[4];
const double *z[4];
/* 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);
@ -308,7 +333,8 @@ int finduniquepoints(const struct grdecl *g,
pt += 3;
}
/* Increment pointer to sparse table of unique zcorn values */
/* Increment pointer to sparse table of unique zcorn
* values */
zout = zout + len;
zptr[pos++] = zout - zlist;
@ -318,25 +344,25 @@ int finduniquepoints(const struct grdecl *g,
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 */
/* 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 */
int pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2);
pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2);
/* cell column position */
int cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]);
cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]);
/* zcorn column position */
int zix = 2*g->dims[2]*(i+2*g->dims[0]*j);
const int *a = g->actnum + cix;
const double *z = g->zcorn + zix;
zix = 2*g->dims[2]*(i+2*g->dims[0]*j);
if (!assignPointNumbers(zptr[pix], zptr[pix+1], zlist,
2*g->dims[2], z, a, p, tolerance)){
2*g->dims[2],
g->zcorn + zix, g->actnum + cix,
p, tolerance)){
fprintf(stderr, "Something went wrong in assignPointNumbers");
return 0;
}
@ -345,9 +371,8 @@ int finduniquepoints(const struct grdecl *g,
}
}
free_sparse_table(ztab);
free(zptr);
free(zlist);
return 1;
}
@ -355,4 +380,3 @@ int finduniquepoints(const struct grdecl *g,
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

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

@ -32,15 +32,15 @@ You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_UNIQUEPOINTS_HEADER
#define OPENRS_UNIQUEPOINTS_HEADER
#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 /* OPENRS_UNIQUEPOINTS_HEADER */
#endif /* OPM_UNIQUEPOINTS_HEADER */
/* Local Variables: */
/* c-basic-offset:4 */

10
opm/core/linalg/sparse_sys.h
View File

@ -37,13 +37,13 @@ extern "C" {
*/
struct CSRMatrix
{
size_t m; /** Number of rows */
size_t nnz; /** Number of structurally non-zero elements */
size_t m; /**< Number of rows */
size_t nnz; /**< Number of structurally non-zero elements */
int *ia; /** Row pointers */
int *ja; /** Column indices */
int *ia; /**< Row pointers */
int *ja; /**< Column indices */
double *sa; /** Matrix elements */
double *sa; /**< Matrix elements */
};

2
tests/Makefile.am
View File

@ -52,8 +52,6 @@ test_lapack_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_readpolymer_SOURCES = test_readpolymer.cpp
test_read_vag_SOURCES = test_read_vag.cpp
test_readvector_SOURCES = test_readvector.cpp
test_sf2p_SOURCES = test_sf2p.cpp
test_writeVtkData_SOURCES = test_writeVtkData.cpp

67
tests/test_readvector.cpp
View File

@ -1,67 +0,0 @@
#include <iostream>
#include <string>
#include <vector>
#include <opm/core/grid.h>
#include <opm/core/grid/cpgpreprocess/preprocess.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/grid/cpgpreprocess/readvector.hpp>
static struct UnstructuredGrid*
read_grid(const std::string& dir)
{
std::string fn;
fn = dir + '/' + "zcorn.txt";
std::vector<double> zcorn;
read_vector_from_file(fn, zcorn);
fn = dir + '/' + "coord.txt";
::std::vector<double> coord;
read_vector_from_file(fn, coord);
fn = dir + '/' + "actnum.txt";
std::vector<int> actnum;
read_vector_from_file(fn, actnum);
fn = dir + '/' + "dimens.txt";
::std::vector<int> dimens;
read_vector_from_file(fn, dimens);
struct grdecl grdecl;
grdecl.zcorn = &zcorn[0];
grdecl.coord = &coord[0];
grdecl.actnum = &actnum[0];
grdecl.dims[0] = dimens[0];
grdecl.dims[1] = dimens[1];
grdecl.dims[2] = dimens[2];
struct UnstructuredGrid *g = create_grid_cornerpoint(&grdecl, 0.0);
double vol = 0.0;
for (int c = 0; c < g->number_of_cells; c++) {
vol += g->cell_volumes[c];
}
std::cout << "Sum volumes = " << vol << '\n';
for (int c = 0, i = 0; c < g->number_of_cells; c++) {
for (; i < g->cell_facepos[c + 1]; i++) {
std::cout << "(c,i) = (" << c << "," << g->cell_facetag[i] << ")\n";
}
}
return g;
}
int main()
{
struct UnstructuredGrid *g;
g = read_grid(std::string("example"));
destroy_grid(g);
return 0;
}