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

Move variable declarations to top of each function. Declare function signatures in the top of each .c file (to avoid warnings).

Browse Source 
Signed-off-by: Bård Skaflestad <Bard.Skaflestad@sintef.no>
This commit is contained in:
Jostein R. Natvig 2012-01-27 12:27:02 +00:00 committed by Bård Skaflestad
parent 1600f1b4ca
commit 51881748d5
5 changed files with 132 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
facetopology.c
View File

@ -69,7 +69,8 @@ static int *computeFaceTopology(int *a1,
int *faces) int *faces)
{ {
int mask[8] = {-1}; int mask[8] = {-1};
int k;
int *f;
/* Which pillar points should we use? */ /* Which pillar points should we use? */
if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; } 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[1] > b2[1]){ mask[2] = b2[1]; } else { mask[2] = a2[1]; }
@ -136,8 +137,7 @@ static int *computeFaceTopology(int *a1,
mask[3] = intersect[2]; mask[3] = intersect[2];
} }
} }
int k; f = faces;
int *f = faces;
for (k=7; k>=0; --k){ for (k=7; k>=0; --k){
if(mask[k] != -1){ if(mask[k] != -1){
*f++ = mask[k]; *f++ = mask[k];
@ -217,6 +217,7 @@ void findconnections(int n, int *pts[4],
int i,j=0; int i,j=0;
int intersect[4]= {-1}; int intersect[4]= {-1};
int *tmp;
/* for (i=0; i<2*n; work[i++]=-1); */ /* for (i=0; i<2*n; work[i++]=-1); */
for (i = 0; i<n-1; ++i){ for (i = 0; i<n-1; ++i){
@ -341,7 +342,7 @@ void findconnections(int n, int *pts[4],
/* 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] */
int *tmp; tmp = itop; itop = ibottom; ibottom = tmp; tmp = itop; itop = ibottom; ibottom = tmp;
/* Zero out the "new" itop */ /* Zero out the "new" itop */
for(j=0;j<n; ++j) itop[j]=-1; for(j=0;j<n; ++j) itop[j]=-1;

9
mxgrdecl.c
View File

@ -1,4 +1,4 @@
//=========================================================================== /*=========================================================================
// //
// File: mxgrdecl.c // File: mxgrdecl.c
// //
@ -10,7 +10,7 @@
// //
// $Revision$ // $Revision$
// //
//=========================================================================== //=======================================================================*/
/* /*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
@ -42,6 +42,7 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
#include "grdecl.h" #include "grdecl.h"
void mx_init_grdecl(struct grdecl *g, const mxArray *s);
/* Get COORD, ZCORN, ACTNUM and DIMS from mxArray. */ /* Get COORD, ZCORN, ACTNUM and DIMS from mxArray. */
/*-------------------------------------------------------*/ /*-------------------------------------------------------*/
@ -50,7 +51,7 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s)
int i,n; int i,n;
mxArray *field; mxArray *field;
int numel; int numel;
double *tmp;
mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL; mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL;
if (!mxIsStruct(s) if (!mxIsStruct(s)
@ -71,7 +72,7 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s)
mexErrMsgTxt("cartDims field must be 3 numbers"); mexErrMsgTxt("cartDims field must be 3 numbers");
} }
double *tmp = mxGetPr(cartdims); tmp = mxGetPr(cartdims);
n = 1; n = 1;
for (i=0; i<3; ++i){ for (i=0; i<3; ++i){
g->dims[i] = tmp[i]; g->dims[i] = tmp[i];

6
mxgrdecl.h
View File

@ -1,4 +1,4 @@
//=========================================================================== /*=========================================================================
// //
// File: mxgrdecl.h // File: mxgrdecl.h
// //
@ -10,7 +10,7 @@
// //
// $Revision$ // $Revision$
// //
//=========================================================================== //=======================================================================*/
/* /*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
@ -37,6 +37,6 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
void mx_init_grdecl (struct grdecl *g, const mxArray *s); void mx_init_grdecl (struct grdecl *g, const mxArray *s);
#endif // OPENRS_MXGRDECL_HEADER #endif /* OPENRS_MXGRDECL_HEADER */

140
preprocess.c
View File

@ -46,7 +46,15 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
#define min(i,j) ((i)<(j) ? (i) : (j)) #define min(i,j) ((i)<(j) ? (i) : (j))
#define max(i,j) ((i)>(j) ? (i) : (j)) #define max(i,j) ((i)>(j) ? (i) : (j))
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);
void process_horizontal_faces(int **intersections,
int *plist,
struct processed_grid *out);
/*----------------------------------------------------------------- /*-----------------------------------------------------------------
Given a vector <field> with k index running faster than i running Given a vector <field> with k index running faster than i running
@ -164,13 +172,17 @@ void process_vertical_faces(int direction,
{ {
int i,j; int i,j;
int *cornerpts[4]; int *cornerpts[4];
int d[3];
int f;
enum face_tag tag[] = { LEFT, BACK }; enum face_tag tag[] = { LEFT, BACK };
int *tmp;
int nx = out->dimensions[0]; int nx = out->dimensions[0];
int ny = out->dimensions[1]; int ny = out->dimensions[1];
int nz = out->dimensions[2]; int nz = out->dimensions[2];
int startface;
int num_intersections;
int *ptr;
int len;
for (j=0; j<ny+direction; ++j) { for (j=0; j<ny+direction; ++j) {
for (i=0; i<nx+1-direction; ++i){ for (i=0; i<nx+1-direction; ++i){
@ -180,7 +192,9 @@ void process_vertical_faces(int direction,
} }
/* Vectors of point numbers */ /* Vectors of point numbers */
int d[] = {2*nx, 2*ny, 2+2*nz}; 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){ if(direction==1){
@ -188,7 +202,7 @@ void process_vertical_faces(int direction,
/* ---> */ /* ---> */
/* 0 2 2 3 */ /* 0 2 2 3 */
/* rotate clockwise */ /* rotate clockwise */
int *tmp = cornerpts[1]; tmp = cornerpts[1];
cornerpts[1] = cornerpts[0]; cornerpts[1] = cornerpts[0];
cornerpts[0] = cornerpts[2]; cornerpts[0] = cornerpts[2];
cornerpts[2] = cornerpts[3]; cornerpts[2] = cornerpts[3];
@ -196,23 +210,23 @@ void process_vertical_faces(int direction,
} }
/* int startface = ftab->position; */ /* int startface = ftab->position; */
int startface = out->number_of_faces; startface = out->number_of_faces;
/* int num_intersections = *npoints - npillarpoints; */ /* int num_intersections = *npoints - npillarpoints; */
int num_intersections = out->number_of_nodes - num_intersections = out->number_of_nodes -
out->number_of_nodes_on_pillars; out->number_of_nodes_on_pillars;
findconnections(2*nz+2, cornerpts, findconnections(2*nz+2, cornerpts,
*intersections+4*num_intersections, *intersections+4*num_intersections,
work, out); work, out);
int *ptr = out->face_neighbors + 2*startface; ptr = out->face_neighbors + 2*startface;
int len = 2*out->number_of_faces - 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-1+direction, j-direction, ptr, len);
compute_cell_index(out->dimensions, i, j, ptr+1, len); compute_cell_index(out->dimensions, i, j, ptr+1, len);
/* Tag the new faces */ /* Tag the new faces */
int f = startface; f = startface;
for (; f < out->number_of_faces; ++f) { for (; f < out->number_of_faces; ++f) {
out->face_tag[f] = tag[direction]; out->face_tag[f] = tag[direction];
} }
@ -249,9 +263,15 @@ void process_horizontal_faces(int **intersections,
int *cell = out->local_cell_index; int *cell = out->local_cell_index;
int cellno = 0; int cellno = 0;
int *f, *n, *c[4];
int prevcell, thiscell;
int idx;
/* dimensions of plist */ /* dimensions of plist */
int d[] = {2*nx, 2*ny, 2+2*nz}; int d[3];
d[0] = 2*nx;
d[1] = 2*ny;
d[2] = 2+2*nz;
for(j=0; j<ny; ++j) { for(j=0; j<ny; ++j) {
@ -264,15 +284,14 @@ void process_horizontal_faces(int **intersections,
} }
int *f = out->face_nodes + out->face_ptr[out->number_of_faces]; f = out->face_nodes + out->face_ptr[out->number_of_faces];
int *n = out->face_neighbors + 2*out->number_of_faces; n = out->face_neighbors + 2*out->number_of_faces;
/* Vectors of point numbers */ /* Vectors of point numbers */
int *c[4];
igetvectors(d, 2*i+1, 2*j+1, plist, c); igetvectors(d, 2*i+1, 2*j+1, plist, c);
int prevcell = -1; prevcell = -1;
for (k = 1; k<nz*2+1; ++k){ for (k = 1; k<nz*2+1; ++k){
@ -285,7 +304,7 @@ void process_horizontal_faces(int **intersections,
/* If the pinch is a cell: */ /* If the pinch is a cell: */
if (k%2){ 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; cell[idx] = -1;
} }
} }
@ -301,7 +320,7 @@ void process_horizontal_faces(int **intersections,
out->face_tag[ out->number_of_faces] = TOP; out->face_tag[ out->number_of_faces] = TOP;
out->face_ptr[++out->number_of_faces] = f - out->face_nodes; out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
int thiscell = linearindex(out->dimensions, i,j,(k-1)/2); thiscell = linearindex(out->dimensions, i,j,(k-1)/2);
*n++ = prevcell; *n++ = prevcell;
*n++ = prevcell = thiscell; *n++ = prevcell = thiscell;
@ -371,8 +390,9 @@ compute_intersection_coordinates(int *intersections,
{ {
int n = out->number_of_nodes; int n = out->number_of_nodes;
int np = out->number_of_nodes_on_pillars; 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. */ /* Make sure the space allocated for nodes match the number of node. */
void *p = realloc (out->node_coordinates, 3*n*sizeof(double)); void *p = realloc (out->node_coordinates, 3*n*sizeof(double));
if (p) { if (p) {
@ -384,9 +404,7 @@ compute_intersection_coordinates(int *intersections,
/* Append intersections */ /* Append intersections */
int k; pt = out->node_coordinates + 3*np;
double *pt = out->node_coordinates + 3*np;
int *itsct = intersections;
for (k=np; k<n; ++k){ for (k=np; k<n; ++k){
approximate_intersection_pt(itsct, out->node_coordinates, pt); approximate_intersection_pt(itsct, out->node_coordinates, pt);
@ -397,6 +415,9 @@ compute_intersection_coordinates(int *intersections,
} }
/*----------------------------------------------------------------- /*-----------------------------------------------------------------
Public interface Public interface
*/ */
@ -405,22 +426,43 @@ void process_grdecl(const struct grdecl *in,
struct processed_grid *out) struct processed_grid *out)
{ {
int i,j,k; 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; struct grdecl g;
int *actnum, *iptr;
double *zcorn;
int sign, error;
double z1, z2;
int c1, c2;
double *dptr;
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;
int *global_cell_index;
/* internal */
int *intersections = malloc(BIGNUM* sizeof(*intersections));
/* internal */
int *work = malloc(2* (2*nz+2)* sizeof(*work));
/* Allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */
int *plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist);
for(i=0; i<4*(nz+1); ++i) work[i] = -1;
g.dims[0] = nx; g.dims[0] = nx;
g.dims[1] = ny; g.dims[1] = ny;
g.dims[2] = nz; g.dims[2] = nz;
int *actnum = malloc (nc * sizeof *actnum); actnum = malloc (nc * sizeof *actnum);
double *zcorn = malloc (nc * 8 * sizeof *zcorn); zcorn = malloc (nc * 8 * sizeof *zcorn);
/* Permute actnum */ /* Permute actnum */
int *iptr = actnum; iptr = actnum;
for (j=0; j<ny; ++j){ for (j=0; j<ny; ++j){
for (i=0; i<nx; ++i){ for (i=0; i<nx; ++i){
for (k=0; k<nz; ++k){ for (k=0; k<nz; ++k){
@ -434,7 +476,7 @@ void process_grdecl(const struct grdecl *in,
/* HACK */ /* HACK */
/* Check that ZCORN is strictly nodecreaing along pillars. If */ /* Check that ZCORN is strictly nodecreaing along pillars. If */
/* not, check if -ZCORN is strictly nondecreasing. */ /* not, check if -ZCORN is strictly nondecreasing. */
int sign, error;
for (sign = 1; sign>-2; sign = sign - 2){ for (sign = 1; sign>-2; sign = sign - 2){
error = 0; error = 0;
@ -442,11 +484,11 @@ void process_grdecl(const struct grdecl *in,
for (j=0; j<2*ny; ++j){ for (j=0; j<2*ny; ++j){
for (i=0; i<2*nx; ++i){ for (i=0; i<2*nx; ++i){
for (k=0; k<2*nz-1; ++k){ for (k=0; k<2*nz-1; ++k){
double z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))]; z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))];
double z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))]; z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))];
int c1 = i/2 + nx*(j/2 + ny*k/2); c1 = i/2 + nx*(j/2 + ny*k/2);
int c2 = i/2 + nx*(j/2 + ny*(k+1)/2); c2 = i/2 + nx*(j/2 + ny*(k+1)/2);
if (in->actnum[c1] && in->actnum[c2] && (z2 < z1)){ if (in->actnum[c1] && in->actnum[c2] && (z2 < z1)){
fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n"); fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n");
@ -469,7 +511,7 @@ void process_grdecl(const struct grdecl *in,
} }
/* Permute zcorn */ /* Permute zcorn */
double *dptr = zcorn; dptr = zcorn;
for (j=0; j<2*ny; ++j){ for (j=0; j<2*ny; ++j){
for (i=0; i<2*nx; ++i){ for (i=0; i<2*nx; ++i){
for (k=0; k<2*nz; ++k){ for (k=0; k<2*nz; ++k){
@ -488,7 +530,7 @@ void process_grdecl(const struct grdecl *in,
/* Code below assumes k index runs fastests, ie. that dimensions of /* Code below assumes k index runs fastests, ie. that dimensions of
table is permuted to (dims[2], dims[0], dims[1]) */ table is permuted to (dims[2], dims[0], dims[1]) */
const int BIGNUM = 64;
out->m = BIGNUM/3; out->m = BIGNUM/3;
out->n = BIGNUM; out->n = BIGNUM;
@ -509,16 +551,6 @@ void process_grdecl(const struct grdecl *in,
out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index); out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index);
/* internal */
int *intersections = malloc(BIGNUM* sizeof(*intersections));
/* internal */
int *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 */
int *plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist);
/* Do actual work here:*/ /* Do actual work here:*/
@ -540,15 +572,15 @@ void process_grdecl(const struct grdecl *in,
/* Convert to local cell numbers in face_neighbors */ /* Convert to local cell numbers in face_neighbors */
int *ptr=out->face_neighbors;; iptr=out->face_neighbors;;
for (i=0; i<out->number_of_faces*2; ++i, ++ptr){ for (i=0; i<out->number_of_faces*2; ++i, ++iptr){
if (*ptr != -1){ if (*iptr != -1){
*ptr = out->local_cell_index[*ptr]; *iptr = out->local_cell_index[*iptr];
} }
} }
/* Invert global-to-local map */ /* Invert global-to-local map */
int *global_cell_index = malloc(out->number_of_cells * global_cell_index = malloc(out->number_of_cells *
sizeof (*global_cell_index)); sizeof (*global_cell_index));
for (i=0; i<nx*ny*nz; ++i){ for (i=0; i<nx*ny*nz; ++i){
if(out->local_cell_index[i]!=-1){ if(out->local_cell_index[i]!=-1){

55
uniquepoints.c
View File

@ -85,12 +85,15 @@ static int createSortedList(double *list, int n, int m,
*/ */
static int uniquify(int n, double *list, double tolerance) static int uniquify(int n, double *list, double tolerance)
{ {
int i;
int pos;
double val;
assert (!(tolerance < 0.0)); assert (!(tolerance < 0.0));
if (n<1) return 0; if (n<1) return 0;
int i; pos = 0;
int pos = 0; val = list[pos++];/* Keep first value */
double val = list[pos++];/* Keep first value */
for (i=1; i<n; ++i){ for (i=1; i<n; ++i){
if (val + tolerance < list [i]){ if (val + tolerance < list [i]){
@ -249,9 +252,11 @@ int finduniquepoints(const struct grdecl *g,
struct processed_grid *out) struct processed_grid *out)
{ {
int nx = out->dimensions[0];
int ny = out->dimensions[1]; const int nx = out->dimensions[0];
int nz = out->dimensions[2]; 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 */ /* ztab->data may need extra space temporarily due to simple boundary treatement */
@ -263,31 +268,38 @@ int finduniquepoints(const struct grdecl *g,
int nc = g->dims[0]*g->dims[1]*g->dims[2];
out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates));
double *zlist = ztab->data; /* casting void* to double* */ double *zlist = ztab->data; /* casting void* to double* */
int *zptr = ztab->ptr; int *zptr = ztab->ptr;
int i,j,k; int i,j,k;
int d1[3] = {2*g->dims[0], 2*g->dims[1], 2*g->dims[2]}; int d1[3];
int len = 0; int len = 0;
double *zout = zlist; double *zout = zlist;
int pos = 0; 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];
out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates));
zptr[pos++] = zout - zlist; zptr[pos++] = zout - zlist;
double *coord = (double*)g->coord; pt = out->node_coordinates;
double *pt = out->node_coordinates;
/* Loop over pillars, find unique points on each pillar */ /* Loop over pillars, find unique points on each pillar */
for (j=0; j < g->dims[1]+1; ++j){ for (j=0; j < g->dims[1]+1; ++j){
for (i=0; i < g->dims[0]+1; ++i){ 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 */ /* Get positioned pointers for actnum and zcorn data */
igetvectors(g->dims, i, j, g->actnum, a); igetvectors(g->dims, i, j, g->actnum, a);
dgetvectors(d1, 2*i, 2*j, g->zcorn, z); dgetvectors(d1, 2*i, 2*j, g->zcorn, z);
@ -313,24 +325,23 @@ int finduniquepoints(const struct grdecl *g,
out->number_of_nodes = 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 */
int *p = plist; p = plist;
for (j=0; j < 2*g->dims[1]; ++j){ for (j=0; j < 2*g->dims[1]; ++j){
for (i=0; i < 2*g->dims[0]; ++i){ for (i=0; i < 2*g->dims[0]; ++i){
/* pillar index */ /* 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 */ /* 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 */ /* zcorn column position */
int zix = 2*g->dims[2]*(i+2*g->dims[0]*j); zix = 2*g->dims[2]*(i+2*g->dims[0]*j);
const int *a = g->actnum + cix;
const double *z = g->zcorn + zix;
if (!assignPointNumbers(zptr[pix], zptr[pix+1], zlist, 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"); fprintf(stderr, "Something went wrong in assignPointNumbers");
return 0; return 0;
} }