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

Re-indent to four-spaces. While at it, do whitespace-cleanup.

Browse Source 
Signed-off-by: Bård Skaflestad <Bard.Skaflestad@sintef.no>
This commit is contained in:
Jostein R. Natvig
2012-01-31 08:48:40 +00:00
committed by Bård Skaflestad
parent 53bf30cb43
commit 2be1ecf16c
9 changed files with 1246 additions and 1233 deletions
Download Patch File Download Diff File Expand all files Collapse all files

508
facetopology.c
View File

@@ -13,23 +13,23 @@
//===========================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
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 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.
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/>.
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>
@@ -61,132 +61,132 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
/* 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)
int *a2,
int *b1,
int *b2,
int intersect[4],
int *faces)
{
int mask[8] = {-1};
int k;
int *f;
/* 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]; }
int mask[8] = {-1};
int k;
int *f;
/* 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");
}
/* 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];
/* Partial pinch of face */
if (mask[0] == mask[6]){
mask[6] = -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];
if (mask[2] == mask[4]){
mask[4] = -1;
}
else{
mask[2] = -1;
mask[4] = -1;
mask[3] = intersect[2];
/* 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];
}
}
}
f = faces;
for (k=7; k>=0; --k){
if(mask[k] != -1){
*f++ = mask[k];
/* 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");
}
/* 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");
fprintf(stderr, "\n");
#endif
return f;
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.
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.
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)
{
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]);
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]);
}
@@ -197,156 +197,158 @@ static int faceintersection(int *a1, int *a2, int *b1, int *b2)
/* 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)
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];
/* 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;
/* 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 k1 = 0;
int k2 = 0;
int i,j=0;
int intersect[4]= {-1};
int *tmp;
/* for (i=0; i<2*n; work[i++]=-1); */
for (i = 0; i<n-1; ++i){
/* pinched a-cell */
if (a1[i] == a1[i+1] && a2[i] == a2[i+1]){
continue;
int i,j=0;
int intersect[4]= {-1};
int *tmp;
/* for (i=0; i<2*n; work[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 (lineintersection(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);
}
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 (lineintersection(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: */

36
facetopology.h
View File

@@ -13,35 +13,35 @@
//===========================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
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 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.
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/>.
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],
int *intersectionlist,
int *intersectionlist,
int *work,
struct processed_grid *out);
struct processed_grid *out);
#endif /* OPENRS_FACETOPOLOGY_HEADER */
#endif /* OPM_FACETOPOLOGY_HEADER */
/* Local Variables: */
/* c-basic-offset:4 */

18
grdecl.h
View File

@@ -1,12 +1,16 @@
#ifndef GRDECL_H
#define GRDECL_H
#ifndef GRDECL_H_INCLUDED
#define GRDECL_H_INCLUDED
struct grdecl{
int dims[3];
const double *coord;
const double *zcorn;
const int *actnum;
int dims[3];
const double *coord;
const double *zcorn;
const int *actnum;
};
#endif
#endif /* GRDECL_H_INCLUDED */
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

132
mxgrdecl.c
View File

