From 5d4e216cecb0ee1ecd19f829a92f2168fc8afdf1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 26 Jun 2012 11:44:30 +0200 Subject: [PATCH 01/72] Anchor Doxygen comments to current line. Existing mark-up (/** ... */) would erroneously apply the documentation of one field to the one below. Using /**< ... */ avoids this problem. Found by reading the Doxygen manual more carefully... --- opm/core/linalg/sparse_sys.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/core/linalg/sparse_sys.h b/opm/core/linalg/sparse_sys.h index 1485c2f4..ce34bc51 100644 --- a/opm/core/linalg/sparse_sys.h +++ b/opm/core/linalg/sparse_sys.h @@ -37,13 +37,13 @@ extern "C" { */ struct CSRMatrix { - size_t m; /** Number of rows */ - size_t nnz; /** Number of structurally non-zero elements */ + size_t m; /**< Number of rows */ + size_t nnz; /**< Number of structurally non-zero elements */ - int *ia; /** Row pointers */ - int *ja; /** Column indices */ + int *ia; /**< Row pointers */ + int *ja; /**< Column indices */ - double *sa; /** Matrix elements */ + double *sa; /**< Matrix elements */ }; From a566bc834666f130bec2d7dd21fbe8bc18ec27e7 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Thu, 14 Oct 2010 05:37:00 +0000 Subject: [PATCH 02/72] Comment. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/facetopology.c b/facetopology.c index 897eba61..66eb08b6 100644 --- a/facetopology.c +++ b/facetopology.c @@ -59,7 +59,7 @@ along with OpenRS. If not, see . /*------------------------------------------------------*/ -/* Determine face geometry first, then compute intersections. */ +/* Determine face topology first, then compute intersection. */ /* All intersections that occur are present in the final face geometry.*/ static int *computeFaceTopology(int *a1, int *a2, @@ -169,11 +169,14 @@ static int *computeFaceTopology(int *a1, -/* a) If we assume that the index increase when z increase for - each pillar (but only separately), we can use only the point indices. +/* a) If we assume that the index increase when z increase for each + pillar (but only separately), we can use only the point indices, + since we only need to compare z-values on one pillar at a time. - b) We assume no intersections occur on the first and last lines. - This is convenient in the identification of (unique) intersections. + b) We assume input is preprocessed such that no intersections occur + on the first and last lines, for instance by padding the grid + with extra cells. This is convenient in the identification of + (unique) intersections. */ From fbe81573329f7f19d26e43d1d5cc9b32cab9b85e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 20 Dec 2010 10:39:20 +0000 Subject: [PATCH 03/72] Export the cell-face tags (i.e., cells.faces(:,2)) back to M in a manner consistent with the traditional semantics of the grid_structure. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/processgrid.c b/processgrid.c index f15ee344..1e43ce94 100644 --- a/processgrid.c +++ b/processgrid.c @@ -184,16 +184,26 @@ void fill_grid(mxArray **out, struct processed_grid *grid) } } - mxArray *cellfaces = mxCreateNumericMatrix(num_half_faces, 1, + mxArray *cellfaces = mxCreateNumericMatrix(num_half_faces, 2, mxINT32_CLASS, mxREAL); { int *iptr = mxGetData(cellfaces); + int c1, c2, cf_tag; for (i=0; inumber_of_faces; ++i) { - int c1 = grid->face_neighbors[2*i]; - int c2 = grid->face_neighbors[2*i+1]; - if(c1 != -1) iptr[counter[c1]++] = i+1; - if(c2 != -1) iptr[counter[c2]++] = i+1; + cf_tag = 2 * grid->face_tag[i] + 1; /* [1, 3, 5] */ + c1 = grid->face_neighbors[2*i+0]; + c2 = grid->face_neighbors[2*i+1]; + if(c1 != -1) { + iptr[counter[c1] + 0*num_half_faces] = i+1; + iptr[counter[c1] + 1*num_half_faces] = cf_tag + 1; /* out */ + counter[c1] += 1; + } + if(c2 != -1) { + iptr[counter[c2] + 0*num_half_faces] = i+1; + iptr[counter[c2] + 1*num_half_faces] = cf_tag + 0; /* in */ + counter[c2] += 1; + } } } mxSetField(cells, 0, "faces", cellfaces); From b23ad048a296baefd8df7c7d72942686a8e4a8ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 20 Dec 2010 23:29:10 +0000 Subject: [PATCH 04/72] Suggest an alternative implementation of the 'processgrid' function. This, disabled, implementation follows a two-stage allocate+fill algorithm that is slightly easier to control in memory-tight environments. We still have to maintain two copies of the grid, though. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 408 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 408 insertions(+) diff --git a/processgrid.c b/processgrid.c index 1e43ce94..0719d2aa 100644 --- a/processgrid.c +++ b/processgrid.c @@ -46,6 +46,412 @@ along with OpenRS. If not, see . #include "mxgrdecl.h" +#if 0 +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_nodes(size_t nnodes) +/* ---------------------------------------------------------------------- */ +{ + size_t nflds; + const char *fields[] = { "num", "coords" }; + + mxArray *nodes, *num, *coords; + + nflds = sizeof(fields) / sizeof(fields[0]); + + nodes = mxCreateStructMatrix(1, 1, nflds, fields); + + num = mxCreateDoubleScalar(nnodes); + coords = mxCreateDoubleMatrix(nnodes, 3, mxREAL); + + if ((nodes != NULL) && (num != NULL) && (coords != NULL)) { + mxSetField(nodes, 0, "num" , num); + mxSetField(nodes, 0, "coords", coords); + } else { + if (coords != NULL) { mxDestroyArray(coords); } + if (num != NULL) { mxDestroyArray(num); } + if (nodes != NULL) { mxDestroyArray(nodes); } + + nodes = NULL; + } + + return nodes; +} + + +/* ---------------------------------------------------------------------- */ +static void +fill_nodes(mxArray *nodes, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + size_t i, j, nnodes; + double *coords; + + nnodes = grid->number_of_nodes; + coords = mxGetPr(mxGetField(nodes, 0, "coords")); + + for (i = 0; i < nnodes; ++i) { + for (j = 0; j < 3; ++j) { + coords[i + j*nnodes] = grid->node_coordinates[3*i + j]; + } + } +} + + +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_faces(size_t nf, size_t nfacenodes) +/* ---------------------------------------------------------------------- */ +{ + size_t nflds; + const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" }; + + mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag; + + nflds = sizeof(fields) / sizeof(fields[0]); + + faces = mxCreateStructMatrix(1, 1, nflds, fields); + + num = mxCreateDoubleScalar (nf); + neighbors = mxCreateNumericMatrix(nf , 2, mxINT32_CLASS, mxREAL); + nodePos = mxCreateNumericMatrix(nf + 1 , 1, mxINT32_CLASS, mxREAL); + nodes = mxCreateNumericMatrix(nfacenodes, 1, mxINT32_CLASS, mxREAL); + tag = mxCreateNumericMatrix(nf , 1, mxINT32_CLASS, mxREAL); + + if ((faces != NULL) && (num != NULL) && (neighbors != NULL) && + (nodePos != NULL) && (nodes != NULL) && (tag != NULL)) { + mxSetField(faces, 0, "num" , num); + mxSetField(faces, 0, "neighbors", neighbors); + mxSetField(faces, 0, "nodePos" , nodePos); + mxSetField(faces, 0, "nodes" , nodes); + mxSetField(faces, 0, "tag" , tag); + } else { + if (tag != NULL) { mxDestroyArray(tag); } + if (nodes != NULL) { mxDestroyArray(nodes); } + if (nodePos != NULL) { mxDestroyArray(nodePos); } + if (neighbors != NULL) { mxDestroyArray(neighbors); } + if (num != NULL) { mxDestroyArray(num); } + if (faces != NULL) { mxDestroyArray(faces); } + + faces = NULL; + } + + return faces; +} + + +/* ---------------------------------------------------------------------- */ +static void +fill_faces(mxArray *faces, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + size_t i, f, nf, nfn; + + int *pi; + + nf = grid->number_of_faces; + + /* Fill faces.neighbors */ + pi = mxGetData(mxGetField(faces, 0, "neighbors")); + for (f = 0; f < nf; f++) { + /* Add one for one-based indexing in M */ + pi[f + 0*nf] = grid->face_neighbors[2*f + 0] + 1; + pi[f + 1*nf] = grid->face_neighbors[2*f + 1] + 1; + } + + /* Fill faces.nodePos */ + pi = mxGetData(mxGetField(faces, 0, "nodePos")); + for (i = 0; i <= nf; i++) { pi[i] = grid->face_ptr[i] + 1; } + + /* Fill faces.nodes */ + pi = mxGetData(mxGetField(faces, 0, "nodes")); + nfn = grid->face_ptr[nf]; /* Total number of face nodes */ + for (i = 0; i < nfn; i++) { pi[i] = grid->face_nodes[i] + 1; } + + /* Fill faces.tag */ + pi = mxGetData(mxGetField(faces, 0, "tag")); + for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[i] + 1; } +} + + +/* ---------------------------------------------------------------------- */ +static size_t +count_halffaces(size_t nf, const int *neighbors) +/* ---------------------------------------------------------------------- */ +{ + int c1, c2; + size_t nhf, f; + + for (f = nhf = 0; f < nf; f++) { + c1 = neighbors[2*f + 0]; + c2 = neighbors[2*f + 1]; + + nhf += c1 >= 0; + nhf += c2 >= 0; + } + + return nhf; +} + + +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_cells(size_t nc, size_t ncf) +/* ---------------------------------------------------------------------- */ +{ + size_t nflds; + const char *fields[] = { "num", "facePos", "faces", "indexMap" }; + + mxArray *cells, *num, *facePos, *faces, *indexMap; + + nflds = sizeof(fields) / sizeof(fields[0]); + + cells = mxCreateStructMatrix(1, 1, nflds, fields); + + num = mxCreateDoubleScalar (nc); + facePos = mxCreateNumericMatrix(nc + 1, 1, mxINT32_CLASS, mxREAL); + faces = mxCreateNumericMatrix(ncf , 2, mxINT32_CLASS, mxREAL); + indexMap = mxCreateNumericMatrix(nc , 1, mxINT32_CLASS, mxREAL); + + if ((cells != NULL) && (num != NULL) && (facePos != NULL) && + (faces != NULL) && (indexMap != NULL)) { + mxSetField(cells, 0, "num" , num ); + mxSetField(cells, 0, "facePos" , facePos ); + mxSetField(cells, 0, "faces" , faces ); + mxSetField(cells, 0, "indexMap", indexMap); + } else { + if (indexMap != NULL) { mxDestroyArray(indexMap); } + if (faces != NULL) { mxDestroyArray(faces); } + if (facePos != NULL) { mxDestroyArray(facePos); } + if (num != NULL) { mxDestroyArray(num); } + if (cells != NULL) { mxDestroyArray(cells); } + + cells = NULL; + } + + return cells; +} + + +/* ---------------------------------------------------------------------- */ +static void +fill_cells(mxArray *cells, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + size_t c, nc, f, nf, i; + + int c1, c2, cf_tag, nhf; + int *pi1, *pi2; + + nc = grid->number_of_cells; + nf = grid->number_of_faces; + + /* Simultaneously fill cells.facePos and cells.faces by transposing the + * neighbours mapping. */ + pi1 = mxGetData(mxGetField(cells, 0, "facePos")); + pi2 = mxGetData(mxGetField(cells, 0, "faces" )); + for (i = 0; i < nc + 1; i++) { pi1[i] = 0; } + + /* 1) Count connections (i.e., faces per cell). */ + for (f = 0; f < nf; f++) { + c1 = grid->face_neighbors[2*f + 0]; + c2 = grid->face_neighbors[2*f + 1]; + + if (c1 >= 0) { pi1[c1 + 1] += 1; } + if (c2 >= 0) { pi1[c2 + 1] += 1; } + } + + /* 2) Define start pointers (really, position *end* pointers at start). */ + for (c = 1; c <= nc; c++) { + pi1[0] += pi1[c]; + pi1[c] = pi1[0] - pi1[c]; + } + + /* 3) Fill connection structure whilst advancing end pointers. */ + nhf = pi1[0]; + pi1[0] = 0; + + mxAssert (((size_t) nhf) == mxGetM(mxGetField(cells, 0, "faces")), + "Number of half faces (SIZE(cells.faces,1)) incorrectly " + "determined earlier."); + + for (f = 0; f < nf; f++) { + cf_tag = 2*grid->face_tag[f] + 1; /* [1, 3, 5] */ + c1 = grid->face_neighbors[2*f + 0]; + c2 = grid->face_neighbors[2*f + 1]; + + if (c1 >= 0) { + pi2[ pi1[ c1 + 1 ] + 0*nhf ] = f + 1; + pi2[ pi1[ c1 + 1 ] + 1*nhf ] = cf_tag + 1; /* out */ + + pi1[ c1 + 1 ] += 1; + } + if (c2 >= 0) { + pi2[ pi1[ c2 + 1 ] + 0*nhf ] = f + 1; + pi2[ pi1[ c2 + 1 ] + 1*nhf ] = cf_tag + 0; /* in */ + + pi1[ c2 + 1 ] += 1; + } + } + + /* Finally, adjust pointer array for one-based indexing in M. */ + for (i = 0; i < nc + 1; i++) { pi1[i] += 1; } + + /* Fill cells.indexMap. Note that despite the name, 'local_cell_index' + * really *is* the (zero-based) indexMap of the 'processed_grid'. */ + pi1 = mxGetData(mxGetField(cells, 0, "indexMap")); + for (c = 0; c < nc; c++) { pi1[c] = grid->local_cell_index[c] + 1; } +} + + +/* ---------------------------------------------------------------------- */ +static mxArray * +allocate_grid(struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + size_t nflds, nhf; + const char *fields[] = { "nodes", "faces", "cells", "type", "cartDims" }; + + mxArray *G, *nodes, *faces, *cells, *type, *cartDims; + + nflds = sizeof(fields) / sizeof(fields[0]); + nhf = count_halffaces(grid->number_of_faces, grid->face_neighbors); + + G = mxCreateStructMatrix(1, 1, nflds, fields); + + nodes = allocate_nodes(grid->number_of_nodes); + faces = allocate_faces(grid->number_of_faces, + grid->face_ptr[ grid->number_of_faces ]); + cells = allocate_cells(grid->number_of_cells, nhf); + type = mxCreateString("processgrid"); + cartDims = mxCreateDoubleMatrix(1, 3, mxREAL); + + if ((G != NULL) && (nodes != NULL) && (faces != NULL) && + (cells != NULL) && (type != NULL) && (cartDims != NULL)) { + 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); + } else { + if (cartDims != NULL) { mxDestroyArray(cartDims); } + 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; + } + + return G; +} + + +/* ---------------------------------------------------------------------- */ +static void +fill_grid(mxArray *G, struct processed_grid *grid) +/* ---------------------------------------------------------------------- */ +{ + double *pr; + + pr = mxGetPr(mxGetField(G, 0, "cartDims")); + pr[0] = grid->dimensions[0]; + pr[1] = grid->dimensions[1]; + pr[2] = grid->dimensions[2]; + + fill_nodes(mxGetField(G, 0, "nodes"), grid); + fill_faces(mxGetField(G, 0, "faces"), grid); + fill_cells(mxGetField(G, 0, "cells"), grid); +} + + +/* ---------------------------------------------------------------------- */ +static int +args_ok(int nlhs, int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int ok; + + ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2)); + + ok = ok && !mxIsEmpty(prhs[0]); + ok = ok && mxIsStruct(prhs[0]); + + ok = ok && (mxGetFieldNumber(prhs[0], "cartDims") >= 0); + ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0); + ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0); + ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0); + + if (nrhs == 2) { + ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1); + } + + return ok; +} + + +/* ---------------------------------------------------------------------- */ +static double +define_tolerance(int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + double tol; + + tol = 0.0; + + if (nrhs == 2) { + tol = mxGetScalar(prhs[1]); + } + + 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; + + if (args_ok(nlhs, nrhs, prhs)) { + mx_init_grdecl(&grdecl, prhs[0]); + tolerance = define_tolerance(nrhs, prhs); + + process_grdecl(&grdecl, tolerance, &g); + + plhs[0] = allocate_grid(&g); + + if (plhs[0] != NULL) { + fill_grid(plhs[0], &g); + } else { + /* Failed to create grid structure. Return empty. */ + plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL); + } + + free_processed_grid(&g); + } else { + sprintf(errmsg, + "Calling sequence is\n\t" + "G = %s(grdecl)\t%%or\n\t" + "G = %s(grdecl, tolerance)\n" + "The 'grdecl' must be a valid structure with fields\n" + "\t'cartDims', 'COORD', 'ZCORN', and 'ACTNUM'", + mexFunctionName(), mexFunctionName()); + mexErrMsgTxt(errmsg); + } +} + +#else + void fill_grid(mxArray **out, struct processed_grid *grid) { const char *names[] = {"nodes", "faces", "cells", "cartDims", "type"}; @@ -247,6 +653,7 @@ void fill_grid(mxArray **out, struct processed_grid *grid) } + /* Gateway routine for Matlab mex function. */ /*-------------------------------------------------------*/ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) @@ -275,3 +682,4 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* Free whatever was allocated in initGrdecl. */ free_processed_grid(&out); } +#endif From d78aae52763b932ddf211e81ffbab881c483a9c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 21 Dec 2010 08:47:17 +0000 Subject: [PATCH 05/72] The '.type' is supposed to be a cell array of strings, not a mere string. Update accordingly. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit While here, infer the '.type' value from the function name rather than hard-coding 'processgrid'. Signed-off-by: Bård Skaflestad --- processgrid.c | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/processgrid.c b/processgrid.c index 0719d2aa..53c34169 100644 --- a/processgrid.c +++ b/processgrid.c @@ -306,13 +306,13 @@ fill_cells(mxArray *cells, struct processed_grid *grid) /* ---------------------------------------------------------------------- */ static mxArray * -allocate_grid(struct processed_grid *grid) +allocate_grid(struct processed_grid *grid, const char *func) /* ---------------------------------------------------------------------- */ { size_t nflds, nhf; const char *fields[] = { "nodes", "faces", "cells", "type", "cartDims" }; - mxArray *G, *nodes, *faces, *cells, *type, *cartDims; + mxArray *G, *nodes, *faces, *cells, *type, *typestr, *cartDims; nflds = sizeof(fields) / sizeof(fields[0]); nhf = count_halffaces(grid->number_of_faces, grid->face_neighbors); @@ -323,11 +323,15 @@ allocate_grid(struct processed_grid *grid) faces = allocate_faces(grid->number_of_faces, grid->face_ptr[ grid->number_of_faces ]); cells = allocate_cells(grid->number_of_cells, nhf); - type = mxCreateString("processgrid"); + type = mxCreateCellMatrix(1, 1); + typestr = mxCreateString(func); cartDims = mxCreateDoubleMatrix(1, 3, mxREAL); - if ((G != NULL) && (nodes != NULL) && (faces != NULL) && - (cells != NULL) && (type != NULL) && (cartDims != NULL)) { + if ((G != NULL) && (nodes != NULL) && (faces != NULL) && + (cells != NULL) && (type != NULL) && (typestr != NULL) && + (cartDims != NULL)) { + mxSetCell(type, 0, typestr); + mxSetField(G, 0, "nodes" , nodes ); mxSetField(G, 0, "faces" , faces ); mxSetField(G, 0, "cells" , cells ); @@ -335,6 +339,7 @@ allocate_grid(struct processed_grid *grid) mxSetField(G, 0, "cartDims", cartDims); } else { 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); } @@ -428,7 +433,7 @@ mexFunction(int nlhs, mxArray *plhs[], process_grdecl(&grdecl, tolerance, &g); - plhs[0] = allocate_grid(&g); + plhs[0] = allocate_grid(&g, mexFunctionName()); if (plhs[0] != NULL) { fill_grid(plhs[0], &g); From 714e0528229072a0708d3a7abfbe44c73516bccb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 21 Dec 2010 09:31:54 +0000 Subject: [PATCH 06/72] Preserve parameter validity status in (nrhs == 2) case. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processgrid.c b/processgrid.c index 53c34169..c03233ed 100644 --- a/processgrid.c +++ b/processgrid.c @@ -388,7 +388,7 @@ args_ok(int nlhs, int nrhs, const mxArray *prhs[]) ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0); ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0); - if (nrhs == 2) { + if (ok && (nrhs == 2)) { ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1); } From 691c3a3b7b46f1d34743fbe4517b08dcc70ebe17 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sun, 9 Jan 2011 15:15:43 +0000 Subject: [PATCH 07/72] By jrn's blessing, enable the alternative implementation of 'processgrid'. Also remove the existing implementation. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit While here, assert copyright for 2011. Signed-off-by: Bård Skaflestad --- processgrid.c | 239 +------------------------------------------------- 1 file changed, 2 insertions(+), 237 deletions(-) diff --git a/processgrid.c b/processgrid.c index c03233ed..276c3fa5 100644 --- a/processgrid.c +++ b/processgrid.c @@ -13,8 +13,8 @@ //===========================================================================*/ /* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. +Copyright 2009, 2010, 2011 SINTEF ICT, Applied Mathematics. +Copyright 2009, 2010, 2011 Statoil ASA. This file is part of The Open Reservoir Simulator Project (OpenRS). @@ -46,7 +46,6 @@ along with OpenRS. If not, see . #include "mxgrdecl.h" -#if 0 /* ---------------------------------------------------------------------- */ static mxArray * allocate_nodes(size_t nnodes) @@ -454,237 +453,3 @@ mexFunction(int nlhs, mxArray *plhs[], mexErrMsgTxt(errmsg); } } - -#else - -void fill_grid(mxArray **out, struct processed_grid *grid) -{ - const char *names[] = {"nodes", "faces", "cells", "cartDims", "type"}; - mxArray *G = mxCreateStructMatrix(1,1,sizeof names / sizeof names[0],names); - - int i,j; - double *ptr; - - - /* nodes */ - const char *n2[] = {"num", "coords"}; - mxArray *nodes = mxCreateStructMatrix(1,1,sizeof n2 / sizeof n2[0],n2); - mxSetField(nodes, 0, "num", mxCreateDoubleScalar(grid->number_of_nodes)); - - mxArray *nodecoords = mxCreateDoubleMatrix(grid->number_of_nodes, 3, mxREAL); - ptr = mxGetPr(nodecoords); - - for (j=0;j<3;++j) - { - for(i=0; inumber_of_nodes; ++i) - { - ptr[i+grid->number_of_nodes*j] = grid->node_coordinates[3*i+j]; - } - } - mxSetField(nodes, 0, "coords", nodecoords); - mxSetField(G, 0, "nodes", nodes); - - - /* faces */ - const char *n3[] = {"num", "neighbors", "nodes", "numNodes", "nodePos", "tag"}; - mxArray *faces = mxCreateStructMatrix(1,1,sizeof n3 / sizeof n3[0], n3); - - - mxSetField(faces, 0, "num", mxCreateDoubleScalar(grid->number_of_faces)); - - mxArray *faceneighbors = mxCreateNumericMatrix(grid->number_of_faces, 2, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(faceneighbors); - for(j=0; j<2; ++j) - { - for (i=0; inumber_of_faces; ++i) - { - int ix = grid->face_neighbors[2*i+j]; - if (ix == -1) - { - iptr[i+grid->number_of_faces*j] = 0; - } - else - { - iptr[i+grid->number_of_faces*j] = ix+1; - } - } - } - } - mxSetField(faces, 0, "neighbors", faceneighbors); - mxArray *numnodes = mxCreateNumericMatrix(grid->number_of_faces, 1, - mxINT32_CLASS, mxREAL); - mxArray *nodepos = mxCreateNumericMatrix(grid->number_of_faces+1, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr1 = mxGetData(numnodes); - int *iptr2 = mxGetData(nodepos); - iptr2[0] = 1; - for (i=0; inumber_of_faces; ++i) - { - iptr1[i] = grid->face_ptr[i+1] - grid->face_ptr[i]; - iptr2[i+1] = iptr2[i] + iptr1[i]; - } - } - mxSetField(faces, 0, "numNodes", numnodes); - mxSetField(faces, 0, "nodePos", nodepos); - - mxArray *tags = mxCreateNumericMatrix(grid->number_of_faces, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(tags); - for (i = 0; i < grid->number_of_faces; ++i) - { - iptr[i] = grid->face_tag[i] + 1; - } - } - mxSetField(faces, 0, "tag", tags); - - - const char *n4[] = {"num", "faces", "facePos", "indexMap"}; - mxArray *cells = mxCreateStructMatrix(1,1,sizeof n4 / sizeof n4[0], n4); - - mxSetField(cells, 0, "num", mxCreateDoubleScalar(grid->number_of_cells)); - - mxArray *map = mxCreateNumericMatrix(grid->number_of_cells, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(map); - for(i=0; inumber_of_cells; ++i) - { - iptr[i] = grid->local_cell_index[i]+1; - } - } - mxSetField(cells, 0, "indexMap", map); - - mxArray *facepos = mxCreateNumericMatrix(grid->number_of_cells+1, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(facepos); - for(i=0; inumber_of_cells+1; ++i) - { - iptr[i] = 0; - } - for (i=0; i<2*grid->number_of_faces; ++i) - { - int c=grid->face_neighbors[i]; - if(c != -1) - { - iptr[c+1]++; - } - } - iptr[0] = 1; - for(i=0; inumber_of_cells; ++i) - { - iptr[i+1] += iptr[i]; - } - } - mxSetField(cells, 0, "facePos", facepos); - - - - int *counter = calloc(grid->number_of_cells, sizeof(*counter)); - int num_half_faces = 0; - { - int *iptr = mxGetData(facepos); - for(i=0; inumber_of_cells; ++i) - { - counter[i] = num_half_faces; - num_half_faces += iptr[i+1]-iptr[i]; - } - } - - mxArray *cellfaces = mxCreateNumericMatrix(num_half_faces, 2, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(cellfaces); - int c1, c2, cf_tag; - for (i=0; inumber_of_faces; ++i) - { - cf_tag = 2 * grid->face_tag[i] + 1; /* [1, 3, 5] */ - c1 = grid->face_neighbors[2*i+0]; - c2 = grid->face_neighbors[2*i+1]; - if(c1 != -1) { - iptr[counter[c1] + 0*num_half_faces] = i+1; - iptr[counter[c1] + 1*num_half_faces] = cf_tag + 1; /* out */ - counter[c1] += 1; - } - if(c2 != -1) { - iptr[counter[c2] + 0*num_half_faces] = i+1; - iptr[counter[c2] + 1*num_half_faces] = cf_tag + 0; /* in */ - counter[c2] += 1; - } - } - } - mxSetField(cells, 0, "faces", cellfaces); - - mxSetField(G, 0, "cells", cells); - - int n = grid->face_ptr[grid->number_of_faces]; - - mxArray *facenodes = mxCreateNumericMatrix(n, 1, - mxINT32_CLASS, mxREAL); - { - int *iptr = mxGetData(facenodes); - for (i=0; iface_nodes[i]+1; - } - } - mxSetField(faces, 0, "nodes", facenodes); - - mxSetField(G, 0, "faces", faces); - - - free(counter); - - /* cartDims */ - mxArray *cartDims = mxCreateDoubleMatrix(1, 3, mxREAL); - ptr = mxGetPr(cartDims); - ptr[0] = grid->dimensions[0]; - ptr[1] = grid->dimensions[1]; - ptr[2] = grid->dimensions[2]; - - mxSetField(G, 0, "cartDims", cartDims); - - /* type */ - mxArray *type = mxCreateCellMatrix(1, 1); - mxSetCell(type, 0, mxCreateString("processgrid")); - - mxSetField(G, 0, "type", type); - - out[0] = G; -} - - - -/* Gateway routine for Matlab mex function. */ -/*-------------------------------------------------------*/ -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ - /* Set up data passed from Matlab */ - struct grdecl g; - struct processed_grid out; - double tolerance = 0.0; - - mx_init_grdecl(&g, prhs[0]); - if (nrhs == 2) - { - tolerance = mxGetScalar (prhs[1]); - } - - process_grdecl(&g, tolerance, &out); - - - if (plhs >0) - { - /* write to matlab struct */ - fill_grid(plhs, &out); - } - - - /* Free whatever was allocated in initGrdecl. */ - free_processed_grid(&out); -} -#endif From ba347a925781e79e7855aeb97a09072eddb3b8ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 8 Feb 2011 21:35:50 +0000 Subject: [PATCH 08/72] GCC 4.4 is perfectly fine. Don't bother using a non-default compiler that might not be available (or differently named) on another system. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- mexopts.sh | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 mexopts.sh diff --git a/mexopts.sh b/mexopts.sh deleted file mode 100644 index 24995736..00000000 --- a/mexopts.sh +++ /dev/null @@ -1,2 +0,0 @@ -. $MATLAB/bin/mexopts.sh -CC='gcc-4.1' From f1ebf2e2cfce4f666a8b3d22106e43952b22a355 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Fri, 8 Apr 2011 12:46:19 +0000 Subject: [PATCH 09/72] Fix stupid typo. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processgrid.c b/processgrid.c index 276c3fa5..a7fe553f 100644 --- a/processgrid.c +++ b/processgrid.c @@ -169,7 +169,7 @@ fill_faces(mxArray *faces, struct processed_grid *grid) /* Fill faces.tag */ pi = mxGetData(mxGetField(faces, 0, "tag")); - for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[i] + 1; } + for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[f] + 1; } } From 518904652891aff4fb39ec177355896e1f971625 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Wed, 18 May 2011 08:10:43 +0000 Subject: [PATCH 10/72] Separate build script and calling interface for processgrid. This enables simple Matlab postprocessing of output from C-code. This revision adds the 'griddim' field. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.m | 9 ++------- processgrid_mex.m | 42 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 7 deletions(-) create mode 100644 processgrid_mex.m diff --git a/processgrid.m b/processgrid.m index af3a1bcf..67ad1790 100644 --- a/processgrid.m +++ b/processgrid.m @@ -32,10 +32,5 @@ function varargout = processgrid(varargin) % $Date$ % $Revision$ -% Build MEX edition of same. -% -buildmex CFLAGS='$CFLAGS -Wall -fPIC' processgrid.c preprocess.c ... - uniquepoints.c facetopology.c sparsetable.c mxgrdecl.c - -% Call MEX edition. -[varargout{1:nargout}] = processgrid(varargin{:}); +[varargout{1:nargout}] = processgrid_mex(varargin{:}); +varargout{1}.griddim = 3; diff --git a/processgrid_mex.m b/processgrid_mex.m new file mode 100644 index 00000000..1becfb7a --- /dev/null +++ b/processgrid_mex.m @@ -0,0 +1,42 @@ +function varargout = processgrid_mex(varargin) +%Compute grid topology and geometry from pillar grid description. +% +% SYNOPSIS: +% G = processGRDECL(grdecl) +% +% PARAMETERS: +% grdecl - Raw pillar grid structure, as defined by function +% 'readGRDECL', with fields COORDS, ZCORN and, possibly, ACTNUM. +% +% RETURNS: +% G - Valid grid definition containing connectivity, cell +% geometry, face geometry and unique nodes. +% +% EXAMPLE: +% G = processgrid(readGRDECL('small.grdecl')); +% plotGrid(G); view(10,45); +% +% SEE ALSO: +% processGRDECL, readGRDECL, deactivateZeroPoro, removeCells. + +%{ +#COPYRIGHT# +%} + +% $Date$ +% $Revision$ + +% Copyright 2009 SINTEF ICT, Applied Mathematics. +% Mex gateway by Jostein R. Natvig, SINTEF ICT. + +% $Date$ +% $Revision$ + +% Build MEX edition of same. +% +buildmex CFLAGS='$CFLAGS -Wall -fPIC' processgrid.c preprocess.c ... + uniquepoints.c facetopology.c sparsetable.c mxgrdecl.c ... + -output processgrid_mex.mexa64 + +% Call MEX edition. +[varargout{1:nargout}] = processgrid_mex(varargin{:}); From f2b111f504b70d2b16ff49e12c7a0b3f0531b950 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 18 May 2011 16:16:58 +0000 Subject: [PATCH 11/72] Implement new, required field 'G.griddim'. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/processgrid.c b/processgrid.c index a7fe553f..e97dda5c 100644 --- a/processgrid.c +++ b/processgrid.c @@ -309,9 +309,11 @@ allocate_grid(struct processed_grid *grid, const char *func) /* ---------------------------------------------------------------------- */ { size_t nflds, nhf; - const char *fields[] = { "nodes", "faces", "cells", "type", "cartDims" }; + const char *fields[] = { "nodes", "faces", "cells", + "type", "cartDims", "griddim" }; - mxArray *G, *nodes, *faces, *cells, *type, *typestr, *cartDims; + 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); @@ -325,10 +327,11 @@ allocate_grid(struct processed_grid *grid, const char *func) 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)) { + 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 ); @@ -336,7 +339,9 @@ allocate_grid(struct processed_grid *grid, const char *func) 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); } From 1f0b08755c394bb04d8ed7afe34359a134d45abc Mon Sep 17 00:00:00 2001 From: "Halvor M. Nilsen" Date: Wed, 19 Oct 2011 07:35:37 +0000 Subject: [PATCH 12/72] bska fixed integer overflow in malloc call MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- sparsetable.c | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/sparsetable.c b/sparsetable.c index cb301469..9ef36bab 100644 --- a/sparsetable.c +++ b/sparsetable.c @@ -51,6 +51,7 @@ void free_sparse_table (sparse_table_t *tab) sparse_table_t *malloc_sparse_table(int m, int n, int datasz) { + size_t alloc_sz; sparse_table_t *tab = malloc(sizeof *tab); tab->m = m; tab->n = n; @@ -61,8 +62,9 @@ sparse_table_t *malloc_sparse_table(int m, int n, int datasz) return NULL; } - - if(!(tab->data = malloc(n * datasz))){ + alloc_sz = datasz; + alloc_sz *= n; + if(!(tab->data = malloc(alloc_sz))){ fprintf(stderr, "Could not allocate space for sparse data(%d)\n", n); free_sparse_table(tab); return NULL; From f11772b265466f13b25bfaeb0a1f2345a9c03649 Mon Sep 17 00:00:00 2001 From: "Halvor M. Nilsen" Date: Wed, 9 Nov 2011 14:26:38 +0000 Subject: [PATCH 13/72] corrected documentation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/processgrid.m b/processgrid.m index 67ad1790..87e6be58 100644 --- a/processgrid.m +++ b/processgrid.m @@ -2,11 +2,13 @@ function varargout = processgrid(varargin) %Compute grid topology and geometry from pillar grid description. % % SYNOPSIS: -% G = processGRDECL(grdecl) +% G = processgrid(grdecl) +% G = processgrid(grdecl,ztol) % % PARAMETERS: % grdecl - Raw pillar grid structure, as defined by function % 'readGRDECL', with fields COORDS, ZCORN and, possibly, ACTNUM. +% ztol - tolerance for unique points % % RETURNS: % G - Valid grid definition containing connectivity, cell From ca0e2600dbcd2b0ddca36a7ad24edb38597b37f1 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Mon, 16 Jan 2012 12:41:34 +0000 Subject: [PATCH 14/72] Do not check ZCORN monotonicity in inactive cells. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/preprocess.c b/preprocess.c index dbbddd65..3b08e731 100644 --- a/preprocess.c +++ b/preprocess.c @@ -445,7 +445,10 @@ void process_grdecl(const struct grdecl *in, double z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))]; double z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))]; - if (z2 < z1){ + int c1 = i/2 + nx*(j/2 + ny*k/2); + int c2 = i/2 + nx*(j/2 + ny*(k+1)/2); + + if (in->actnum[c1] && in->actnum[c2] && (z2 < z1)){ fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n"); error = 1; goto end; From d711a08b559b1403041e2db45efcdd1513dde68a Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Fri, 27 Jan 2012 10:03:23 +0000 Subject: [PATCH 15/72] To bring behavior of processgrid closer to that of processGRDECL, add post-processing to split grid in connected pieces. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.m | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/processgrid.m b/processgrid.m index 87e6be58..372cdbdc 100644 --- a/processgrid.m +++ b/processgrid.m @@ -1,4 +1,4 @@ -function varargout = processgrid(varargin) +function G = processgrid(varargin) %Compute grid topology and geometry from pillar grid description. % % SYNOPSIS: @@ -34,5 +34,33 @@ function varargout = processgrid(varargin) % $Date$ % $Revision$ -[varargout{1:nargout}] = processgrid_mex(varargin{:}); -varargout{1}.griddim = 3; + G = processgrid_mex(varargin{:}); + G.griddim = 3; + G = splitDisconnectedGrid(G, false); +end + +function G = splitDisconnectedGrid(G, verbose) + % Check if grid is connected + ix = all(G.faces.neighbors~=0, 2); + I = [G.faces.neighbors(ix,1);G.faces.neighbors(ix,2)]; + J = [G.faces.neighbors(ix,2);G.faces.neighbors(ix,1)]; + N = double(max(G.faces.neighbors(:))); + A = sparse(double(I),double(J),1,N,N)+speye(N); + clear ix I J + [a,b,c,d]=dmperm(A); %#ok + clear A b d + if numel(c) > 2, + dispif(verbose, '\nGrid has %d disconnected components\n', ... + numel(c)- 1); + % Partition grid into connected subgrids + for i = 1:numel(c) - 1, + g(i) = extractSubgrid(G, a(c(i):c(i+1)-1)); %#ok + sz(i) = g(i).cells.num; %#ok + g(i).cartDims = G.cartDims; %#ok + end + + % Return largest (in number of cells) grid first + [i,i] = sort(-sz); %#ok + G = g(i); + end +end From e8dd63b5dfaa9b7dc66d0d923decfe6eab2559c0 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Fri, 27 Jan 2012 10:06:11 +0000 Subject: [PATCH 16/72] Fix a couple of typo's. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/preprocess.c b/preprocess.c index 3b08e731..6e64cb2e 100644 --- a/preprocess.c +++ b/preprocess.c @@ -102,7 +102,7 @@ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len /*----------------------------------------------------------------- Ensure there's sufficient memory */ -int checkmemeory(int nz, struct processed_grid *out, int **intersections) +int checkmemory(int nz, struct processed_grid *out, int **intersections) { /* Ensure there is enough space */ @@ -174,7 +174,7 @@ void process_vertical_faces(int direction, for (j=0; jlocal_cell_index); out->local_cell_index = global_cell_index; - /* HACK: undo sign reversal in ZCORN prepprocessing */ + /* HACK: undo sign reversal in ZCORN preprocessing */ for (i=2; i<3*out->number_of_nodes; i = i+3) out->node_coordinates[i]=sign*out->node_coordinates[i]; From 1600f1b4caccb2a66518aa0f3c792c1a702dd93b Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Fri, 27 Jan 2012 12:25:31 +0000 Subject: [PATCH 17/72] Replace automatic build system with more warnings turned on- MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid_mex.m | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/processgrid_mex.m b/processgrid_mex.m index 1becfb7a..87fc65e0 100644 --- a/processgrid_mex.m +++ b/processgrid_mex.m @@ -32,11 +32,31 @@ function varargout = processgrid_mex(varargin) % $Date$ % $Revision$ -% Build MEX edition of same. -% -buildmex CFLAGS='$CFLAGS -Wall -fPIC' processgrid.c preprocess.c ... - uniquepoints.c facetopology.c sparsetable.c mxgrdecl.c ... - -output processgrid_mex.mexa64 + % Build MEX edition of same. + % + + v = version;v = v([1,3]); + + CFLAGS = {'CFLAGS="\$CFLAGS', '-g', '-Wall', '-Wextra', '-ansi', ... + '-pedantic', '-Wformat-nonliteral', '-Wcast-align', ... + '-Wpointer-arith', '-Wbad-function-cast', ... + '-Wmissing-prototypes ', '-Wstrict-prototypes', ... + '-Wmissing-declarations', '-Winline', '-Wundef', ... + '-Wnested-externs', '-Wcast-qual', '-Wshadow', ... + '-Wconversion', '-Wwrite-strings', '-Wno-conversion', ... + '-Wchar-subscripts', '-Wredundant-decls"'}; + + SRC = {'processgrid.c', 'preprocess.c', 'uniquepoints.c', ... + 'facetopology.c', 'sparsetable.c', 'mxgrdecl.c'}; + + INCLUDE = {}; + + OPTS = {'-output', ['processgrid_mex.', mexext], ... + '-largeArrayDims', ['-DMATLABVERSION=', v], '-g'}; + + buildmex(CFLAGS{:}, INCLUDE{:}, OPTS{:}, SRC{:}) + + % Call MEX edition. + [varargout{1:nargout}] = processgrid_mex(varargin{:}); +end -% Call MEX edition. -[varargout{1:nargout}] = processgrid_mex(varargin{:}); From 51881748d58c265d8e3c6ff1c03628552751bf91 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Fri, 27 Jan 2012 12:27:02 +0000 Subject: [PATCH 18/72] Move variable declarations to top of each function. Declare function signatures in the top of each .c file (to avoid warnings). MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 9 ++-- mxgrdecl.c | 9 ++-- mxgrdecl.h | 6 +-- preprocess.c | 140 ++++++++++++++++++++++++++++++------------------- uniquepoints.c | 55 +++++++++++-------- 5 files changed, 132 insertions(+), 87 deletions(-) diff --git a/facetopology.c b/facetopology.c index 66eb08b6..07f8d7d2 100644 --- a/facetopology.c +++ b/facetopology.c @@ -69,7 +69,8 @@ static int *computeFaceTopology(int *a1, 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]; } @@ -136,8 +137,7 @@ static int *computeFaceTopology(int *a1, mask[3] = intersect[2]; } } - int k; - int *f = faces; + f = faces; for (k=7; k>=0; --k){ if(mask[k] != -1){ *f++ = mask[k]; @@ -217,6 +217,7 @@ void findconnections(int n, int *pts[4], int i,j=0; int intersect[4]= {-1}; + int *tmp; /* for (i=0; i<2*n; work[i++]=-1); */ for (i = 0; i. #include "grdecl.h" +void mx_init_grdecl(struct grdecl *g, const mxArray *s); /* 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; mxArray *field; int numel; - + double *tmp; mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL; if (!mxIsStruct(s) @@ -71,7 +72,7 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s) mexErrMsgTxt("cartDims field must be 3 numbers"); } - double *tmp = mxGetPr(cartdims); + tmp = mxGetPr(cartdims); n = 1; for (i=0; i<3; ++i){ g->dims[i] = tmp[i]; diff --git a/mxgrdecl.h b/mxgrdecl.h index 54fa2267..e46c6973 100644 --- a/mxgrdecl.h +++ b/mxgrdecl.h @@ -1,4 +1,4 @@ -//=========================================================================== +/*========================================================================= // // File: mxgrdecl.h // @@ -10,7 +10,7 @@ // // $Revision$ // -//=========================================================================== +//=======================================================================*/ /* Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. @@ -37,6 +37,6 @@ along with OpenRS. If not, see . void mx_init_grdecl (struct grdecl *g, const mxArray *s); -#endif // OPENRS_MXGRDECL_HEADER +#endif /* OPENRS_MXGRDECL_HEADER */ diff --git a/preprocess.c b/preprocess.c index 6e64cb2e..cf667cb9 100644 --- a/preprocess.c +++ b/preprocess.c @@ -46,7 +46,15 @@ along with OpenRS. If not, see . #define min(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 with k index running faster than i running @@ -164,13 +172,17 @@ void process_vertical_faces(int direction, { 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 */ /* 0 2 2 3 */ /* rotate clockwise */ - int *tmp = cornerpts[1]; + tmp = cornerpts[1]; cornerpts[1] = cornerpts[0]; cornerpts[0] = cornerpts[2]; cornerpts[2] = cornerpts[3]; @@ -196,23 +210,23 @@ void process_vertical_faces(int direction, } /* int startface = ftab->position; */ - int startface = out->number_of_faces; + startface = out->number_of_faces; /* int num_intersections = *npoints - npillarpoints; */ - int num_intersections = out->number_of_nodes - + num_intersections = out->number_of_nodes - out->number_of_nodes_on_pillars; findconnections(2*nz+2, cornerpts, *intersections+4*num_intersections, work, out); - int *ptr = out->face_neighbors + 2*startface; - int 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); /* Tag the new faces */ - int f = startface; + f = startface; for (; f < out->number_of_faces; ++f) { out->face_tag[f] = tag[direction]; } @@ -249,9 +263,15 @@ void process_horizontal_faces(int **intersections, int *cell = out->local_cell_index; int cellno = 0; + int *f, *n, *c[4]; + int prevcell, thiscell; + int idx; /* 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; jface_nodes + out->face_ptr[out->number_of_faces]; - int *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 */ - int *c[4]; igetvectors(d, 2*i+1, 2*j+1, plist, c); - int prevcell = -1; + prevcell = -1; for (k = 1; kdimensions, i,j,(k-1)/2); + idx = linearindex(out->dimensions, i,j,(k-1)/2); cell[idx] = -1; } } @@ -301,7 +320,7 @@ void process_horizontal_faces(int **intersections, out->face_tag[ out->number_of_faces] = TOP; 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 = thiscell; @@ -371,8 +390,9 @@ compute_intersection_coordinates(int *intersections, { 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) { @@ -384,9 +404,7 @@ compute_intersection_coordinates(int *intersections, /* Append intersections */ - int k; - double *pt = out->node_coordinates + 3*np; - int *itsct = intersections; + pt = out->node_coordinates + 3*np; for (k=np; knode_coordinates, pt); @@ -397,6 +415,9 @@ compute_intersection_coordinates(int *intersections, } + + + /*----------------------------------------------------------------- Public interface */ @@ -405,22 +426,43 @@ void process_grdecl(const struct grdecl *in, struct processed_grid *out) { int i,j,k; - - int nx = in->dims[0]; - int ny = in->dims[1]; - int nz = in->dims[2]; - - int nc = nx*ny*nz; struct grdecl g; + int *actnum, *iptr; + 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[1] = ny; g.dims[2] = nz; - int *actnum = malloc (nc * sizeof *actnum); - double *zcorn = malloc (nc * 8 * sizeof *zcorn); + actnum = malloc (nc * sizeof *actnum); + zcorn = malloc (nc * 8 * sizeof *zcorn); /* Permute actnum */ - int *iptr = actnum; + iptr = actnum; for (j=0; j-2; sign = sign - 2){ error = 0; @@ -442,11 +484,11 @@ void process_grdecl(const struct grdecl *in, for (j=0; j<2*ny; ++j){ for (i=0; i<2*nx; ++i){ for (k=0; k<2*nz-1; ++k){ - double z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))]; - double z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))]; + z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))]; + z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))]; - int c1 = i/2 + nx*(j/2 + ny*k/2); - int c2 = i/2 + nx*(j/2 + ny*(k+1)/2); + c1 = i/2 + nx*(j/2 + ny*k/2); + c2 = i/2 + nx*(j/2 + ny*(k+1)/2); if (in->actnum[c1] && in->actnum[c2] && (z2 < z1)){ fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n"); @@ -469,7 +511,7 @@ void process_grdecl(const struct grdecl *in, } /* Permute zcorn */ - double *dptr = zcorn; + dptr = zcorn; for (j=0; j<2*ny; ++j){ for (i=0; i<2*nx; ++i){ 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 table is permuted to (dims[2], dims[0], dims[1]) */ - const int BIGNUM = 64; + out->m = BIGNUM/3; 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); - /* 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:*/ @@ -540,15 +572,15 @@ void process_grdecl(const struct grdecl *in, /* Convert to local cell numbers in face_neighbors */ - int *ptr=out->face_neighbors;; - for (i=0; inumber_of_faces*2; ++i, ++ptr){ - if (*ptr != -1){ - *ptr = out->local_cell_index[*ptr]; + iptr=out->face_neighbors;; + for (i=0; inumber_of_faces*2; ++i, ++iptr){ + if (*iptr != -1){ + *iptr = out->local_cell_index[*iptr]; } } /* 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)); for (i=0; ilocal_cell_index[i]!=-1){ diff --git a/uniquepoints.c b/uniquepoints.c index 11caecef..e00bb3f5 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -85,12 +85,15 @@ static int createSortedList(double *list, int n, int m, */ static int uniquify(int n, double *list, double tolerance) { + int i; + int pos; + double val; + assert (!(tolerance < 0.0)); if (n<1) return 0; - int i; - int pos = 0; - double val = list[pos++];/* Keep first value */ + pos = 0; + val = list[pos++];/* Keep first value */ for (i=1; idimensions[0]; - int ny = out->dimensions[1]; - int nz = out->dimensions[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 */ @@ -263,30 +268,37 @@ 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* */ int *zptr = ztab->ptr; 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; 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]; + + out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates)); zptr[pos++] = zout - zlist; - double *coord = (double*)g->coord; - double *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){ - - const int *a[4]; - const double *z[4]; /* Get positioned pointers for actnum and zcorn data */ igetvectors(g->dims, i, j, g->actnum, a); @@ -313,24 +325,23 @@ int finduniquepoints(const struct grdecl *g, out->number_of_nodes = zptr[pos-1]; /* 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 (i=0; i < 2*g->dims[0]; ++i){ /* pillar index */ - int pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2); + pix = (i+1)/2 + (g->dims[0]+1)*((j+1)/2); /* cell column position */ - int cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]); + cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]); /* zcorn column position */ - int zix = 2*g->dims[2]*(i+2*g->dims[0]*j); + 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, - 2*g->dims[2], z, a, p, tolerance)){ + 2*g->dims[2], + g->zcorn + zix, g->actnum + cix, + p, tolerance)){ fprintf(stderr, "Something went wrong in assignPointNumbers"); return 0; } From 982c04bcdb401b1043725929139301c1306c0d8e Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Fri, 27 Jan 2012 15:00:04 +0000 Subject: [PATCH 19/72] Order cells lexicographically. While here, fix misprint and adjust whitespace. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/preprocess.c b/preprocess.c index cf667cb9..2360e8e5 100644 --- a/preprocess.c +++ b/preprocess.c @@ -434,6 +434,7 @@ void process_grdecl(const struct grdecl *in, int c1, c2; double *dptr; const int BIGNUM = 64; + int cellnum; const int nx = in->dims[0]; const int ny = in->dims[1]; @@ -474,7 +475,7 @@ void process_grdecl(const struct grdecl *in, /* HACK */ - /* Check that ZCORN is strictly nodecreaing along pillars. If */ + /* Check that ZCORN is strictly nodecreasing along pillars. If */ /* not, check if -ZCORN is strictly nondecreasing. */ for (sign = 1; sign>-2; sign = sign - 2){ @@ -561,7 +562,7 @@ void process_grdecl(const struct grdecl *in, process_vertical_faces (0, &intersections, plist, work, out); process_vertical_faces (1, &intersections, plist, work, out); - process_horizontal_faces (&intersections, plist, out); + process_horizontal_faces ( &intersections, plist, out); free (plist); free (work); @@ -570,23 +571,31 @@ void process_grdecl(const struct grdecl *in, 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; ilocal_cell_index[i]!=-1){ + global_cell_index[cellnum] = i; + out->local_cell_index[i] = cellnum; + cellnum++; + } + } - /* Convert to local cell numbers in face_neighbors */ - iptr=out->face_neighbors;; + /* Remap out->face_neighbors */ + iptr=out->face_neighbors; for (i=0; inumber_of_faces*2; ++i, ++iptr){ if (*iptr != -1){ *iptr = out->local_cell_index[*iptr]; } } - /* Invert global-to-local map */ - global_cell_index = malloc(out->number_of_cells * - sizeof (*global_cell_index)); - for (i=0; ilocal_cell_index[i]!=-1){ - global_cell_index[out->local_cell_index[i]] = i; - } - } + free(out->local_cell_index); out->local_cell_index = global_cell_index; From a740a15d35cb8d06e2ea4f884f452e2f96b2381b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 27 Jan 2012 22:25:35 +0000 Subject: [PATCH 20/72] Fill faces.neighbors and nodes.coordinates along columns (i.e., unit stride in MATLAB). This means non-unit (but constant, 2 or 3) stride in C whence memory pre-fetching will typically be helpful. The cost is a (large) memory "rewind" operation at the end of each column. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/processgrid.c b/processgrid.c index e97dda5c..d1510f40 100644 --- a/processgrid.c +++ b/processgrid.c @@ -89,9 +89,9 @@ fill_nodes(mxArray *nodes, struct processed_grid *grid) nnodes = grid->number_of_nodes; coords = mxGetPr(mxGetField(nodes, 0, "coords")); - for (i = 0; i < nnodes; ++i) { - for (j = 0; j < 3; ++j) { - coords[i + j*nnodes] = 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]; } } } @@ -152,10 +152,11 @@ fill_faces(mxArray *faces, struct processed_grid *grid) /* Fill faces.neighbors */ pi = mxGetData(mxGetField(faces, 0, "neighbors")); - for (f = 0; f < nf; f++) { - /* Add one for one-based indexing in M */ - pi[f + 0*nf] = grid->face_neighbors[2*f + 0] + 1; - pi[f + 1*nf] = grid->face_neighbors[2*f + 1] + 1; + 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 */ From 88f59b1a9c750c0963c1521bc2222f97fc6bbbe3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 27 Jan 2012 22:27:39 +0000 Subject: [PATCH 21/72] Assert copyright for 2012. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/processgrid.c b/processgrid.c index d1510f40..03bc0942 100644 --- a/processgrid.c +++ b/processgrid.c @@ -13,8 +13,8 @@ //===========================================================================*/ /* -Copyright 2009, 2010, 2011 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010, 2011 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). From 24c4335a9a56671a82699185b5701fca7d197171 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Mon, 30 Jan 2012 09:51:35 +0000 Subject: [PATCH 22/72] Move some code from "process_grdecl" to helper functions, to shorten function- MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 148 +++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 131 insertions(+), 17 deletions(-) diff --git a/preprocess.c b/preprocess.c index 2360e8e5..362a41c0 100644 --- a/preprocess.c +++ b/preprocess.c @@ -415,7 +415,114 @@ compute_intersection_coordinates(int *intersections, } +static int* +copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out) +{ + int i,j,k; + int *ptr = out; + /* Permute actnum such that values of each vertical stack of cells + * are adjacent in memory, i.e., + out = [in(0,0,:), in(1,0,:),..., in(nx-1, ny-1,:)] + + in Matlab pseudo-code. + */ + for (j=0; jzcorn[i+2*nx*(j+2*ny*k)]; */ + /* } */ + /* } */ + /* } */ + +static int +get_zcorn_sign(int nx, int ny, int nz, const int *actnum, + const double *zcorn, int *error) +{ + /* Ensure that zcorn (i.e., depth) is strictly nondecreasing in + the k-direction. This is required by the processign algorithm. + + 1) if z(i,j,k) <= z(i,j,k+1) for all (i,j,k), return 1.0 + + 2) if -z(i,j,k) <=-z(i,j,k+1) for all (i,j,k), return -1.0 + + 3) if (1) and (2) fails, return -1.0, and set *error = 1. + + */ + int sign; + int i, j, k; + int c1, c2; + double z1, z2; + + for (sign = 1; sign>-2; sign = sign - 2) + { + *error = 0; + + for (j=0; j<2*ny; ++j){ + for (i=0; i<2*nx; ++i){ + for (k=0; k<2*nz-1; ++k){ + z1 = sign*zcorn[i+2*nx*(j+2*ny*(k))]; + z2 = sign*zcorn[i+2*nx*(j+2*ny*(k+1))]; + + c1 = i/2 + nx*(j/2 + ny*k/2); + c2 = i/2 + nx*(j/2 + ny*(k+1)/2); + + if (actnum[c1] && actnum[c2] && (z2 < z1)){ + fprintf(stderr, "\nZCORN should be strictly " + "nondecreasing along pillars!\n"); + *error = 1; + goto end; + } + } + } + } + + end: + if (!*error){ + break; + } + } + + if (*error){ + fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n" + "Grid definition may be broken\n"); + } + + return sign; +} /*----------------------------------------------------------------- @@ -425,22 +532,22 @@ void process_grdecl(const struct grdecl *in, double tolerance, struct processed_grid *out) { - int i,j,k; struct grdecl g; - int *actnum, *iptr; - double *zcorn; - int sign, error; - double z1, z2; - int c1, c2; - double *dptr; - const int BIGNUM = 64; - int cellnum; + int i; + int sign, error; + int cellnum; + + int *actnum, *iptr; + int *global_cell_index; + + 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; - int *global_cell_index; /* internal */ int *intersections = malloc(BIGNUM* sizeof(*intersections)); @@ -451,7 +558,7 @@ void process_grdecl(const struct grdecl *in, /* 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; + for(i=0; i<4*(nz+1); ++i) { work[i] = -1; } @@ -462,6 +569,7 @@ void process_grdecl(const struct grdecl *in, actnum = malloc (nc * sizeof *actnum); zcorn = malloc (nc * 8 * sizeof *zcorn); +#if 0 /* Permute actnum */ iptr = actnum; for (j=0; jactnum, actnum); + sign = get_zcorn_sign(nx, ny, nz, in->actnum, in->zcorn, &error); +#if 0 /* HACK */ /* Check that ZCORN is strictly nodecreasing along pillars. If */ /* not, check if -ZCORN is strictly nondecreasing. */ - for (sign = 1; sign>-2; sign = sign - 2){ + for (sign = 1; sign>-2; sign = sign - 2) + { error = 0; /* Ensure that zcorn is strictly nondecreasing in k-direction */ @@ -510,7 +623,9 @@ void process_grdecl(const struct grdecl *in, fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n" "Grid definition may be broken\n"); } - +#endif + +#if 0 /* Permute zcorn */ dptr = zcorn; for (j=0; j<2*ny; ++j){ @@ -520,11 +635,10 @@ void process_grdecl(const struct grdecl *in, } } } - - - - g.zcorn = zcorn; +#endif + g.zcorn = copy_and_permute_zcorn(nx, ny, nz, in->zcorn, sign, zcorn); + g.coord = in->coord; From fdc8af50b6d9b81afbe3e322a73c6f505d09a6fb Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Mon, 30 Jan 2012 10:17:51 +0000 Subject: [PATCH 23/72] 1) Remove disabled code 2) Restructure code in attempt to improve readability MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 178 ++++++++++++++++++++------------------------------- 1 file changed, 69 insertions(+), 109 deletions(-) diff --git a/preprocess.c b/preprocess.c index 362a41c0..9a1d8608 100644 --- a/preprocess.c +++ b/preprocess.c @@ -415,8 +415,10 @@ compute_intersection_coordinates(int *intersections, } +/* ------------------------------------------------------------------ */ static int* copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out) +/* ------------------------------------------------------------------ */ { int i,j,k; int *ptr = out; @@ -437,8 +439,11 @@ copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out) return out; } +/* ------------------------------------------------------------------ */ static double* -copy_and_permute_zcorn(int nx, int ny, int nz, const double *in, double sign, double *out) +copy_and_permute_zcorn(int nx, int ny, int nz, const double *in, + double sign, double *out) +/* ------------------------------------------------------------------ */ { int i,j,k; double *ptr = out; @@ -458,19 +463,12 @@ copy_and_permute_zcorn(int nx, int ny, int nz, const double *in, double sign, do } return out; } - /* /\* Permute zcorn *\/ */ - /* dptr = zcorn; */ - /* for (j=0; j<2*ny; ++j){ */ - /* for (i=0; i<2*nx; ++i){ */ - /* for (k=0; k<2*nz; ++k){ */ - /* *dptr++ = sign*in->zcorn[i+2*nx*(j+2*ny*k)]; */ - /* } */ - /* } */ - /* } */ +/* ------------------------------------------------------------------ */ static int get_zcorn_sign(int nx, int ny, int nz, const int *actnum, const double *zcorn, int *error) +/* ------------------------------------------------------------------ */ { /* Ensure that zcorn (i.e., depth) is strictly nondecreasing in the k-direction. This is required by the processign algorithm. @@ -549,103 +547,19 @@ void process_grdecl(const struct grdecl *in, const int nz = in->dims[2]; const int nc = nx*ny*nz; - /* 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; } + /* internal work arrays */ + int *work; + int *plist; + int *intersections; - - g.dims[0] = nx; - g.dims[1] = ny; - g.dims[2] = nz; - actnum = malloc (nc * sizeof *actnum); - zcorn = malloc (nc * 8 * sizeof *zcorn); -#if 0 - /* Permute actnum */ - iptr = actnum; - for (j=0; jactnum[i+nx*(j+ny*k)]; - } - } - } - g.actnum = actnum; -#endif - - g.actnum = copy_and_permute_actnum(nx, ny, nz, in->actnum, actnum); - sign = get_zcorn_sign(nx, ny, nz, in->actnum, in->zcorn, &error); - -#if 0 - /* HACK */ - /* Check that ZCORN is strictly nodecreasing along pillars. If */ - /* not, check if -ZCORN is strictly nondecreasing. */ - - for (sign = 1; sign>-2; sign = sign - 2) - { - error = 0; - - /* Ensure that zcorn is strictly nondecreasing in k-direction */ - for (j=0; j<2*ny; ++j){ - for (i=0; i<2*nx; ++i){ - for (k=0; k<2*nz-1; ++k){ - z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))]; - z2 = sign*in->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 (in->actnum[c1] && in->actnum[c2] && (z2 < z1)){ - fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n"); - error = 1; - goto end; - } - } - } - } - - end: - if (!error){ - break; - } - } - - if (error){ - fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n" - "Grid definition may be broken\n"); - } -#endif - -#if 0 - /* Permute zcorn */ - dptr = zcorn; - for (j=0; j<2*ny; ++j){ - for (i=0; i<2*nx; ++i){ - for (k=0; k<2*nz; ++k){ - *dptr++ = sign*in->zcorn[i+2*nx*(j+2*ny*k)]; - } - } - } - g.zcorn = zcorn; -#endif - g.zcorn = copy_and_permute_zcorn(nx, ny, nz, in->zcorn, sign, zcorn); - - g.coord = in->coord; - - - /* Code below assumes k index runs fastests, ie. that dimensions of - table is permuted to (dims[2], dims[0], dims[1]) */ - - + /* -----------------------------------------------------------------*/ + /* 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; @@ -655,9 +569,9 @@ void process_grdecl(const struct grdecl *in, out->face_tag = malloc( out->m * sizeof *out->face_tag); out->face_ptr[0] = 0; - out->dimensions[0] = g.dims[0]; - out->dimensions[1] = g.dims[1]; - out->dimensions[2] = g.dims[2]; + out->dimensions[0] = nx; + out->dimensions[1] = ny; + out->dimensions[2] = nz; out->number_of_faces = 0; out->number_of_nodes = 0; out->number_of_cells = 0; @@ -669,11 +583,49 @@ void process_grdecl(const struct grdecl *in, /* Do actual work here:*/ + /* -----------------------------------------------------------------*/ + /* For each pillar, compare zcorn values for adjacent cells to find + * the unique node z-coordinates specified by the input. While + * here, enumerate unique points and assign point numbers (in plist) + * for each cornerpoint cell. In other words, plist has 8 node + * numbers for each cornerpoint cell.*/ + + /* initialize grdecl structure "g" that will be processd by + * "finduniquepoints" */ + g.dims[0] = 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); + + zcorn = malloc (nc * 8 * sizeof *zcorn); + sign = get_zcorn_sign(nx, ny, nz, in->actnum, in->zcorn, &error); + g.zcorn = copy_and_permute_zcorn(nx, ny, nz, in->zcorn, sign, zcorn); + + g.coord = in->coord; + + + /* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */ + plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist); + finduniquepoints(&g, plist, tolerance, out); free (zcorn); free (actnum); + /* -----------------------------------------------------------------*/ + /* 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 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); @@ -681,10 +633,14 @@ void process_grdecl(const struct grdecl *in, free (plist); free (work); + /* -----------------------------------------------------------------*/ + /* (re)allocate space for and compute coordinates of nodes that + * arise from intesecting cells (faults) */ compute_intersection_coordinates(intersections, out); free (intersections); + /* -----------------------------------------------------------------*/ /* Enumerate compressed cells: -make array [0...#cells-1] of global cell numbers -make [0...nx*ny*nz-1] array of local cell numbers, @@ -713,9 +669,13 @@ void process_grdecl(const struct grdecl *in, free(out->local_cell_index); out->local_cell_index = global_cell_index; - /* HACK: undo sign reversal in ZCORN preprocessing */ - for (i=2; i<3*out->number_of_nodes; i = i+3) - out->node_coordinates[i]=sign*out->node_coordinates[i]; + /* if sign==-1 in ZCORN preprocessing, the sign of the z-coordinate + * need to change before we finish */ + if (sign == -1) + { + for (i=2; i<3*out->number_of_nodes; i = i+3) + out->node_coordinates[i]=sign*out->node_coordinates[i]; + } } From 0f2cbf66a2f0fecc2a7fdf90437b95eacf32333e Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Mon, 30 Jan 2012 11:41:50 +0000 Subject: [PATCH 24/72] Remove references to sparse_table_t, as this is only used to allocate and free two data vectors. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- uniquepoints.c | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/uniquepoints.c b/uniquepoints.c index e00bb3f5..bf589c26 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -41,7 +41,6 @@ #include -#include "sparsetable.h" #include "preprocess.h" #include "uniquepoints.h" @@ -262,16 +261,13 @@ int finduniquepoints(const struct grdecl *g, /* 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); - sparse_table_t *ztab = malloc_sparse_table(npillars, - npillarpoints, - sizeof(double)); + + double *zlist = malloc(npillarpoints*sizeof *zlist); + int *zptr = malloc((npillars+1)*sizeof *zptr); - double *zlist = ztab->data; /* casting void* to double* */ - int *zptr = ztab->ptr; - int i,j,k; int d1[3]; @@ -350,9 +346,8 @@ int finduniquepoints(const struct grdecl *g, } } - - free_sparse_table(ztab); - + free(zptr); + free(zlist); return 1; } From 03c752013d5982ccd309e4c6e3282980d79226b8 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Mon, 30 Jan 2012 11:51:09 +0000 Subject: [PATCH 25/72] Remove last references to sparse_table_t. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 1 - processgrid_mex.m | 2 +- sparsetable.c | 100 ---------------------------------------------- sparsetable.h | 53 ------------------------ 4 files changed, 1 insertion(+), 155 deletions(-) delete mode 100644 sparsetable.c delete mode 100644 sparsetable.h diff --git a/facetopology.c b/facetopology.c index 07f8d7d2..21e7e7e6 100644 --- a/facetopology.c +++ b/facetopology.c @@ -40,7 +40,6 @@ along with OpenRS. If not, see . #include #include "preprocess.h" -#include "sparsetable.h" #include "facetopology.h" /* No checking of input arguments in this code! */ diff --git a/processgrid_mex.m b/processgrid_mex.m index 87fc65e0..52ec3fb8 100644 --- a/processgrid_mex.m +++ b/processgrid_mex.m @@ -47,7 +47,7 @@ function varargout = processgrid_mex(varargin) '-Wchar-subscripts', '-Wredundant-decls"'}; SRC = {'processgrid.c', 'preprocess.c', 'uniquepoints.c', ... - 'facetopology.c', 'sparsetable.c', 'mxgrdecl.c'}; + 'facetopology.c', 'mxgrdecl.c'}; INCLUDE = {}; diff --git a/sparsetable.c b/sparsetable.c deleted file mode 100644 index 9ef36bab..00000000 --- a/sparsetable.c +++ /dev/null @@ -1,100 +0,0 @@ -/*=========================================================================== -// -// File: sparsetable.c -// -// Created: Fri Jun 19 08:48:04 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//==========================================================================*/ - -/* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. - -This file is part of The Open Reservoir Simulator Project (OpenRS). - -OpenRS is free software: you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -OpenRS is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with OpenRS. If not, see . -*/ - -#include -#include -#include -#include -#include - - -#include "sparsetable.h" - -void free_sparse_table (sparse_table_t *tab) -{ - if(tab->ptr) free(tab->ptr); - if(tab->data) free(tab->data); - - free(tab); -} - -sparse_table_t *malloc_sparse_table(int m, int n, int datasz) -{ - size_t alloc_sz; - sparse_table_t *tab = malloc(sizeof *tab); - tab->m = m; - tab->n = n; - tab->position = 0; - if (!(tab->ptr = malloc((m+1) * sizeof (*tab->ptr)))){ - fprintf(stderr, "Could not allocate space for sparse ptr\n"); - free_sparse_table(tab); - return NULL; - } - - alloc_sz = datasz; - alloc_sz *= n; - if(!(tab->data = malloc(alloc_sz))){ - fprintf(stderr, "Could not allocate space for sparse data(%d)\n", n); - free_sparse_table(tab); - return NULL; - } - - return tab; -} - -sparse_table_t *realloc_sparse_table(sparse_table_t *tab, int m, int n, int datasz) -{ - void *p = realloc(tab->ptr, (m+1) * sizeof (*tab->ptr)); - if (p){ - tab->ptr = p; - }else{ - fprintf(stderr, "Could not reallocate space for sparse ptr\n"); - free_sparse_table(tab); - return NULL; - } - - p = realloc(tab->data, n * datasz); - if(p){ - tab->data = p; - }else{ - fprintf(stderr, "Could not reallocate space for sparse data(%d)\n", n); - free_sparse_table(tab); - return NULL; - } - - tab->m = m; - tab->n = n; - - return tab; -} diff --git a/sparsetable.h b/sparsetable.h deleted file mode 100644 index f28df362..00000000 --- a/sparsetable.h +++ /dev/null @@ -1,53 +0,0 @@ -/*=========================================================================== -// -// File: sparsetable.h -// -// Created: Fri Jun 19 08:47:45 2009 -// -// Author: Jostein R. Natvig -// -// $Date$ -// -// $Revision$ -// -//===========================================================================*/ - -/* -Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010 Statoil ASA. - -This file is part of The Open Reservoir Simulator Project (OpenRS). - -OpenRS is free software: you can redistribute it and/or modify -it under the terms of the GNU General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -OpenRS is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU General Public License for more details. - -You should have received a copy of the GNU General Public License -along with OpenRS. If not, see . -*/ - -#ifndef OPENRS_SPARSETABLE_HEADER -#define OPENRS_SPARSETABLE_HEADER - -typedef struct{ - int m; /* number of rows */ - int *ptr; /* row pointer of size m+1 */ - int position; /* first position in ptr that is not filled. */ - - int n; /* size of data */ - void *data; /* sparse table data */ - - -} sparse_table_t; - -void free_sparse_table (sparse_table_t *tab); -sparse_table_t *malloc_sparse_table (int m, int n, int datasz); -sparse_table_t *realloc_sparse_table (sparse_table_t *tab, int m, int n, int datasz); - -#endif /* OPENRS_SPARSETABLE_HEADER */ From fb46be25c21d0f76e76eb791e26716738a3352c0 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Mon, 30 Jan 2012 11:51:28 +0000 Subject: [PATCH 26/72] Remove unused files. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- Makefile.am | 17 ------- libpreprocess.mk | 39 ---------------- newinterface.c | 114 ----------------------------------------------- newinterface.h | 56 ----------------------- 4 files changed, 226 deletions(-) delete mode 100644 Makefile.am delete mode 100644 libpreprocess.mk delete mode 100644 newinterface.c delete mode 100644 newinterface.h diff --git a/Makefile.am b/Makefile.am deleted file mode 100644 index a9285ce1..00000000 --- a/Makefile.am +++ /dev/null @@ -1,17 +0,0 @@ -# $Id: duneproject 5489 2009-03-25 11:19:24Z sander $ - -#SUBDIRS = test - -preprocessdir = $(includedir)/dune/grid/preprocess - -preprocess_HEADERS = facetopology.h preprocess.h sparsetable.h \ - uniquepoints.h - -lib_LTLIBRARIES = libpreprocess.la - -libpreprocess_la_SOURCES = facetopology.c uniquepoints.c \ - preprocess.c sparsetable.c - - -include $(top_srcdir)/am/global-rules - diff --git a/libpreprocess.mk b/libpreprocess.mk deleted file mode 100644 index 9664f7bc..00000000 --- a/libpreprocess.mk +++ /dev/null @@ -1,39 +0,0 @@ -PEDANTIC=-W -Wall -pedantic # -W -Wformat-nonliteral \ - #-Wcast-align -Wpointer-arith -Wbad-function-cast \ - #-Wmissing-prototypes -Wstrict-prototypes \ - #-Wmissing-declarations -Winline -Wundef -Wnested-externs\ - #-Wcast-qual -Wshadow -Wconversion -Wwrite-strings\ - #-Wno-conversion -Wchar-subscripts -Wredundant-decls\ - - - -GCC = gcc -fPIC -CC = $(GCC) -INC = -I. - -#rememeber export LD_LIBRARY_PATH=/home/jrn/devel/c/repo/preprocess -LIB = -L. -CFLAGS = -g $(INC) -Wall $(PEDANTIC) -SHELL = /bin/bash -DEFAULT_LIB_INSTALL_PATH = /home/jrn/devel/c/repo/opm-gridprocessing -all: tags depend libpreprocess - -OBJ = preprocess.o uniquepoints.o facetopology.o sparsetable.o newinterface.o ../reorder-utils/grid.o - -libpreprocess: libpreprocess.so.1.0.1 -libpreprocess.so.1.0.1: $(OBJ) - $(CC) -shared -Wl,-soname,libpreprocess.so.1,\ - -rpath,$(DEFAULT_LIB_INSTALL_PATH) -o $@ $(OBJ) -lc $(LIB) - ln -s libpreprocess.so.1.0.1 libpreprocess.so.1 - ln -s libpreprocess.so.1 libpreprocess.so -.PHONY: clean depend all - -clean:; rm -f *~ $(OBJ) libpreprocess.so.1.0.1 libpreprocess.so.1 \ - libpreprocess.so TAGS; makedepend - -tags : $(OBJ:.o=.c) $(wildcard *.h) - etags *.c *.h - -depend :; makedepend $(INC) -f makefile *.c 2>/dev/null - -# DO NOT DELETE THIS LINE -- make depend depends on it. diff --git a/newinterface.c b/newinterface.c deleted file mode 100644 index 6cbb81dd..00000000 --- a/newinterface.c +++ /dev/null @@ -1,114 +0,0 @@ -#include -#include "newinterface.h" - - - -static int *compute_cell_facepos(grid_t *g) -{ - int i,j,k; - int *facepos = malloc((g->number_of_cells + 1) * sizeof *facepos); - int *fcells = g->face_cells; - - for (i=0; inumber_of_cells; ++i) { - facepos [i] = 0; - } - - for (i=0; i<2*g->number_of_faces; ++i) { - if (*fcells != -1) { - (facepos[*fcells])++; - } - fcells++; - } - - /* cumsum */ - j=0; - for (i=0; inumber_of_cells; ++i) { - k = j + facepos[i]; - facepos[i] = j; - j = k; - } - facepos[i] = j; - - return facepos; -} - - -static int *compute_cell_faces(grid_t *g) -{ - int *cfaces = malloc(g->cell_facepos[g->number_of_cells] * sizeof *cfaces); - int *work = malloc(g->number_of_cells * sizeof *work); - int *fcells = g->face_cells; - int i,k,cell; - for(i=0; inumber_of_cells; ++i) { - work[i] = 0; - } - - for (i=0; inumber_of_faces; ++i) { - for (k=0;k<2; ++k) { - - if (*fcells != -1) { - cell = *fcells; - cfaces[g->cell_facepos[cell] + work[cell]] = i; - work[cell]++; - } - fcells++; - } - } - free(work); - - return cfaces; -} - -void preprocess (const struct grdecl *in, - double tol, - cornerpoint_grid_t *out) -{ - struct processed_grid pg; - process_grdecl(in, tol, &pg); - - /* - * General grid interface - */ - out->dimensions = 3; - - out->number_of_nodes = pg.number_of_nodes; - out->number_of_faces = pg.number_of_faces; - out->number_of_cells = pg.number_of_cells; - - out->node_coordinates = pg.node_coordinates; - - out->face_nodes = pg.face_nodes; - out->face_nodepos = pg.face_ptr; - out->face_cells = pg.face_neighbors; - - out->face_centroids = NULL; - out->face_normals = NULL; - out->face_areas = NULL; - - /* NB: compute_cell_facepos must be called before compute_cell_faces */ - out->cell_facepos = compute_cell_facepos((grid_t*) out); - out->cell_faces = compute_cell_faces ((grid_t*) out); - out->cell_centroids = NULL; - out->cell_volumes = NULL; - - - /* - * Cornerpoint grid interface - */ - out->cartdims[0] = pg.dimensions[0]; - out->cartdims[1] = pg.dimensions[1]; - out->cartdims[2] = pg.dimensions[2]; - - out->face_tag = pg.face_tag; - out->number_of_nodes_on_pillars = pg.number_of_nodes_on_pillars; - out->cartesian_cell_index = pg.local_cell_index; -} - -void free_cornerpoint_grid(cornerpoint_grid_t *g) -{ - free_grid((grid_t*) g); - - free(g->face_tag); - free(g->cartesian_cell_index); -} - diff --git a/newinterface.h b/newinterface.h deleted file mode 100644 index dbc2ea42..00000000 --- a/newinterface.h +++ /dev/null @@ -1,56 +0,0 @@ -/*=========================================================================== -// -// File: preprocess.h -// -// Created: Fri Jun 19 08:43:04 2009 -// -// Author: Jostein R. Natvig -// -// $Date: 2010-08-27 19:12:16 +0200 (Fri, 27 Aug 2010) $ -// -// $Revision: 930 $ -// -//==========================================================================*/ - -/* - Copyright 2010 SINTEF ICT, Applied Mathematics. -*/ - -#ifndef NEWINTERFACE_H -#define NEWINTERFACE_H -#include "preprocess.h" -#include "../reorder-utils/grid.h" - -#ifdef __cplusplus -extern "C" { -#endif - - typedef struct { - /* - * Common grid definitions - */ - GRID_TOPOLOGY - GRID_GEOMETRY - - /* - * Special cornerpoint definitions - */ - int cartdims[3]; - enum face_tag *face_tag; - int number_of_nodes_on_pillars; - int *cartesian_cell_index; - - } cornerpoint_grid_t; - - - void preprocess (const struct grdecl *in, - double tol, - cornerpoint_grid_t *out); - - void free_cornerpoint_grid(cornerpoint_grid_t *g); - -#ifdef __cplusplus -} -#endif - -#endif From 828141662e40b4d6fc5bc29b89987f2ead492204 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Mon, 30 Jan 2012 12:15:51 +0000 Subject: [PATCH 27/72] Add braces. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/facetopology.c b/facetopology.c index 21e7e7e6..b4fc8a75 100644 --- a/facetopology.c +++ b/facetopology.c @@ -220,7 +220,7 @@ void findconnections(int n, int *pts[4], /* for (i=0; i<2*n; work[i++]=-1); */ for (i = 0; iface_ptr[++out->number_of_faces] = f - out->face_nodes; @@ -332,8 +332,8 @@ void findconnections(int n, int *pts[4], } /* Update candidates for restart of j for in next i-iteration */ - if (b1[j] < a1[i+1]) k1 = j; - if (b2[j] < a2[i+1]) k2 = j; + if (b1[j] < a1[i+1]) { k1 = j; } + if (b2[j] < a2[i+1]) { k2 = j; } j = j+1; } @@ -344,7 +344,7 @@ void findconnections(int n, int *pts[4], tmp = itop; itop = ibottom; ibottom = tmp; /* Zero out the "new" itop */ - for(j=0;j Date: Mon, 30 Jan 2012 13:25:19 +0000 Subject: [PATCH 28/72] Parametrize macro, adjust braces. Remove commented-out warning message, that only confuse the reader. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 40 +++++++++++++++++++--------------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/facetopology.c b/facetopology.c index b4fc8a75..a819fa3a 100644 --- a/facetopology.c +++ b/facetopology.c @@ -190,8 +190,10 @@ static int faceintersection(int *a1, int *a2, int *b1, int *b2) } -#define meaningful_face !((a1[i] ==INT_MIN) && (b1[j] ==INT_MIN)) && \ - !((a1[i+1]==INT_MAX) && (b1[j+1]==INT_MAX)) +#define meaningful_face(i,j) \ + !((a1[i] ==INT_MIN) && (b1[j] ==INT_MIN)) && \ + !((a1[i+1]==INT_MAX) && (b1[j+1]==INT_MAX)) + /* work should be pointer to 2n ints initialised to zero . */ void findconnections(int n, int *pts[4], @@ -220,13 +222,18 @@ void findconnections(int n, int *pts[4], /* for (i=0; i<2*n; work[i++]=-1); */ for (i = 0; i Date: Tue, 31 Jan 2012 08:48:40 +0000 Subject: [PATCH 29/72] Re-indent to four-spaces. While at it, do whitespace-cleanup. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 508 ++++++++++++++++----------------- facetopology.h | 36 +-- grdecl.h | 18 +- mxgrdecl.c | 132 ++++----- mxgrdecl.h | 36 +-- preprocess.c | 752 ++++++++++++++++++++++++------------------------- processgrid.c | 526 +++++++++++++++++----------------- uniquepoints.c | 433 ++++++++++++++-------------- uniquepoints.h | 38 +-- 9 files changed, 1246 insertions(+), 1233 deletions(-) diff --git a/facetopology.c b/facetopology.c index a819fa3a..e9e1ec5f 100644 --- a/facetopology.c +++ b/facetopology.c @@ -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 . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ #include @@ -61,132 +61,132 @@ along with OpenRS. If not, see . /* 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)&&(a2b2))) 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; iface_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;jface_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. + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ -#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 */ diff --git a/grdecl.h b/grdecl.h index e9a636ae..f9045f32 100644 --- a/grdecl.h +++ b/grdecl.h @@ -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: */ diff --git a/mxgrdecl.c b/mxgrdecl.c index ab356ef7..fdb633e9 100644 --- a/mxgrdecl.c +++ b/mxgrdecl.c @@ -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 . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ #include @@ -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: */ diff --git a/mxgrdecl.h b/mxgrdecl.h index e46c6973..4621c202 100644 --- a/mxgrdecl.h +++ b/mxgrdecl.h @@ -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 . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ -#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: */ diff --git a/preprocess.c b/preprocess.c index 9a1d8608..9c8aa69f 100644 --- a/preprocess.c +++ b/preprocess.c @@ -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 . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ #include @@ -49,31 +49,30 @@ along with OpenRS. If not, see . 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 with k index running faster than i running faster than j, and Cartesian dimensions , find pointers to the (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. - */ + 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=dims[0] || j<0 || j >= dims[1]){ + for(k=0; km; - 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; jdimensions[0]; + int ny = out->dimensions[1]; + int nz = out->dimensions[2]; + int startface; + int num_intersections; + int *ptr; + int len; + for (j=0; j */ - /* 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, ..., #. 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; jface_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; kdimensions, 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; knode_coordinates, pt); - pt += 3; - itsct += 4; + for (k=np; knode_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; jdims[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; ilocal_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; ilocal_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; inumber_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; inumber_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: */ diff --git a/processgrid.c b/processgrid.c index 03bc0942..0cbfb814 100644 --- a/processgrid.c +++ b/processgrid.c @@ -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 . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ /* 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: */ diff --git a/uniquepoints.c b/uniquepoints.c index bf589c26..628c179e 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -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 apart in of increasing - doubles. -*/ + Remove points less than apart in 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 - 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 + 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 with k index running faster than i running - faster than j, and Cartesian dimensions , find pointers to the - (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. - */ -static void igetvectors(const int dims[3], int i, int j, - const int *field, const int *v[]) + Given a vector with k index running faster than i running + faster than j, and Cartesian dimensions , find pointers to the + (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of + field. */ +static void igetvectors(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 with k index running faster than i running - faster than j, and Cartesian dimensions , find pointers to the - (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of - field. - */ -static void dgetvectors(const int dims[3], int i, int j, - const double *field, const double *v[]) + Given a vector with k index running faster than i running + faster than j, and Cartesian dimensions , find pointers to the + (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of + field. */ +static void dgetvectors(const int dims[3], int i, int j, + const double *field, const double *v[]) { - - 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; knumber_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: */ - diff --git a/uniquepoints.h b/uniquepoints.h index ac3574bf..b9d698ac 100644 --- a/uniquepoints.h +++ b/uniquepoints.h @@ -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 . + You should have received a copy of the GNU General Public License + along with OpenRS. If not, see . */ -#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 */ From 08f1e2e3e67ec36a2d188d2f8c181fc1a9f84134 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Wed, 1 Feb 2012 11:48:04 +0000 Subject: [PATCH 30/72] Change code that computes new node coordinates on bilinear surfaces defined by pillar pairs. Old code computes *rough* approximation, new code computes exact value. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 51 ++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 46 insertions(+), 5 deletions(-) diff --git a/preprocess.c b/preprocess.c index 9c8aa69f..5b07097c 100644 --- a/preprocess.c +++ b/preprocess.c @@ -358,24 +358,65 @@ void process_horizontal_faces(int **intersections, pt holds coordinates to intersection between lines given by point numbers L[0]-L[1] and L[2]-L[3]. */ +#define OLD_INTERSECTION 0 static void approximate_intersection_pt(int *L, double *c, double *pt) { + double a; + double z0, z1, z2, z3; +#if OLD_INTERSECTION +#else + double b1, b2; + double x1, y1; + double x2, y2; + double z; +#endif - double z0 = c[3*L[0]+2]; - double z1 = c[3*L[1]+2]; - double z2 = c[3*L[2]+2]; - double z3 = c[3*L[3]+2]; + /* no intersection on pillars expected here! */ + assert(L[0]!=L[2]); + assert(L[1]!=L[3]); - double a = (z2-z0)/(z1-z0 - (z3-z2)); + z0 = c[3*L[0]+2]; + z1 = c[3*L[1]+2]; + z2 = c[3*L[2]+2]; + z3 = c[3*L[3]+2]; + + /* find parameter a where lines L0L1 and L2L3 have same + * z-coordinate */ + a = (z2-z0)/(z1-z0 - (z3-z2)); if (isinf(a) || isnan(a)){ a = 0; } +#if OLD_INTERSECTION 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; +#else + /* the corresponding z-coordinate is */ + z = z0* (1.0-a) + z1* a; + + + /* find point (x1, y1, z) on pillar 1 */ + b1 = (z2-z)/(z2-z0); + b2 = (z-z0)/(z2-z0); + x1 = c[3*L[0]+0]*b1 + c[3*L[2]+0]*b2; + y1 = c[3*L[0]+1]*b1 + c[3*L[2]+1]*b2; + + /* find point (x2, y2, z) on pillar 2 */ + b1 = (z-z3)/(z1-z3); + b2 = (z1-z)/(z1-z3); + x2 = c[3*L[1]+0]*b1 + c[3*L[3]+0]*b2; + y2 = c[3*L[1]+1]*b1 + c[3*L[3]+1]*b2; + + /* horizontal lines are by definition ON the bilinear surface + spanned by L0, L1, L2 and L3. find point (x, y, z) on + horizontal line between point (x1, y1, z) and (x2, y2, z).*/ + pt[0] = x1* (1.0-a) + x2* a; + pt[1] = y1* (1.0-a) + y2* a; + pt[2] = z; +#endif } /*----------------------------------------------------------------- From ea6c01dd95c9ceaf0fff5ad5eb9104d608ecd15a Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Wed, 1 Feb 2012 11:49:52 +0000 Subject: [PATCH 31/72] Remove references to old intersection implementation. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/preprocess.c b/preprocess.c index 5b07097c..734f7fb3 100644 --- a/preprocess.c +++ b/preprocess.c @@ -358,18 +358,14 @@ void process_horizontal_faces(int **intersections, pt holds coordinates to intersection between lines given by point numbers L[0]-L[1] and L[2]-L[3]. */ -#define OLD_INTERSECTION 0 static void approximate_intersection_pt(int *L, double *c, double *pt) { double a; double z0, z1, z2, z3; -#if OLD_INTERSECTION -#else double b1, b2; double x1, y1; double x2, y2; double z; -#endif /* no intersection on pillars expected here! */ assert(L[0]!=L[2]); @@ -387,12 +383,6 @@ static void approximate_intersection_pt(int *L, double *c, double *pt) a = 0; } - -#if OLD_INTERSECTION - 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; -#else /* the corresponding z-coordinate is */ z = z0* (1.0-a) + z1* a; @@ -409,14 +399,14 @@ static void approximate_intersection_pt(int *L, double *c, double *pt) b2 = (z1-z)/(z1-z3); x2 = c[3*L[1]+0]*b1 + c[3*L[3]+0]*b2; y2 = c[3*L[1]+1]*b1 + c[3*L[3]+1]*b2; - + /* horizontal lines are by definition ON the bilinear surface spanned by L0, L1, L2 and L3. find point (x, y, z) on horizontal line between point (x1, y1, z) and (x2, y2, z).*/ pt[0] = x1* (1.0-a) + x2* a; pt[1] = y1* (1.0-a) + y2* a; pt[2] = z; -#endif + } /*----------------------------------------------------------------- From e5ecded0a8c88f914121cea8b7bc2a6b04412d97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 20:24:02 +0000 Subject: [PATCH 32/72] Add convenience: Interpret ACTNUM==NULL as "all cells active". MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- mxgrdecl.c | 71 +++++++++++++++++++++++++++------------------------ preprocess.c | 33 ++++++++++++++++-------- processgrid.c | 4 ++- 3 files changed, 64 insertions(+), 44 deletions(-) diff --git a/mxgrdecl.c b/mxgrdecl.c index fdb633e9..6782f970 100644 --- a/mxgrdecl.c +++ b/mxgrdecl.c @@ -32,37 +32,33 @@ along with OpenRS. If not, see . */ -#include -#include -#include #include +#include +#include +#include +#include + #include - #include "grdecl.h" - - -void mx_init_grdecl(struct grdecl *g, const mxArray *s); +#include "mxgrdecl.h" /* Get COORD, ZCORN, ACTNUM and DIMS from mxArray. */ /*-------------------------------------------------------*/ void mx_init_grdecl(struct grdecl *g, const mxArray *s) { int i,n; - mxArray *field; - int numel; - double *tmp; + size_t numel; mxArray *cartdims=NULL, *actnum=NULL, *coord=NULL, *zcorn=NULL; if (!mxIsStruct(s) || !(cartdims = mxGetField(s, 0, "cartDims")) - || !(actnum = mxGetField(s, 0, "ACTNUM")) || !(coord = mxGetField(s, 0, "COORD")) - || !(zcorn = mxGetField(s, 0, "ZCORN")) + || !(zcorn = mxGetField(s, 0, "ZCORN")) ) { - char str[]="Input must be a single Matlab struct with fields\n" - "cartDims, ACTNUM, COORD and ZCORN\n"; + char str[]="Input must be a single MATLAB struct with fields\n" + "'cartDims', 'COORD' and 'ZCORN'. ACTNUM may be included.\n"; mexErrMsgTxt(str); } @@ -72,35 +68,44 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s) 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]; + if (mxIsDouble(cartdims)) { + double *tmp = mxGetPr(cartdims); + for (i = 0; i < 3; ++i) { + g->dims[i] = (int) tmp[i]; + } + } + else if (mxIsInt32(cartdims)) { + int *tmp = mxGetData(cartdims); + memcpy(g->dims, tmp, 3 * sizeof *g->dims); + } + + n = g->dims[0]; + for (i = 1; i < 3; i++) { n *= g->dims[ i ]; } + + + if ((actnum = mxGetField(s, 0, "ACTNUM")) != NULL) { + numel = mxGetNumberOfElements(actnum); + if ((! mxIsInt32(actnum)) || (numel != (size_t)(n))) { + mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32"); + } + g->actnum = mxGetData(actnum); + } + else { + g->actnum = NULL; } - 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); - - - - field = mxGetField(s, 0, "COORD"); numel = mxGetNumberOfElements(coord); - if (mxGetClassID(coord) != mxDOUBLE_CLASS || - numel != 6*(g->dims[0]+1)*(g->dims[1]+1)){ + if ((! mxIsDouble(coord)) || + numel != (size_t)(6*(g->dims[0]+1)*(g->dims[1]+1))) { mexErrMsgTxt("COORD field must have 6*(nx+1)*(ny+1) doubles."); } g->coord = mxGetPr(coord); numel = mxGetNumberOfElements(zcorn); - if (mxGetClassID(zcorn) != mxDOUBLE_CLASS || - numel != 8*g->dims[0]*g->dims[1]*g->dims[2]){ + if ((! mxIsDouble(zcorn)) || + numel != (size_t)(8*g->dims[0]*g->dims[1]*g->dims[2])) { mexErrMsgTxt("ZCORN field must have 8*nx*ny*nz doubles."); } g->zcorn = mxGetPr(zcorn); diff --git a/preprocess.c b/preprocess.c index 734f7fb3..14641ed6 100644 --- a/preprocess.c +++ b/preprocess.c @@ -451,20 +451,30 @@ copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out) { int i,j,k; int *ptr = out; + /* Permute actnum such that values of each vertical stack of cells * are adjacent in memory, i.e., - - out = [in(0,0,:), in(1,0,:),..., in(nx-1, ny-1,:)] - - in Matlab pseudo-code. - */ - for (j=0; j= 0); ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0); ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0); +#if 0 ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0); +#endif if (ok && (nrhs == 2)) { ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1); @@ -454,7 +456,7 @@ mexFunction(int nlhs, mxArray *plhs[], "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'", + "\t'cartDims', 'COORD', 'ZCORN'", mexFunctionName(), mexFunctionName()); mexErrMsgTxt(errmsg); } From ee9db9493c58fdaa4fe42578e4aae0b5f6f16cf0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 15 May 2012 21:12:25 +0000 Subject: [PATCH 33/72] Remove code that was disabled in -r967. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 3 --- 1 file changed, 3 deletions(-) diff --git a/processgrid.c b/processgrid.c index 47dfd376..cc66b6ec 100644 --- a/processgrid.c +++ b/processgrid.c @@ -391,9 +391,6 @@ args_ok(int nlhs, int nrhs, const mxArray *prhs[]) ok = ok && (mxGetFieldNumber(prhs[0], "cartDims") >= 0); ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0); ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0); -#if 0 - ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0); -#endif if (ok && (nrhs == 2)) { ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1); From de91035e3abb1e78d9208f0bdbe061cbc7a589fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 8 Jun 2012 11:40:43 +0000 Subject: [PATCH 34/72] Replace comparison function with one that honours qsort() requirements. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- uniquepoints.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/uniquepoints.c b/uniquepoints.c index 628c179e..b83847f4 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -53,8 +53,13 @@ Compare function passed to qsortx */ static int compare(const void *a, const void *b) { - if (*(double*)a < *(double*) b) return -1; - else return 1; + const double a0 = *(const double*) a; + const double b0 = *(const double*) b; + + /* { -1, a < b + * compare(a,b) = { 0, a = b + * { 1, a > b */ + return (a0 > b0) - (a0 < b0); } /*----------------------------------------------------------------- From 8884e4fc598da1b5c9438f197dcc418a9c3da9b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 8 Jun 2012 17:25:09 +0000 Subject: [PATCH 35/72] computeFaceTopology(): Initialise all elements of the mask array to -1. The ={-1} syntax just sets mask[0], leaving [1..7] at the default value (zero). MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/facetopology.c b/facetopology.c index e9e1ec5f..b5a10193 100644 --- a/facetopology.c +++ b/facetopology.c @@ -67,9 +67,12 @@ static int *computeFaceTopology(int *a1, int intersect[4], int *faces) { - int mask[8] = {-1}; + int mask[8]; int k; int *f; + + for (k = 0; k < 8; k++) { mask[k] = -1; } + /* Which pillar points should we use? */ if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; } if (a2[1] > b2[1]){ mask[2] = b2[1]; } else { mask[2] = a2[1]; } From 52c322047bef6665d18cb5870cfa0f1e9132fde4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 18 Jun 2012 14:23:55 +0000 Subject: [PATCH 36/72] Staticise several internal functions to reduce impact on libraries embedding this software. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 47 +++++++++++++++++++++++++++++------------------ 1 file changed, 29 insertions(+), 18 deletions(-) diff --git a/preprocess.c b/preprocess.c index 14641ed6..4c7324f1 100644 --- a/preprocess.c +++ b/preprocess.c @@ -46,15 +46,22 @@ #define min(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); +static void +compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len); + +static int +checkmemory(int nz, struct processed_grid *out, int **intersections); + +static void +process_vertical_faces(int direction, + int **intersections, + int *plist, int *work, + struct processed_grid *out); + +static void +process_horizontal_faces(int **intersections, + int *plist, + struct processed_grid *out); /*----------------------------------------------------------------- Given a vector with k index running faster than i running @@ -87,7 +94,8 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[]) other element in neighbors. */ -void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len) +static void +compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len) { int k; @@ -108,7 +116,8 @@ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len /*----------------------------------------------------------------- Ensure there's sufficient memory */ -int checkmemory(int nz, struct processed_grid *out, int **intersections) +static int +checkmemory(int nz, struct processed_grid *out, int **intersections) { /* Ensure there is enough space */ @@ -163,10 +172,11 @@ int checkmemory(int nz, struct processed_grid *out, int **intersections) 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) +static void +process_vertical_faces(int direction, + int **intersections, + int *plist, int *work, + struct processed_grid *out) { int i,j; int *cornerpts[4]; @@ -249,9 +259,10 @@ static int linearindex(const int dims[3], int i, int j, int k) ACTNUM==0) */ -void process_horizontal_faces(int **intersections, - int *plist, - struct processed_grid *out) +static void +process_horizontal_faces(int **intersections, + int *plist, + struct processed_grid *out) { int i,j,k; From 16169c57c525d00cba4e65beafbfc633f1b8e963 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 18 Jun 2012 15:39:46 +0000 Subject: [PATCH 37/72] Sort includes. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- uniquepoints.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/uniquepoints.c b/uniquepoints.c index b83847f4..bc73f8ab 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -33,12 +33,12 @@ */ #include -#include -#include -#include -#include #include +#include +#include #include +#include +#include #include "preprocess.h" From b21134fede7097c21579d7abe690b5bec7692f74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 18 Jun 2012 15:57:07 +0000 Subject: [PATCH 38/72] Remove dependency on C99's isinf() and isnan(). MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- uniquepoints.c | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/uniquepoints.c b/uniquepoints.c index bc73f8ab..ef6046ad 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -226,10 +226,18 @@ static void dgetvectors(const int dims[3], int i, int j, */ static void interpolate_pillar(const double *coord, double *pt) { - double a = (pt[2]-coord[2])/(coord[5]-coord[2]); - if (isinf(a) || isnan(a)){ + double a; + + if (fabs(coord[5] - coord[2]) > 0) { + + a = (pt[2] - coord[2]) / (coord[5] - coord[2]); + + } else { + a = 0; + } + #if 0 pt[0] = coord[0] + a*(coord[3]-coord[0]); pt[1] = coord[1] + a*(coord[4]-coord[1]); From 7dfb478e2e561d82e203f0861da4f21d171b23c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 18 Jun 2012 16:00:17 +0000 Subject: [PATCH 39/72] Remove dependency on C99's isinf() and isnan() functions. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/preprocess.c b/preprocess.c index 4c7324f1..1edf28e7 100644 --- a/preprocess.c +++ b/preprocess.c @@ -389,14 +389,18 @@ static void approximate_intersection_pt(int *L, double *c, double *pt) /* find parameter a where lines L0L1 and L2L3 have same * z-coordinate */ - a = (z2-z0)/(z1-z0 - (z3-z2)); - if (isinf(a) || isnan(a)){ + if (fabs((z1 - z0) - (z3 - z2)) > 0.0) { + + a = (z2 - z0) / ((z1 - z0) - (z3 - z2)); + + } else { + a = 0; + } /* the corresponding z-coordinate is */ - z = z0* (1.0-a) + z1* a; - + z = z0*(1.0 - a) + z1*a; /* find point (x1, y1, z) on pillar 1 */ From 140d4ade25ec8fb0e7a1ed4defaa4a2a2fb3f530 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 18 Jun 2012 16:03:11 +0000 Subject: [PATCH 40/72] Make a few concessions to readability by inserting some white-space around binary operators where appropriate. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/preprocess.c b/preprocess.c index 1edf28e7..573af101 100644 --- a/preprocess.c +++ b/preprocess.c @@ -244,7 +244,7 @@ process_vertical_faces(int direction, static int linearindex(const int dims[3], int i, int j, int k) { - return i+dims[0]*(j+dims[1]*k); + return i + dims[0]*(j + dims[1]*k); } @@ -379,13 +379,13 @@ static void approximate_intersection_pt(int *L, double *c, double *pt) double z; /* no intersection on pillars expected here! */ - assert(L[0]!=L[2]); - assert(L[1]!=L[3]); + assert (L[0] != L[2]); + assert (L[1] != L[3]); - z0 = c[3*L[0]+2]; - z1 = c[3*L[1]+2]; - z2 = c[3*L[2]+2]; - z3 = c[3*L[3]+2]; + z0 = c[3*L[0] + 2]; + z1 = c[3*L[1] + 2]; + z2 = c[3*L[2] + 2]; + z3 = c[3*L[3] + 2]; /* find parameter a where lines L0L1 and L2L3 have same * z-coordinate */ @@ -404,24 +404,23 @@ static void approximate_intersection_pt(int *L, double *c, double *pt) /* find point (x1, y1, z) on pillar 1 */ - b1 = (z2-z)/(z2-z0); - b2 = (z-z0)/(z2-z0); - x1 = c[3*L[0]+0]*b1 + c[3*L[2]+0]*b2; - y1 = c[3*L[0]+1]*b1 + c[3*L[2]+1]*b2; + b1 = (z2 - z) / (z2 - z0); + b2 = (z - z0) / (z2 - z0); + x1 = c[3*L[0] + 0]*b1 + c[3*L[2] + 0]*b2; + y1 = c[3*L[0] + 1]*b1 + c[3*L[2] + 1]*b2; /* find point (x2, y2, z) on pillar 2 */ - b1 = (z-z3)/(z1-z3); - b2 = (z1-z)/(z1-z3); - x2 = c[3*L[1]+0]*b1 + c[3*L[3]+0]*b2; - y2 = c[3*L[1]+1]*b1 + c[3*L[3]+1]*b2; + b1 = (z - z3) / (z1 - z3); + b2 = (z1 - z) / (z1 - z3); + x2 = c[3*L[1] + 0]*b1 + c[3*L[3] + 0]*b2; + y2 = c[3*L[1] + 1]*b1 + c[3*L[3] + 1]*b2; /* horizontal lines are by definition ON the bilinear surface spanned by L0, L1, L2 and L3. find point (x, y, z) on horizontal line between point (x1, y1, z) and (x2, y2, z).*/ - pt[0] = x1* (1.0-a) + x2* a; - pt[1] = y1* (1.0-a) + y2* a; + pt[0] = x1*(1.0 - a) + x2*a; + pt[1] = y1*(1.0 - a) + y2*a; pt[2] = z; - } /*----------------------------------------------------------------- From e9931d8c9c4eec8c0f37ffe5b57f3a101f494e51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 18 Jun 2012 16:19:41 +0000 Subject: [PATCH 41/72] Make 'coord' a pointer-to-const rather than a pointer-to-modifiable. We don't need to change the coord array. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- uniquepoints.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/uniquepoints.c b/uniquepoints.c index ef6046ad..c82817a8 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -283,7 +283,6 @@ int finduniquepoints(const struct grdecl *g, int len = 0; double *zout = zlist; int pos = 0; - double *coord = (double*)g->coord; double *pt; const double *z[4]; const int *a[4]; @@ -291,6 +290,8 @@ int finduniquepoints(const struct grdecl *g, int pix, cix; int zix; + const double *coord = g->coord; + d1[0] = 2*g->dims[0]; d1[1] = 2*g->dims[1]; d1[2] = 2*g->dims[2]; From 37351d2adb3942c3127323576bf9b815d935d1dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 18 Jun 2012 16:34:11 +0000 Subject: [PATCH 42/72] Initialise all "intersect" elements to -1 to avoid confusion surrounding ={-1} syntax. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/facetopology.c b/facetopology.c index b5a10193..a63f4206 100644 --- a/facetopology.c +++ b/facetopology.c @@ -220,10 +220,12 @@ void findconnections(int n, int *pts[4], int k2 = 0; int i,j=0; - int intersect[4]= {-1}; + int intersect[4]; int *tmp; /* for (i=0; i<2*n; work[i++]=-1); */ + for (i = 0; i < 4; i++) { intersect[i] = -1; } + for (i = 0; i Date: Mon, 18 Jun 2012 17:24:47 +0000 Subject: [PATCH 43/72] process_grdecl(): Prefer "size_t" when defining allocation sizes. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 53 ++++++++++++++++++++++++++-------------------------- 1 file changed, 26 insertions(+), 27 deletions(-) diff --git a/preprocess.c b/preprocess.c index 573af101..4700cc14 100644 --- a/preprocess.c +++ b/preprocess.c @@ -588,7 +588,7 @@ void process_grdecl(const struct grdecl *in, { struct grdecl g; - int i; + size_t i; int sign, error; int cellnum; @@ -597,11 +597,11 @@ void process_grdecl(const struct grdecl *in, 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 size_t BIGNUM = 64; + const int nx = in->dims[0]; + const int ny = in->dims[1]; + const int nz = in->dims[2]; + const size_t nc = ((size_t) nx) * ((size_t) ny) * ((size_t) nz); /* internal work arrays */ int *work; @@ -617,8 +617,8 @@ void process_grdecl(const struct grdecl *in, increased) 2) set Cartesian imensions */ - out->m = BIGNUM/3; - out->n = BIGNUM; + out->m = (int) (BIGNUM / 3); + out->n = (int) BIGNUM; out->face_neighbors = malloc( BIGNUM * sizeof *out->face_neighbors); out->face_nodes = malloc( out->n * sizeof *out->face_nodes); @@ -626,15 +626,15 @@ void process_grdecl(const struct grdecl *in, out->face_tag = malloc( out->m * sizeof *out->face_tag); out->face_ptr[0] = 0; - out->dimensions[0] = nx; - out->dimensions[1] = ny; - out->dimensions[2] = nz; + out->dimensions[0] = in->dims[0]; + out->dimensions[1] = in->dims[1]; + out->dimensions[2] = in->dims[2]; out->number_of_faces = 0; out->number_of_nodes = 0; out->number_of_cells = 0; out->node_coordinates = NULL; - out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index); + out->local_cell_index = malloc(nc * sizeof *out->local_cell_index); @@ -649,9 +649,9 @@ void process_grdecl(const struct grdecl *in, /* initialize grdecl structure "g" that will be processd by * "finduniquepoints" */ - g.dims[0] = nx; - g.dims[1] = ny; - g.dims[2] = nz; + g.dims[0] = in->dims[0]; + g.dims[1] = in->dims[1]; + g.dims[2] = in->dims[2]; actnum = malloc (nc * sizeof *actnum); g.actnum = copy_and_permute_actnum(nx, ny, nz, in->actnum, actnum); @@ -665,7 +665,7 @@ void process_grdecl(const struct grdecl *in, /* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) * padding */ - plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist); + plist = malloc(8 * (nc + ((size_t)nx)*((size_t)ny)) * sizeof *plist); finduniquepoints(&g, plist, tolerance, out); @@ -676,8 +676,8 @@ void process_grdecl(const struct grdecl *in, /* 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; } + work = malloc(2 * ((size_t) (2*nz + 2)) * sizeof *work); + for(i = 0; i < ((size_t)4) * (nz + 1); ++i) { work[i] = -1; } /* internal array to store intersections */ intersections = malloc(BIGNUM* sizeof(*intersections)); @@ -704,20 +704,19 @@ void process_grdecl(const struct grdecl *in, -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)); + global_cell_index = malloc(nc * sizeof *global_cell_index); cellnum = 0; - for (i=0; ilocal_cell_index[i]!=-1){ - global_cell_index[cellnum] = i; + for (i = 0; i < nc; ++i) { + if (out->local_cell_index[i] != -1) { + global_cell_index[cellnum] = (int) i; out->local_cell_index[i] = cellnum; cellnum++; } } /* Remap out->face_neighbors */ - iptr=out->face_neighbors; - for (i=0; inumber_of_faces*2; ++i, ++iptr){ + iptr = out->face_neighbors; + for (i = 0; i < (((size_t) 2) * out->number_of_faces; ++i, ++iptr) { if (*iptr != -1){ *iptr = out->local_cell_index[*iptr]; } @@ -731,8 +730,8 @@ void process_grdecl(const struct grdecl *in, * 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]; + for (i = 2; i < (((size_t) 3) * out->number_of_nodes); i += 3) + out->node_coordinates[i] *= sign; } } From 2864de1ddb88af96aac7b284d13ffda779943f12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 19 Jun 2012 07:42:49 +0000 Subject: [PATCH 44/72] Sort and group includes. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- processgrid.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/processgrid.c b/processgrid.c index cc66b6ec..9ef848ef 100644 --- a/processgrid.c +++ b/processgrid.c @@ -36,10 +36,11 @@ /* Mex gateway by Jostein R. Natvig, SINTEF ICT. */ -#include -#include -#include #include +#include +#include +#include + #include #include "preprocess.h" From ce073c03189c25768472416126a4c3960202d2e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 19 Jun 2012 07:51:15 +0000 Subject: [PATCH 45/72] Fix unmatched parentheses. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/preprocess.c b/preprocess.c index 4700cc14..426599d5 100644 --- a/preprocess.c +++ b/preprocess.c @@ -716,7 +716,7 @@ void process_grdecl(const struct grdecl *in, /* Remap out->face_neighbors */ iptr = out->face_neighbors; - for (i = 0; i < (((size_t) 2) * out->number_of_faces; ++i, ++iptr) { + for (i = 0; i < ((size_t) 2) * out->number_of_faces; ++i, ++iptr) { if (*iptr != -1){ *iptr = out->local_cell_index[*iptr]; } @@ -730,7 +730,7 @@ void process_grdecl(const struct grdecl *in, * z-coordinate need to change before we finish */ if (sign == -1) { - for (i = 2; i < (((size_t) 3) * out->number_of_nodes); i += 3) + for (i = 2; i < ((size_t) 3) * out->number_of_nodes; i += 3) out->node_coordinates[i] *= sign; } From fd3b5783fc02b9808d0855d925474f78c9132801 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 19 Jun 2012 07:51:26 +0000 Subject: [PATCH 46/72] Sort includes. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/facetopology.c b/facetopology.c index a63f4206..8617eb13 100644 --- a/facetopology.c +++ b/facetopology.c @@ -32,12 +32,12 @@ along with OpenRS. If not, see . */ -#include -#include -#include #include -#include #include +#include +#include +#include +#include #include "preprocess.h" #include "facetopology.h" From eaaad6424211b0950cee57203cd2d4ca3e5b2950 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 19 Jun 2012 12:58:51 +0000 Subject: [PATCH 47/72] ?getvectors(): Re-factor field offset calculations out to new helper, vector_positions(). MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- uniquepoints.c | 47 +++++++++++++++++++++++++++++++---------------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/uniquepoints.c b/uniquepoints.c index c82817a8..0dde8c20 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -177,6 +177,27 @@ static int assignPointNumbers(int begin, } +/* ---------------------------------------------------------------------- */ +static void +vector_positions(const int dims[3] , + const int i , + const int j , + size_t start[4]) +/* ---------------------------------------------------------------------- */ +{ + size_t im, ip, jm, jp; + + im = max(1, i ) - 1; + ip = min(dims[0], i+1) - 1; + jm = max(1, j ) - 1; + jp = min(dims[1], j+1) - 1; + + start[ 0 ] = dims[2] * (im + dims[0]*jm); + start[ 1 ] = dims[2] * (im + dims[0]*jp); + start[ 2 ] = dims[2] * (ip + dims[0]*jm); + start[ 3 ] = dims[2] * (ip + dims[0]*jp); +} + /*----------------------------------------------------------------- Given a vector with k index running faster than i running @@ -186,16 +207,13 @@ static int assignPointNumbers(int begin, static void igetvectors(const int dims[3], int i, int j, const int *field, const int *v[]) { + size_t p, start[4]; - int im = max(1, i ) - 1; - int ip = min(dims[0], i+1) - 1; - int jm = max(1, j ) - 1; - int jp = min(dims[1], j+1) - 1; + vector_positions(dims, i, j, start); - v[0] = field + dims[2]*(im + dims[0]* jm); - v[1] = field + dims[2]*(im + dims[0]* jp); - v[2] = field + dims[2]*(ip + dims[0]* jm); - v[3] = field + dims[2]*(ip + dims[0]* jp); + for (p = 0; p < 4; p++) { + v[p] = field + start[p]; + } } @@ -207,16 +225,13 @@ static void igetvectors(const int dims[3], int i, int j, static void dgetvectors(const int dims[3], int i, int j, const double *field, const double *v[]) { + size_t p, start[4]; - int im = max(1, i ) - 1; - int ip = min(dims[0], i+1) - 1; - int jm = max(1, j ) - 1; - int jp = min(dims[1], j+1) - 1; + vector_positions(dims, i, j, start); - v[0] = field + dims[2]*(im + dims[0]* jm); - v[1] = field + dims[2]*(im + dims[0]* jp); - v[2] = field + dims[2]*(ip + dims[0]* jm); - v[3] = field + dims[2]*(ip + dims[0]* jp); + for (p = 0; p < 4; p++) { + v[p] = field + start[p]; + } } From 1ee27a4484d3e423f0df7e55ed43cb2e31ed5fb2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 19 Jun 2012 13:22:05 +0000 Subject: [PATCH 48/72] Up-case "min" and "max" macros to provide visual cues that there's macro expansion afoot. While here, remove unused "overlap" macro. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- uniquepoints.c | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/uniquepoints.c b/uniquepoints.c index 0dde8c20..4e6779d9 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -44,10 +44,8 @@ #include "preprocess.h" #include "uniquepoints.h" -#define min(i,j) ((i)<(j) ? (i) : (j)) -#define max(i,j) ((i)>(j) ? (i) : (j)) -#define overlap(a1,a2,b1,b2) max(a1,b1) < min(a2,b2) - +#define MIN(i,j) (((i) < (j)) ? (i) : (j)) +#define MAX(i,j) (((i) > (j)) ? (i) : (j)) /*----------------------------------------------------------------- Compare function passed to qsortx */ @@ -187,10 +185,10 @@ vector_positions(const int dims[3] , { size_t im, ip, jm, jp; - im = max(1, i ) - 1; - ip = min(dims[0], i+1) - 1; - jm = max(1, j ) - 1; - jp = min(dims[1], j+1) - 1; + im = MAX(1, i ) - 1; + jm = MAX(1, j ) - 1; + ip = MIN(dims[0], i+1) - 1; + jp = MIN(dims[1], j+1) - 1; start[ 0 ] = dims[2] * (im + dims[0]*jm); start[ 1 ] = dims[2] * (im + dims[0]*jp); From 7d68712136d49faac9b6c1595098d5877279cd6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 09:43:16 +0000 Subject: [PATCH 49/72] Up-case all macro names. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit While here, insert white-space for readability. Signed-off-by: Bård Skaflestad --- facetopology.c | 36 +++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/facetopology.c b/facetopology.c index 8617eb13..d70d5025 100644 --- a/facetopology.c +++ b/facetopology.c @@ -43,8 +43,8 @@ #include "facetopology.h" /* No checking of input arguments in this code! */ -#define min(i,j) ((i)<(j) ? (i) : (j)) -#define max(i,j) ((i)>(j) ? (i) : (j)) +#define MIN(i,j) ((i)<(j) ? (i) : (j)) +#define MAX(i,j) ((i)>(j) ? (i) : (j)) #define DEBUG 1 @@ -182,18 +182,23 @@ static int *computeFaceTopology(int *a1, */ -#define lineintersection(a1,a2,b1,b2)(((a1>b1)&&(a2b2))) -static int faceintersection(int *a1, int *a2, int *b1, int *b2) +#define LINE_INTERSECTION(a1, a2, b1, b2) \ + (((a1 > b1) && (a2 < b2)) || \ + ((a1 < b1) && (a2 > b2))) + +static int +faceintersection(const int *a1, const int *a2, + const int *b1, const int *b2) { return - max(a1[0],b1[0]) < min(a1[1],b1[1]) || - max(a2[0],b2[0]) < min(a2[1],b2[1]) || - lineintersection(a1[0], a2[0], b1[0], b2[0]) || - lineintersection(a1[1], a2[1], b1[1], b2[1]); + MAX(a1[0], b1[0]) < MIN(a1[1], b1[1]) || + MAX(a2[0], b2[0]) < MIN(a2[1], b2[1]) || + LINE_INTERSECTION(a1[0], a2[0], b1[0], b2[0]) || + LINE_INTERSECTION(a1[1], a2[1], b1[1], b2[1]); } -#define meaningful_face(i,j) \ +#define MEANINGFUL_FACE(i,j) \ !((a1[i] ==INT_MIN) && (b1[j] ==INT_MIN)) && \ !((a1[i+1]==INT_MAX) && (b1[j+1]==INT_MAX)) @@ -253,10 +258,10 @@ void findconnections(int n, int *pts[4], /* Completely matching faces */ - if (a1[i]==b1[j] && a1[i+1]==b1[j+1] && - a2[i]==b2[j] && a2[i+1]==b2[j+1]){ + if (((a1[i] == b1[j]) && (a1[i+1] == b1[j+1])) && + ((a2[i] == b2[j]) && (a2[i+1] == b2[j+1]))) { - if (meaningful_face(i,j)){ + 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; @@ -286,7 +291,8 @@ void findconnections(int n, int *pts[4], else{ /* Find new intersection */ - if (lineintersection(a1[i+1],a2[i+1],b1[j+1],b2[j+1])) { + if (LINE_INTERSECTION(a1[i+1], a2[i+1], + b1[j+1], b2[j+1])) { itop[j+1] = out->number_of_nodes++; /* store point numbers of intersecting lines */ @@ -309,7 +315,7 @@ void findconnections(int n, int *pts[4], /* Add face to list of faces if no INT_MIN or * INT_MAX appear in a or b. */ - if (meaningful_face(i,j)){ + if (MEANINGFUL_FACE(i,j)) { /* Even indices refer to space between cells, @@ -352,7 +358,7 @@ void findconnections(int n, int *pts[4], for(j=0;j Date: Wed, 20 Jun 2012 09:51:20 +0000 Subject: [PATCH 50/72] Fully parenthesise the MEANINGFUL_FACE macro. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/facetopology.c b/facetopology.c index d70d5025..0b60a515 100644 --- a/facetopology.c +++ b/facetopology.c @@ -198,9 +198,9 @@ faceintersection(const int *a1, const int *a2, } -#define MEANINGFUL_FACE(i,j) \ - !((a1[i] ==INT_MIN) && (b1[j] ==INT_MIN)) && \ - !((a1[i+1]==INT_MAX) && (b1[j+1]==INT_MAX)) +#define MEANINGFUL_FACE(i, j) \ + (! ((a1[i] == INT_MIN) && (b1[j] == INT_MIN)) && \ + ! ((a1[i+1] == INT_MAX) && (b1[j+1] == INT_MAX))) /* work should be pointer to 2n ints initialised to zero . */ From 406016998bac2190c9cf5a738aaacad17f020c4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 13:02:57 +0000 Subject: [PATCH 51/72] Continue reorganisation to promote readability. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Specifically: * Move linearindex() ahead of all existing functions and use it in compute_cell_index(). * compute_cell_index(): Insert white-space and comments. * process_vertical_faces(): Add comments to describe stages in the process of computing new connections (faces). Signed-off-by: Bård Skaflestad --- preprocess.c | 53 ++++++++++++++++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/preprocess.c b/preprocess.c index 426599d5..e626b37e 100644 --- a/preprocess.c +++ b/preprocess.c @@ -63,6 +63,13 @@ process_horizontal_faces(int **intersections, int *plist, struct processed_grid *out); +static int +linearindex(const int dims[3], int i, int j, int k) +{ + return i + dims[0]*(j + dims[1]*k); +} + + /*----------------------------------------------------------------- Given a vector with k index running faster than i running faster than j, and Cartesian dimensions , find pointers to the @@ -95,19 +102,22 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[]) */ static void -compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len) +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= dims[0])) || /* 'i' outside [0, dims[0]) */ + ((j < 0) || (j >= dims[1]))) { /* 'j' outside [0, dims[1]) */ + + for (k = 0; k < len; k += 2) { + neighbors[k] = -1; /* Neighbour is outside domain */ } - }else{ - for(k=0; k */ /* 0 2 2 3 */ @@ -223,15 +233,23 @@ process_vertical_faces(int direction, num_intersections = out->number_of_nodes - out->number_of_nodes_on_pillars; - findconnections(2*nz+2, cornerpts, - *intersections+4*num_intersections, + /* Establish new connections (faces) along pillar pair. */ + findconnections(2*nz + 2, cornerpts, + *intersections + 4*num_intersections, work, out); + /* Start of ->face_neighbors[] for this set of connections. */ ptr = out->face_neighbors + 2*startface; + + /* Total number of cells (both sides) connected by this + * set of connections (faces). */ len = 2*out->number_of_faces - 2*startface; - compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr, len); - compute_cell_index(out->dimensions, i, j, ptr+1, len); + /* Derive inter-cell connectivity (i.e. ->face_neighbors) + * of global (uncompressed) cells for this set of + * connections (faces). */ + compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr , len); + compute_cell_index(out->dimensions, i , j , ptr + 1, len); /* Tag the new faces */ f = startface; @@ -242,11 +260,6 @@ process_vertical_faces(int direction, } } -static int linearindex(const int dims[3], int i, int j, int k) -{ - return i + dims[0]*(j + dims[1]*k); -} - /*----------------------------------------------------------------- For each horizontal face (i.e. k constant), From 736b623d4b623ab51f90f7a4fd9aa91794ef362a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 13:03:51 +0000 Subject: [PATCH 52/72] findconnections(): Insert white-space for readability. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/facetopology.c b/facetopology.c index 0b60a515..87131179 100644 --- a/facetopology.c +++ b/facetopology.c @@ -231,20 +231,23 @@ void findconnections(int n, int *pts[4], for (i = 0; i < 4; i++) { intersect[i] = -1; } - for (i = 0; i Date: Wed, 20 Jun 2012 13:25:11 +0000 Subject: [PATCH 53/72] linearindex(): Assert component indices in range. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit While here, split a long line and sort includes. Signed-off-by: Bård Skaflestad --- preprocess.c | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/preprocess.c b/preprocess.c index e626b37e..b61ed474 100644 --- a/preprocess.c +++ b/preprocess.c @@ -32,12 +32,13 @@ along with OpenRS. If not, see . */ -#include -#include -#include #include -#include #include +#include +#include +#include +#include + #include "preprocess.h" #include "uniquepoints.h" #include "facetopology.h" @@ -66,6 +67,14 @@ process_horizontal_faces(int **intersections, static int linearindex(const int dims[3], int i, int j, int k) { + assert (0 <= i); + assert (0 <= j); + assert (0 <= k); + + assert (i < dims[0]); + assert (j < dims[1]); + assert (k < dims[2]); + return i + dims[0]*(j + dims[1]*k); } @@ -75,9 +84,9 @@ linearindex(const int dims[3], int i, int j, int k) faster than j, and Cartesian dimensions , find pointers to the (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of field. */ -static void igetvectors(int dims[3], int i, int j, int *field, int *v[]) +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; From 40848dc0b0836c23f82f0df868cd0c855775a08c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 16:38:07 +0000 Subject: [PATCH 54/72] checkmemory(): Group (re)allocations according to the size-parameter used to effect the allocations. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Furthermore: Capture all allocations that succeed even if some other allocation fails. Signed-off-by: Bård Skaflestad --- preprocess.c | 59 ++++++++++++++++++++++++++++------------------------ 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/preprocess.c b/preprocess.c index b61ed474..0e461d7b 100644 --- a/preprocess.c +++ b/preprocess.c @@ -138,48 +138,53 @@ compute_cell_index(const int dims[3], int i, int j, static int checkmemory(int nz, struct processed_grid *out, int **intersections) { + int r, m, n, ok; /* Ensure there is enough space */ - int r = (2*nz+2)*(2*nz+2); - int m = out->m; - int n = out->n; + r = (2*nz + 2) * (2*nz + 2); + m = out->m; + n = out->n; - if(out->number_of_faces + r > m){ + 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){ + 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; - } + ok = m == out->m; + if (! ok) { + void *p1, *p2, *p3, *p4; + + p1 = realloc(*intersections , 4*m * sizeof **intersections); + p2 = realloc(out->face_neighbors, 2*m * sizeof *out->face_neighbors); + p3 = realloc(out->face_ptr , (m+1) * sizeof *out->face_ptr); + p4 = realloc(out->face_tag , 1*m * sizeof *out->face_tag); + + if (p1 != NULL) { *intersections = p1; } + if (p2 != NULL) { out->face_neighbors = p2; } + if (p3 != NULL) { out->face_ptr = p3; } + if (p4 != NULL) { out->face_tag = p4; } + + ok = (p1 != NULL) && (p2 != NULL) && (p3 != NULL) && (p4 != NULL); + + if (ok) { out->m = m; } } - if (m != out->m || n != out->n) { - void *p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes); - void *p2 = realloc(out->face_ptr, (m+1) * sizeof *out->face_ptr); - void *p3 = realloc(out->face_tag, m * sizeof *out->face_tag); + if (ok && (n != out->n)) { + void *p1; - if (p1 && p2 && p3) { + p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes); + + ok = p1 != NULL; + + if (ok) { out->face_nodes = p1; - out->face_ptr = p2; - out->face_tag = p3; - } else { - return 0; + out->n = n; } - - out->m = m; - out->n = n; } - return 1; + return ok; } /*----------------------------------------------------------------- From a606cd1e76b3cf1e59db24e7fa37563891bc8bb2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 17:45:49 +0000 Subject: [PATCH 55/72] Up-case macros (min->MIN, max->MAX). MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/preprocess.c b/preprocess.c index 0e461d7b..178485fa 100644 --- a/preprocess.c +++ b/preprocess.c @@ -44,8 +44,8 @@ #include "facetopology.h" -#define min(i,j) ((i)<(j) ? (i) : (j)) -#define max(i,j) ((i)>(j) ? (i) : (j)) +#define MIN(i,j) ((i)<(j) ? (i) : (j)) +#define MAX(i,j) ((i)>(j) ? (i) : (j)) static void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len); @@ -87,10 +87,10 @@ linearindex(const int dims[3], int i, int j, int k) 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); @@ -146,10 +146,10 @@ checkmemory(int nz, struct processed_grid *out, int **intersections) n = out->n; if (out->number_of_faces + r > m) { - m += max(m*0.5, 2*r); + 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); + n += MAX(n*0.5, 12*r); } ok = m == out->m; From 1d1bf4b4f7477451483c425d8a258e222e97a2fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 17:53:53 +0000 Subject: [PATCH 56/72] Fix misprint in diagnostic message and split two long lines. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/preprocess.c b/preprocess.c index 178485fa..c3591db0 100644 --- a/preprocess.c +++ b/preprocess.c @@ -215,11 +215,14 @@ process_vertical_faces(int direction, int num_intersections; int *ptr; int len; + for (j=0; j Date: Wed, 20 Jun 2012 17:58:44 +0000 Subject: [PATCH 57/72] process_vertical_faces(): Hoist invariant assignments above loop. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/preprocess.c b/preprocess.c index c3591db0..b455540a 100644 --- a/preprocess.c +++ b/preprocess.c @@ -216,6 +216,10 @@ process_vertical_faces(int direction, int *ptr; int len; + d[0] = 2 * (nx + 0); + d[1] = 2 * (ny + 0); + d[2] = 2 * (nz + 1); + for (j=0; j Date: Wed, 20 Jun 2012 18:07:32 +0000 Subject: [PATCH 58/72] checkmemory(): Expand on the role of "r" size estimate for the number of additional connections (faces) along a faulted stack. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/preprocess.c b/preprocess.c index b455540a..f01d9d26 100644 --- a/preprocess.c +++ b/preprocess.c @@ -140,7 +140,10 @@ checkmemory(int nz, struct processed_grid *out, int **intersections) { int r, m, n, ok; - /* Ensure there is enough space */ + /* Ensure there is enough space to manage the (pathological) case + * of every single cell on one side of a fault connecting to all + * cells on the other side of the fault (i.e., an all-to-all cell + * connectivity pairing). */ r = (2*nz + 2) * (2*nz + 2); m = out->m; n = out->n; From ddb1ae57e086109030cb2656c3b4771373fc3aa6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 18:13:13 +0000 Subject: [PATCH 59/72] Document the (internal) role of the "m" and "n" fields of the "processed_grid". MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/preprocess.h b/preprocess.h index e0453cc4..7ab85928 100644 --- a/preprocess.h +++ b/preprocess.h @@ -53,8 +53,8 @@ extern "C" { /* Output structure holding grid topology */ struct processed_grid{ - int m; - int n; + int m; /** Upper bound on "number_of_faces" */ + int n; /** Upper bound on "number_of_nodes" */ int dimensions[3]; /* Cartesian dimension */ int number_of_faces; int *face_nodes; /* Nodes numbers of each face sequentially. */ From 7ebf14e2999c8997fc0e0c32e4b223bf17300620 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 18:21:02 +0000 Subject: [PATCH 60/72] Document the (slightly opaque) "number_of_nodes_on_pillars" field of the "processed_grid". MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocess.h b/preprocess.h index 7ab85928..28635b12 100644 --- a/preprocess.h +++ b/preprocess.h @@ -63,7 +63,7 @@ extern "C" { enum face_tag *face_tag; int number_of_nodes; - int number_of_nodes_on_pillars; + int number_of_nodes_on_pillars; /** Total number of unique cell vertices that lie on pillars. */ double *node_coordinates; /* 3 doubles per node, sequentially */ int number_of_cells; /* number of active cells */ From e4b24a32628ecc8b3f24b5ddcb806cd875a91bbd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 18:22:20 +0000 Subject: [PATCH 61/72] Remove trailing white-space. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/preprocess.h b/preprocess.h index 28635b12..9592ce92 100644 --- a/preprocess.h +++ b/preprocess.h @@ -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, @@ -61,18 +61,18 @@ extern "C" { int *face_ptr; /* Start position for each face in face_nodes. */ int *face_neighbors; /* Global cell numbers. 2 ints per face sequentially */ enum face_tag *face_tag; - - int number_of_nodes; + + int number_of_nodes; int number_of_nodes_on_pillars; /** Total number of unique cell vertices that lie on pillars. */ double *node_coordinates; /* 3 doubles per node, sequentially */ - + int number_of_cells; /* number of active cells */ int *local_cell_index; /* Global to local map */ }; - void process_grdecl (const struct grdecl *g, - double tol, + void process_grdecl (const struct grdecl *g, + double tol, struct processed_grid *out); void free_processed_grid(struct processed_grid *g); From a3dd542ea4f806f2e0d7772b49f4ea351d501852 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 18:41:29 +0000 Subject: [PATCH 62/72] checkmemory(): Use integer-only operations when calculating array sizes. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/preprocess.c b/preprocess.c index f01d9d26..b38b2039 100644 --- a/preprocess.c +++ b/preprocess.c @@ -149,10 +149,10 @@ checkmemory(int nz, struct processed_grid *out, int **intersections) n = out->n; if (out->number_of_faces + r > m) { - m += MAX(m*0.5, 2*r); + m += MAX(m / 2, 2 * r); } if (out->face_ptr[out->number_of_faces] + 6*r > n) { - n += MAX(n*0.5, 12*r); + n += MAX(n / 2, 12 * r); } ok = m == out->m; From 18930081cef661cd1ec576aac712e21432220c2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 18:42:19 +0000 Subject: [PATCH 63/72] process_vertical_faces(): Explicitly assert that "direction" is zero or one only. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/preprocess.c b/preprocess.c index b38b2039..60ef6720 100644 --- a/preprocess.c +++ b/preprocess.c @@ -219,6 +219,8 @@ process_vertical_faces(int direction, int *ptr; int len; + assert ((direction == 0) || (direction == 1)); + d[0] = 2 * (nx + 0); d[1] = 2 * (ny + 0); d[2] = 2 * (nz + 1); From e72cc0d60052ddbdeb6271e257e31a7986e2b25c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 20 Jun 2012 18:52:49 +0000 Subject: [PATCH 64/72] process_vertical_faces(): Insert a bit of white-space in loop-condition for readability purposes. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- preprocess.c | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/preprocess.c b/preprocess.c index 60ef6720..b5c74347 100644 --- a/preprocess.c +++ b/preprocess.c @@ -225,8 +225,8 @@ process_vertical_faces(int direction, d[1] = 2 * (ny + 0); d[2] = 2 * (nz + 1); - for (j=0; j Date: Mon, 25 Jun 2012 20:09:39 +0000 Subject: [PATCH 65/72] Make special provisions to guarantee positive cell volumes in models featuring a left-handed coordinate system. Implementation modelled on the "processGRDECL" function of MRST. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Specifically: * Use a simple triple-product characterisation of left-handed coordinate systems. * Reflect Y coordinates about XZ plane to guarantee right-handed systems during facetopology() processing. * Restore original Y coordinates upon completion of process_grdecl() lest subsequent geometry processing fail to produce correct results. * Reverse face-node order if an odd number of sign changes have been applied to the node coordinates (i.e., if *either* the Y *or* the Z--but not both--coordinates have been flipped). Signed-off-by: Bård Skaflestad --- preprocess.c | 160 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 158 insertions(+), 2 deletions(-) diff --git a/preprocess.c b/preprocess.c index b5c74347..6e0e5db2 100644 --- a/preprocess.c +++ b/preprocess.c @@ -220,7 +220,7 @@ process_vertical_faces(int direction, int len; assert ((direction == 0) || (direction == 1)); - + d[0] = 2 * (nx + 0); d[1] = 2 * (ny + 0); d[2] = 2 * (nz + 1); @@ -618,6 +618,136 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum, } +static void +ind2sub(const size_t nx, + const size_t ny, + const size_t nz, + size_t c , + size_t *i, size_t *j, size_t *k) +{ + assert (c < (nx * ny * nz)); + +#if !defined(NDEBUG) + (void) nz; +#endif + + *i = c % nx; c /= nx; + *j = c % ny; + *k = c / ny; +} + + +/* ---------------------------------------------------------------------- */ +static double +vert_size(const struct grdecl *in, + const size_t c , + const size_t off[8]) +/* ---------------------------------------------------------------------- */ +{ + size_t i, j, k, start; + double dz; + const double *zcorn; + + ind2sub(in->dims[0], in->dims[1], in->dims[2], c, &i, &j, &k); + + zcorn = in->zcorn; + start = (k * off[4]) + (j * off[2]) + (i * 2); + + for (k = 0, dz = 0.0; (! (fabs(dz) > 0)) && (k < 4); k++) { + dz = zcorn[start + off[k + 4]] - zcorn[start + off[k]]; + } + + return dz; +} + + +/* ---------------------------------------------------------------------- */ +static int +is_lefthanded(const struct grdecl *in) +/* ---------------------------------------------------------------------- */ +{ + int active, searching; + size_t nx, ny, nz, c; + size_t origin, imax, jmax; + size_t off[8]; + double dx[2], dy[2], dz, triple; + const double *pt_coord; + + nx = in->dims[0]; + ny = in->dims[1]; + nz = in->dims[2]; + + off[0] = 0; + off[1] = off[0] + 1; + off[2] = off[0] + (2 * nx); + off[3] = off[2] + 1; + off[4] = off[0] + ((2 * nx) * (2 * ny)); + off[5] = off[4] + 1; + off[6] = off[4] + (2 * nx); + off[7] = off[6] + 1; + + pt_coord = in->coord; + + origin = 0; + imax = (nx + 0) * 1 * (2 * 3); + jmax = (nx + 1) * (ny + 0) * (2 * 3); + + dx[0] = pt_coord[imax + 0] - pt_coord[origin + 0]; + dy[0] = pt_coord[imax + 1] - pt_coord[origin + 1]; + + dx[1] = pt_coord[jmax + 0] - pt_coord[origin + 0]; + dy[1] = pt_coord[jmax + 1] - pt_coord[origin + 1]; + + c = 0; dz = 0.0; + do { + active = (in->actnum == NULL) || (in->actnum[c] != 0); + + if (active) { + dz = vert_size(in, c, off); + } + + searching = ! (active && (fabs(dz) > 0.0)); + + c += 1; + } while (searching && (c < (nx * ny * nz))); + + assert (! searching); /* active && (fabs(dz) > 0) */ + + /* Compute vector triple product to distinguish left-handed (<0) + * from right-handed (>0) coordinate systems. */ + triple = dz * (dx[0]*dy[1] - dx[1]*dy[0]); + + assert (fabs(triple) > 0.0); + + return triple < 0.0; +} + + +/* ---------------------------------------------------------------------- */ +static void +reverse_face_nodes(struct processed_grid *out) +/* ---------------------------------------------------------------------- */ +{ + int f, t, *i, *j; + + for (f = 0; f < out->number_of_faces; f++) { + i = out->face_nodes + (out->face_ptr[f + 0] + 0); + j = out->face_nodes + (out->face_ptr[f + 1] - 1); + + assert (i <= j); + + while (i < j) { + t = *i; + *i = *j; + *j = t; + + i += 1; + j -= 1; + } + } +} + + /*----------------------------------------------------------------- Public interface */ @@ -628,7 +758,7 @@ void process_grdecl(const struct grdecl *in, struct grdecl g; size_t i; - int sign, error; + int sign, error, left_handed; int cellnum; int *actnum, *iptr; @@ -711,6 +841,16 @@ void process_grdecl(const struct grdecl *in, free (zcorn); free (actnum); + /* Determine if coordinate system is left handed or not. */ + left_handed = is_lefthanded(in); + if (left_handed) { + /* Reflect Y coordinates about XZ plane to create right-handed + * coordinate system whilst processing intersections. */ + for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) { + out->node_coordinates[i] = -out->node_coordinates[i]; + } + } + /* -----------------------------------------------------------------*/ /* Find face topology and face-to-cell connections */ @@ -765,6 +905,14 @@ void process_grdecl(const struct grdecl *in, free(out->local_cell_index); out->local_cell_index = global_cell_index; + /* Reflect Y coordinate back to original position if left-handed + * coordinate system was detected and handled earlier. */ + if (left_handed) { + for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) { + out->node_coordinates[i] = -out->node_coordinates[i]; + } + } + /* if sign==-1 in ZCORN preprocessing, the sign of the * z-coordinate need to change before we finish */ if (sign == -1) @@ -773,6 +921,14 @@ void process_grdecl(const struct grdecl *in, out->node_coordinates[i] *= sign; } + /* If an odd number of coordinate reflections were applied, the + * processing routines--especially facetopology()--will produce + * node orderings that lead to normals pointing from 2 to 1. + * Reverse nodes to reestablish expected normal direction (and + * positive cell volumes). */ + if (left_handed ^ (sign == -1)) { + reverse_face_nodes(out); + } } /*-------------------------------------------------------*/ From 71a1f93bf9fb703cdfa89cafd85b6a9bd5630213 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 25 Jun 2012 20:52:07 +0000 Subject: [PATCH 66/72] Reformat computeFaceTopology() declaration. No functional changes. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Bård Skaflestad --- facetopology.c | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/facetopology.c b/facetopology.c index 87131179..15c71c9a 100644 --- a/facetopology.c +++ b/facetopology.c @@ -60,12 +60,10 @@ /* Determine face topology first, then compute intersection. */ /* All intersections that occur are present in the final face geometry.*/ -static int *computeFaceTopology(int *a1, - int *a2, - int *b1, - int *b2, - int intersect[4], - int *faces) +static int * +computeFaceTopology(const int *a1, const int *a2, + const int *b1, const int *b2, + int intersect[4], int *faces) { int mask[8]; int k; From f6a7a488315823856058c3f3860aa3928c2f9130 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 26 Jun 2012 20:54:37 +0200 Subject: [PATCH 67/72] Compute_cell_geometry(): Don't print misleading diagnostic message. We don't know if we should emit a diagnostic concerning non-positive cell volumes until we have accumulated all tetrahedral contributions. While here, amend the diagnostic message to also report the cell in which non-positive volumes are detected. --- geometry.c | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/geometry.c b/geometry.c index 27e4cfdd..c4c06394 100644 --- a/geometry.c +++ b/geometry.c @@ -131,6 +131,12 @@ compute_cell_geometry(int ndims, double *coords, double ccell[3]; double cface[3] = {0}; const double twothirds = 0.666666666666666666666666666667; + + int ndigits; + + ndigits = ((int) (log(ncells) / log(10.0))) + 1; + + for (c=0; c0)){ - fprintf(stderr, "Internal error in mex_compute_geometry: negative volume\n"); - } /* face centroid of triangle */ for (i=0; i 0.0)) { + fprintf(stderr, + "Internal error in mex_compute_geometry(%*d): " + "negative volume\n", ndigits, c); + } + for (i=0; i Date: Tue, 26 Jun 2012 20:57:00 +0200 Subject: [PATCH 68/72] Remove some long-disabled code that only served to confuse the reader. --- geometry.c | 4 ---- 1 file changed, 4 deletions(-) diff --git a/geometry.c b/geometry.c index c4c06394..5af49087 100644 --- a/geometry.c +++ b/geometry.c @@ -198,14 +198,10 @@ compute_cell_geometry(int ndims, double *coords, cross(u,v,w); - - tet_volume = 0; for(i=0; i Date: Tue, 26 Jun 2012 20:59:19 +0200 Subject: [PATCH 69/72] Lower constant tet_volume multiplication below accumulation loop. --- geometry.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/geometry.c b/geometry.c index 5af49087..37e11008 100644 --- a/geometry.c +++ b/geometry.c @@ -200,8 +200,9 @@ compute_cell_geometry(int ndims, double *coords, tet_volume = 0; for(i=0; i Date: Mon, 25 Jun 2012 19:59:10 +0000 Subject: [PATCH 70/72] Use more traditional comparison: "a!=b" rather than "!(a==b)" MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit git-svn-id: http://svn.sintef.no/simmatlab/branches/mrst-reorg/mex/libgeometry@9310 ea62cf4e-aff8-4dd6-bb37-af94872dbd4c Signed-off-by: Bård Skaflestad --- geometry.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geometry.c b/geometry.c index 37e11008..21c1dafd 100644 --- a/geometry.c +++ b/geometry.c @@ -212,7 +212,7 @@ compute_cell_geometry(int ndims, double *coords, if(subnormal_sign<0){ tet_volume*=-1.0; } - if(!(neighbors[2*face+0]==c)){ + if (neighbors[2*face + 0] != c) { tet_volume*=-1.0; } From 37d7f4d8e2331fb8fde712d265e7080ce39bf5e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 25 Jun 2012 20:00:25 +0000 Subject: [PATCH 71/72] Negate quantity by unary minus rather than multiplication by "-1". MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit git-svn-id: http://svn.sintef.no/simmatlab/branches/mrst-reorg/mex/libgeometry@9311 ea62cf4e-aff8-4dd6-bb37-af94872dbd4c Signed-off-by: Bård Skaflestad --- geometry.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geometry.c b/geometry.c index 21c1dafd..4d47c307 100644 --- a/geometry.c +++ b/geometry.c @@ -210,10 +210,10 @@ compute_cell_geometry(int ndims, double *coords, } if(subnormal_sign<0){ - tet_volume*=-1.0; + tet_volume = - tet_volume; } if (neighbors[2*face + 0] != c) { - tet_volume*=-1.0; + tet_volume = - tet_volume; } volume += tet_volume; From 619ce079385a8d6ef215950d0c11c8cd87654335 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 25 Jun 2012 20:01:02 +0000 Subject: [PATCH 72/72] Insert white-space for readability. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit git-svn-id: http://svn.sintef.no/simmatlab/branches/mrst-reorg/mex/libgeometry@9312 ea62cf4e-aff8-4dd6-bb37-af94872dbd4c Signed-off-by: Bård Skaflestad --- geometry.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geometry.c b/geometry.c index 4d47c307..b4dfc211 100644 --- a/geometry.c +++ b/geometry.c @@ -209,7 +209,7 @@ compute_cell_geometry(int ndims, double *coords, subnormal_sign += w[i]*fnormals[3*face+i]; } - if(subnormal_sign<0){ + if (subnormal_sign < 0.0) { tet_volume = - tet_volume; } if (neighbors[2*face + 0] != c) {