@@ -13,23 +13,23 @@
//=======================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
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 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.
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/>.
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>
@@ -48,60 +48,64 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s);
/*-------------------------------------------------------*/
void mx_init_grdecl(struct grdecl *g, const mxArray *s)
{
int i,n;
mxArray *field;
int numel;
double *tmp;
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";
mexErrMsgTxt(str);
}
numel = mxGetNumberOfElements(cartdims);
if (!mxIsNumeric(cartdims) || numel != 3){
mexErrMsgTxt("cartDims field must be 3 numbers");
}
tmp = mxGetPr(cartdims);
n = 1;
for (i=0; i<3; ++i){
g->dims[i] = tmp[i];
n *= tmp[i];
}
int i,n;
mxArray *field;
int numel;
double *tmp;
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";
mexErrMsgTxt(str);
}
numel = mxGetNumberOfElements(actnum);
if (mxGetClassID(actnum) != mxINT32_CLASS ||
numel != g->dims[0]*g->dims[1]*g->dims[2] ){
mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32");
}
g->actnum = mxGetData(actnum);
numel = mxGetNumberOfElements(cartdims);
if (!mxIsNumeric(cartdims) || numel != 3){
mexErrMsgTxt("cartDims field must be 3 numbers");
}
tmp = mxGetPr(cartdims);
n = 1;
for (i=0; i<3; ++i){
g->dims[i] = tmp[i];
n *= tmp[i];
}
field = mxGetField(s, 0, "COORD");
numel = mxGetNumberOfElements(coord);
if (mxGetClassID(coord) != mxDOUBLE_CLASS ||
numel != 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(actnum);
if (mxGetClassID(actnum) != mxINT32_CLASS ||
numel != g->dims[0]*g->dims[1]*g->dims[2] ){
mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32");
}
g->actnum = mxGetData(actnum);
numel = mxGetNumberOfElements(zcorn);
if (mxGetClassID(zcorn) != mxDOUBLE_CLASS ||
numel != 8*g->dims[0]*g->dims[1]*g->dims[2]){
mexErrMsgTxt("ZCORN field must have 8*nx*ny*nz doubles.");
}
g->zcorn = mxGetPr(zcorn);
field = mxGetField(s, 0, "COORD");
numel = mxGetNumberOfElements(coord);
if (mxGetClassID(coord) != mxDOUBLE_CLASS ||
numel != 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]){
mexErrMsgTxt("ZCORN field must have 8*nx*ny*nz doubles.");
}
g->zcorn = mxGetPr(zcorn);
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

36
mxgrdecl.h
View File

@@ -13,30 +13,32 @@
//=======================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
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 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.
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/>.
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: */

752
preprocess.c
View File

@@ -13,23 +13,23 @@
//==========================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
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 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.
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/>.
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>
@@ -49,31 +49,30 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len);
int checkmemory(int nz, struct processed_grid *out, int **intersections);
void process_vertical_faces(int direction,
int **intersections,
int *plist, int *work,
struct processed_grid *out);
int **intersections,
int *plist, int *work,
struct processed_grid *out);
void process_horizontal_faces(int **intersections,
int *plist,
struct processed_grid *out);
int *plist,
struct processed_grid *out);
/*-----------------------------------------------------------------
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(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);
v[2] = field + dims[2]*(ip + dims[0]* jm);
v[3] = field + dims[2]*(ip + dims[0]* jp);
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);
}
@@ -84,269 +83,268 @@ 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.
*/
*/
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]){
for(k=0; k<len; k+=2){
neighbors[k] = -1;
int k;
if (i<0 || i>=dims[0] || j<0 || j >= dims[1]){
for(k=0; k<len; k+=2){
neighbors[k] = -1;
}
}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;
}
}
}
}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;
}
}
}
}
/*-----------------------------------------------------------------
Ensure there's sufficient memory
*/
Ensure there's sufficient memory */
int checkmemory(int nz, struct processed_grid *out, int **intersections)
{
/* 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 */
int r = (2*nz+2)*(2*nz+2);
int m = out->m;
int n = out->n;
if(out->number_of_faces + r > m){
m += max(m*0.5, 2*r);
}
if (out->face_ptr[out->number_of_faces] + 6*r > n){
n += max(n*0.5, 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;
if(out->number_of_faces + r > m){
m += max(m*0.5, 2*r);
}
}
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 (p1 && p2 && p3) {
out->face_nodes = p1;
out->face_ptr = p2;
out->face_tag = p3;
} else {
return 0;
if (out->face_ptr[out->number_of_faces] + 6*r > n){
n += max(n*0.5, 12*r);
}
out->m = m;
out->n = n;
}
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;
}
}
return 1;
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 (p1 && p2 && p3) {
out->face_nodes = p1;
out->face_ptr = p2;
out->face_tag = p3;
} else {
return 0;
}
out->m = m;
out->n = n;
}
return 1;
}
/*-----------------------------------------------------------------
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.
-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.
*/
*/
void process_vertical_faces(int direction,
int **intersections,
int *plist, int *work,
struct processed_grid *out)
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;
for (j=0; j<ny+direction; ++j) {
for (i=0; i<nx+1-direction; ++i){
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;
for (j=0; j<ny+direction; ++j) {
for (i=0; i<nx+1-direction; ++i){
if (!checkmemory(nz, out, intersections)){
fprintf(stderr, "Could not allocat enough space in process_vertical_faces\n");
exit(1);
}
if (!checkmemory(nz, out, intersections)){
fprintf(stderr, "Could not allocat 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);
/* 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);
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;
}
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;
/* 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;
findconnections(2*nz+2, cornerpts,
*intersections+4*num_intersections,
work, out);
findconnections(2*nz+2, cornerpts,
*intersections+4*num_intersections,
work, out);
ptr = out->face_neighbors + 2*startface;
len = 2*out->number_of_faces - 2*startface;
ptr = out->face_neighbors + 2*startface;
len = 2*out->number_of_faces - 2*startface;
compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr, len);
compute_cell_index(out->dimensions, i, j, ptr+1, len);
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];
}
/* Tag the new faces */
f = startface;
for (; f < out->number_of_faces; ++f) {
out->face_tag[f] = tag[direction];
}
}
}
}
}
static int linearindex(const int dims[3], int i, int j, int k)
{
return i+dims[0]*(j+dims[1]*k);
return i+dims[0]*(j+dims[1]*k);
}
/*-----------------------------------------------------------------
For each horizontal face (i.e. k constant),
-find point numbers for the corners and
-cell neighbors.
-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)
*/
*/
void process_horizontal_faces(int **intersections,
int *plist,
struct processed_grid *out)
int *plist,
struct processed_grid *out)
{
int i,j,k;
int i,j,k;
int nx = out->dimensions[0];
int ny = out->dimensions[1];
int nz = out->dimensions[2];
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;
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;
/* 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) {
for(j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
if (!checkmemory(nz, out, intersections)){
fprintf(stderr, "Could not allocat enough space in process_horizontal_faces\n");
exit(1);
}
if (!checkmemory(nz, out, intersections)){
fprintf(stderr, "Could not allocat 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;
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);
/* Vectors of point numbers */
igetvectors(d, 2*i+1, 2*j+1, plist, c);
prevcell = -1;
prevcell = -1;
for (k = 1; k<nz*2+1; ++k){
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]){
/* 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 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];
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;
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;
thiscell = linearindex(out->dimensions, i,j,(k-1)/2);
*n++ = prevcell;
*n++ = prevcell = thiscell;
cell[thiscell] = cellno++;
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];
}
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;
out->face_tag[ out->number_of_faces] = TOP;
out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
*n++ = prevcell;
*n++ = prevcell = -1;
}
}
}
}
*n++ = prevcell;
*n++ = prevcell = -1;
}
}
}
}
}
}
}
out->number_of_cells = cellno;
out->number_of_cells = cellno;
}
@@ -363,55 +361,55 @@ void process_horizontal_faces(int **intersections,
static void approximate_intersection_pt(int *L, double *c, double *pt)
{
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];
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];
double a = (z2-z0)/(z1-z0 - (z3-z2));
if (isinf(a) || isnan(a)){
a = 0;
}
double a = (z2-z0)/(z1-z0 - (z3-z2));
if (isinf(a) || isnan(a)){
a = 0;
}
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;
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;
}
/*-----------------------------------------------------------------
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)
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");
}
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;
/* 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;
for (k=np; k<n; ++k){
approximate_intersection_pt(itsct, out->node_coordinates, pt);
pt += 3;
itsct += 4;
}
}
}
@@ -425,9 +423,9 @@ copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *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,:)]
out = [in(0,0,:), in(1,0,:),..., in(nx-1, ny-1,:)]
in Matlab pseudo-code.
in Matlab pseudo-code.
*/
for (j=0; j<ny; ++j){
for (i=0; i<nx; ++i){
@@ -441,7 +439,7 @@ copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out)
/* ------------------------------------------------------------------ */
static double*
copy_and_permute_zcorn(int nx, int ny, int nz, const double *in,
copy_and_permute_zcorn(int nx, int ny, int nz, const double *in,
double sign, double *out)
/* ------------------------------------------------------------------ */
{
@@ -450,9 +448,9 @@ copy_and_permute_zcorn(int nx, int ny, int nz, const double *in,
/* 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,:)]
out = [in(0,0,:), in(1,0,:),..., in(2*nx-1, 2*ny-1,:)]
in Matlab pseudo-code.
in Matlab pseudo-code.
*/
for (j=0; j<2*ny; ++j){
for (i=0; i<2*nx; ++i){
@@ -466,7 +464,7 @@ copy_and_permute_zcorn(int nx, int ny, int nz, const double *in,
/* ------------------------------------------------------------------ */
static int
get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
const double *zcorn, int *error)
/* ------------------------------------------------------------------ */
{
@@ -494,13 +492,13 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
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[c1] && actnum[c2] && (z2 < z1)){
fprintf(stderr, "\nZCORN should be strictly "
"nondecreasing along pillars!\n");
"nondecreasing along pillars!\n");
*error = 1;
goto end;
}
@@ -525,171 +523,173 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
/*-----------------------------------------------------------------
Public interface
*/
*/
void process_grdecl(const struct grdecl *in,
double tolerance,
struct processed_grid *out)
double tolerance,
struct processed_grid *out)
{
struct grdecl g;
struct grdecl g;
int i;
int sign, error;
int cellnum;
int i;
int sign, error;
int cellnum;
int *actnum, *iptr;
int *global_cell_index;
int *actnum, *iptr;
int *global_cell_index;
double *zcorn;
double *zcorn;
const int BIGNUM = 64;
const int nx = in->dims[0];
const int ny = in->dims[1];
const int nz = in->dims[2];
const int nc = nx*ny*nz;
const int BIGNUM = 64;
const int nx = in->dims[0];
const int ny = in->dims[1];
const int nz = in->dims[2];
const int nc = nx*ny*nz;
/* internal work arrays */
int *work;
int *plist;
int *intersections;
/* 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 = 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 = BIGNUM/3;
out->n = 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->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] = nx;
out->dimensions[1] = ny;
out->dimensions[2] = nz;
out->number_of_faces = 0;
out->number_of_nodes = 0;
out->number_of_cells = 0;
out->dimensions[0] = nx;
out->dimensions[1] = ny;
out->dimensions[2] = nz;
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->node_coordinates = NULL;
out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index);
/* Do actual work here:*/
/* 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.*/
/* -----------------------------------------------------------------*/
/* 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] = nx;
g.dims[1] = ny;
g.dims[2] = nz;
/* initialize grdecl structure "g" that will be processd by
* "finduniquepoints" */
g.dims[0] = nx;
g.dims[1] = ny;
g.dims[2] = nz;
actnum = malloc (nc * sizeof *actnum);
g.actnum = copy_and_permute_actnum(nx, ny, nz, in->actnum, actnum);
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);
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;
g.coord = in->coord;
/* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */
plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist);
/* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX)
* padding */
plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist);
finduniquepoints(&g, plist, tolerance, out);
finduniquepoints(&g, plist, tolerance, out);
free (zcorn);
free (actnum);
free (zcorn);
free (actnum);
/* -----------------------------------------------------------------*/
/* Find face topology and face-to-cell connections */
/* -----------------------------------------------------------------*/
/* Find face topology and face-to-cell connections */
/* internal */
work = malloc(2* (2*nz+2)* sizeof(*work));
for(i=0; i<4*(nz+1); ++i) { work[i] = -1; }
/* internal */
work = malloc(2* (2*nz+2)* sizeof(*work));
for(i=0; i<4*(nz+1); ++i) { work[i] = -1; }
/* internal array to store intersections */
intersections = malloc(BIGNUM* sizeof(*intersections));
/* 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);
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);
free (plist);
free (work);
/* -----------------------------------------------------------------*/
/* (re)allocate space for and compute coordinates of nodes that
* arise from intesecting cells (faults) */
compute_intersection_coordinates(intersections, out);
/* -----------------------------------------------------------------*/
/* (re)allocate space for and compute coordinates of nodes that
* arise from intesecting cells (faults) */
compute_intersection_coordinates(intersections, out);
free (intersections);
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(out->number_of_cells *
sizeof (*global_cell_index));
cellnum = 0;
for (i=0; i<nx*ny*nz; ++i){
if(out->local_cell_index[i]!=-1){
global_cell_index[cellnum] = i;
out->local_cell_index[i] = cellnum;
cellnum++;
/* -----------------------------------------------------------------*/
/* 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(out->number_of_cells *
sizeof (*global_cell_index));
cellnum = 0;
for (i=0; i<nx*ny*nz; ++i){
if(out->local_cell_index[i]!=-1){
global_cell_index[cellnum] = i;
out->local_cell_index[i] = cellnum;
cellnum++;
}
}
}
/* Remap out->face_neighbors */
iptr=out->face_neighbors;
for (i=0; i<out->number_of_faces*2; ++i, ++iptr){
if (*iptr != -1){
*iptr = out->local_cell_index[*iptr];
/* Remap out->face_neighbors */
iptr=out->face_neighbors;
for (i=0; i<out->number_of_faces*2; ++i, ++iptr){
if (*iptr != -1){
*iptr = out->local_cell_index[*iptr];
}
}
}
free(out->local_cell_index);
out->local_cell_index = global_cell_index;
/* 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<3*out->number_of_nodes; i = i+3)
out->node_coordinates[i]=sign*out->node_coordinates[i];
}
free(out->local_cell_index);
out->local_cell_index = global_cell_index;
/* 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<3*out->number_of_nodes; i = i+3)
out->node_coordinates[i]=sign*out->node_coordinates[i];
}
}
/*-------------------------------------------------------*/
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 );
}
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: */

526
processgrid.c
View File

@@ -13,23 +13,23 @@
//===========================================================================*/
/*
Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2011, 2012 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).
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 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.
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/>.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
/* Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. */
@@ -51,30 +51,30 @@ static mxArray *
allocate_nodes(size_t nnodes)
/* ---------------------------------------------------------------------- */
{
size_t nflds;
const char *fields[] = { "num", "coords" };
size_t nflds;
const char *fields[] = { "num", "coords" };
mxArray *nodes, *num, *coords;
mxArray *nodes, *num, *coords;
nflds = sizeof(fields) / sizeof(fields[0]);
nflds = sizeof(fields) / sizeof(fields[0]);
nodes = mxCreateStructMatrix(1, 1, nflds, fields);
nodes = mxCreateStructMatrix(1, 1, nflds, fields);
num = mxCreateDoubleScalar(nnodes);
coords = mxCreateDoubleMatrix(nnodes, 3, mxREAL);
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); }
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;
}
nodes = NULL;
}
return nodes;
return nodes;
}
@@ -83,17 +83,17 @@ static void
fill_nodes(mxArray *nodes, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
size_t i, j, nnodes;
double *coords;
size_t i, j, nnodes;
double *coords;
nnodes = grid->number_of_nodes;
coords = mxGetPr(mxGetField(nodes, 0, "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];
for (j = 0; j < 3; ++j) {
for (i = 0; i < nnodes; ++i) {
*coords++ = grid->node_coordinates[3*i + j];
}
}
}
}
@@ -102,40 +102,40 @@ static mxArray *
allocate_faces(size_t nf, size_t nfacenodes)
/* ---------------------------------------------------------------------- */
{
size_t nflds;
const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" };
size_t nflds;
const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" };
mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag;
mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag;
nflds = sizeof(fields) / sizeof(fields[0]);
nflds = sizeof(fields) / sizeof(fields[0]);
faces = mxCreateStructMatrix(1, 1, nflds, fields);
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);
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); }
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;
}
faces = NULL;
}
return faces;
return faces;
}
@@ -144,33 +144,33 @@ static void
fill_faces(mxArray *faces, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
size_t i, f, nf, nfn;
size_t i, f, nf, nfn;
int *pi;
int *pi;
nf = grid->number_of_faces;
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.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.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.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; }
/* Fill faces.tag */
pi = mxGetData(mxGetField(faces, 0, "tag"));
for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[f] + 1; }
}
@@ -179,18 +179,18 @@ static size_t
count_halffaces(size_t nf, const int *neighbors)
/* ---------------------------------------------------------------------- */
{
int c1, c2;
size_t nhf, f;
int c1, c2;
size_t nhf, f;
for (f = nhf = 0; f < nf; f++) {
c1 = neighbors[2*f + 0];
c2 = neighbors[2*f + 1];
for (f = nhf = 0; f < nf; f++) {
c1 = neighbors[2*f + 0];
c2 = neighbors[2*f + 1];
nhf += c1 >= 0;
nhf += c2 >= 0;
}
nhf += c1 >= 0;
nhf += c2 >= 0;
}
return nhf;
return nhf;
}
@@ -199,37 +199,37 @@ static mxArray *
allocate_cells(size_t nc, size_t ncf)
/* ---------------------------------------------------------------------- */
{
size_t nflds;
const char *fields[] = { "num", "facePos", "faces", "indexMap" };
size_t nflds;
const char *fields[] = { "num", "facePos", "faces", "indexMap" };
mxArray *cells, *num, *facePos, *faces, *indexMap;
mxArray *cells, *num, *facePos, *faces, *indexMap;
nflds = sizeof(fields) / sizeof(fields[0]);
nflds = sizeof(fields) / sizeof(fields[0]);
cells = mxCreateStructMatrix(1, 1, nflds, fields);
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);
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); }
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;
}
cells = NULL;
}
return cells;
return cells;
}
@@ -238,69 +238,69 @@ static void
fill_cells(mxArray *cells, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
size_t c, nc, f, nf, i;
size_t c, nc, f, nf, i;
int c1, c2, cf_tag, nhf;
int *pi1, *pi2;
int c1, c2, cf_tag, nhf;
int *pi1, *pi2;
nc = grid->number_of_cells;
nf = grid->number_of_faces;
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; }
/* 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];
/* 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 (c1 >= 0) { pi1[c1 + 1] += 1; }
if (c2 >= 0) { pi1[c2 + 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;
/* 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];
}
}
/* Finally, adjust pointer array for one-based indexing in M. */
for (i = 0; i < nc + 1; i++) { pi1[i] += 1; }
/* 3) Fill connection structure whilst advancing end pointers. */
nhf = pi1[0];
pi1[0] = 0;
/* 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; }
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; }
}
@@ -309,52 +309,52 @@ static mxArray *
allocate_grid(struct processed_grid *grid, const char *func)
/* ---------------------------------------------------------------------- */
{
size_t nflds, nhf;
const char *fields[] = { "nodes", "faces", "cells",
"type", "cartDims", "griddim" };
size_t nflds, nhf;
const char *fields[] = { "nodes", "faces", "cells",
"type", "cartDims", "griddim" };
mxArray *G, *nodes, *faces, *cells;
mxArray *type, *typestr, *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);
nflds = sizeof(fields) / sizeof(fields[0]);
nhf = count_halffaces(grid->number_of_faces, grid->face_neighbors);
G = mxCreateStructMatrix(1, 1, nflds, fields);
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);
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);
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 );
mxSetField(G, 0, "faces" , faces );
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); }
mxSetField(G, 0, "nodes" , nodes );
mxSetField(G, 0, "faces" , faces );
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); }
G = NULL;
}
G = NULL;
}
return G;
return G;
}
@@ -363,16 +363,16 @@ static void
fill_grid(mxArray *G, struct processed_grid *grid)
/* ---------------------------------------------------------------------- */
{
double *pr;
double *pr;
pr = mxGetPr(mxGetField(G, 0, "cartDims"));
pr[0] = grid->dimensions[0];
pr[1] = grid->dimensions[1];
pr[2] = grid->dimensions[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);
fill_nodes(mxGetField(G, 0, "nodes"), grid);
fill_faces(mxGetField(G, 0, "faces"), grid);
fill_cells(mxGetField(G, 0, "cells"), grid);
}
@@ -381,23 +381,23 @@ static int
args_ok(int nlhs, int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
{
int ok;
int ok;
ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2));
ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2));
ok = ok && !mxIsEmpty(prhs[0]);
ok = ok && mxIsStruct(prhs[0]);
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);
ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0);
ok = ok && (mxGetFieldNumber(prhs[0], "cartDims") >= 0);
ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0);
ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0);
ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0);
if (ok && (nrhs == 2)) {
ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1);
}
if (ok && (nrhs == 2)) {
ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1);
}
return ok;
return ok;
}
@@ -406,56 +406,60 @@ static double
define_tolerance(int nrhs, const mxArray *prhs[])
/* ---------------------------------------------------------------------- */
{
double tol;
double tol;
tol = 0.0;
tol = 0.0;
if (nrhs == 2) {
tol = mxGetScalar(prhs[1]);
}
if (nrhs == 2) {
tol = mxGetScalar(prhs[1]);
}
return tol;
return tol;
}
/* 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;
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);
if (args_ok(nlhs, nrhs, prhs)) {
mx_init_grdecl(&grdecl, prhs[0]);
tolerance = define_tolerance(nrhs, prhs);
process_grdecl(&grdecl, tolerance, &g);
process_grdecl(&grdecl, tolerance, &g);
plhs[0] = allocate_grid(&g, mexFunctionName());
plhs[0] = allocate_grid(&g, mexFunctionName());
if (plhs[0] != NULL) {
fill_grid(plhs[0], &g);
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 {
/* Failed to create grid structure. Return empty. */
plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL);
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', and 'ACTNUM'",
mexFunctionName(), mexFunctionName());
mexErrMsgTxt(errmsg);
}
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', and 'ACTNUM'",
mexFunctionName(), mexFunctionName());
mexErrMsgTxt(errmsg);
}
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

433
uniquepoints.c
View File

@@ -20,7 +20,7 @@
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
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,
@@ -50,309 +50,306 @@
/*-----------------------------------------------------------------
Compare function passed to qsort
*/
Compare function passed to qsortx */
static int compare(const void *a, const void *b)
{
if (*(double*)a < *(double*) b) return -1;
else return 1;
if (*(double*)a < *(double*) b) return -1;
else return 1;
}
/*-----------------------------------------------------------------
Creat sorted list of z-values in zcorn with actnum==1
*/
static int createSortedList(double *list, int n, int m,
const double *z[], const int *a[])
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"); */
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;
qsort(list, ptr-list, sizeof(double), compare);
return ptr-list;
}
/*-----------------------------------------------------------------
Remove points less than <tolerance> apart in <list> of increasing
doubles.
*/
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));
int i;
int pos;
double val;
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;
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.
/*
Preserve outer z-boundary.
This operation is a no-op in the case
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];
list[pos-2] + tolerance < list[n-1].
return pos;
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)
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. */
/* 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. */
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;
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){
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;
/* 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;
}
/* 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 */
*p++ = INT_MAX;/* Padding to ease processing of faults */
return 1;
return 1;
}
/*-----------------------------------------------------------------
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[])
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[])
{
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);
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);
}
/*-----------------------------------------------------------------
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[])
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[])
{
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);
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);
}
/*-----------------------------------------------------------------
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 = (pt[2]-coord[2])/(coord[5]-coord[2]);
if (isinf(a) || isnan(a)){
a = 0;
}
double a = (pt[2]-coord[2])/(coord[5]-coord[2]);
if (isinf(a) || isnan(a)){
a = 0;
}
#if 0
pt[0] = coord[0] + a*(coord[3]-coord[0]);
pt[1] = coord[1] + a*(coord[4]-coord[1]);
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];
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)
*/
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)
/* 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];
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];
/* ztab->data may need extra space temporarily due to simple boundary treatement */
int npillarpoints = 8*(nx+1)*(ny+1)*nz;
int npillars = (nx+1)*(ny+1);
/* 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);
double *zlist = malloc(npillarpoints*sizeof *zlist);
int *zptr = malloc((npillars+1)*sizeof *zptr);
int i,j,k;
int i,j,k;
int d1[3];
int len = 0;
double *zout = zlist;
int pos = 0;
double *coord = (double*)g->coord;
double *pt;
const double *z[4];
const int *a[4];
int *p;
int pix, cix;
int zix;
int d1[3];
int len = 0;
double *zout = zlist;
int pos = 0;
double *coord = (double*)g->coord;
double *pt;
const double *z[4];
const int *a[4];
int *p;
int pix, cix;
int zix;
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));
out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates));
zptr[pos++] = zout - zlist;
zptr[pos++] = zout - zlist;
pt = out->node_coordinates;
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);
/* 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){
len = createSortedList( zout, d1[2], 4, z, a);
len = uniquify (len, zout, tolerance);
/* 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);
/* Assign unique points */
for (k=0; k<len; ++k){
pt[2] = zout[k];
interpolate_pillar(coord, pt);
pt += 3;
}
len = createSortedList( zout, d1[2], 4, z, a);
len = uniquify (len, zout, tolerance);
/* Increment pointer to sparse table of unique zcorn values */
zout = zout + len;
zptr[pos++] = zout - zlist;
/* Assign unique points */
for (k=0; k<len; ++k){
pt[2] = zout[k];
interpolate_pillar(coord, pt);
pt += 3;
}
coord += 6;
/* 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];
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]);
/* 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){
/* 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;
}
/* pillar index */
pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2);
p += 2 + 2*g->dims[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);
free(zptr);
free(zlist);
return 1;
return 1;
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

38
uniquepoints.h
View File

@@ -13,34 +13,34 @@
//==========================================================================*/
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
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 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.
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/>.
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);
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 */