Merge pull request #1 from bska/master

Bring in more resilient grid processing and geometry calculation

Makefile.am

@@ -57,8 +57,6 @@ opm/core/grid/cornerpoint_grid.c			\
 opm/core/grid/cpgpreprocess/facetopology.c		\
 opm/core/grid/cpgpreprocess/geometry.c			\
 opm/core/grid/cpgpreprocess/preprocess.c		\
-opm/core/grid/cpgpreprocess/readvector.cpp		\
-opm/core/grid/cpgpreprocess/sparsetable.c		\
 opm/core/grid/cpgpreprocess/uniquepoints.c		\
 opm/core/linalg/LinearSolverFactory.cpp			\
 opm/core/linalg/LinearSolverInterface.cpp		\
@@ -160,8 +158,6 @@ opm/core/grid/cpgpreprocess/facetopology.h		\
 opm/core/grid/cpgpreprocess/geometry.h			\
 opm/core/grid/cpgpreprocess/grdecl.h			\
 opm/core/grid/cpgpreprocess/preprocess.h		\
-opm/core/grid/cpgpreprocess/readvector.hpp		\
-opm/core/grid/cpgpreprocess/sparsetable.h		\
 opm/core/grid/cpgpreprocess/uniquepoints.h		\
 opm/core/linalg/LinearSolverFactory.hpp			\
 opm/core/linalg/LinearSolverInterface.hpp		\

opm/core/grid/cpgpreprocess/facetopology.c

@@ -32,20 +32,19 @@ You should have received a copy of the GNU General Public License
   along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <stdlib.h>
-#include <math.h>
-#include <string.h>
 #include <assert.h>
-#include <stdio.h>
 #include <limits.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 
 #include "preprocess.h"
-#include "sparsetable.h"
 #include "facetopology.h"
 
 /* No checking of input arguments in this code! */
-#define min(i,j) ((i)<(j) ? (i) : (j))
+#define MIN(i,j) ((i)<(j) ? (i) : (j))
-#define max(i,j) ((i)>(j) ? (i) : (j))
+#define MAX(i,j) ((i)>(j) ? (i) : (j))
 
 #define DEBUG 1
 
@@ -59,17 +58,18 @@ along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 /*------------------------------------------------------*/
 
 
-/* Determine face geometry first, then compute intersections. */
+/* Determine face topology first, then compute intersection. */
 /* All intersections that occur are present in the final face geometry.*/
-static int *computeFaceTopology(int *a1,
+static int *
-				int *a2,
+computeFaceTopology(const int *a1, const int *a2,
-				int *b1,
+                    const int *b1, const int *b2,
-				int *b2,
+                    int intersect[4], int *faces)
-				int intersect[4],
-				int *faces)
 {
-  int mask[8] = {-1};
+    int mask[8];
-  int k, *f;
+    int k;
+    int *f;
+
+    for (k = 0; k < 8; k++) { mask[k] = -1; }
 
     /* Which pillar points should we use? */
     if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; }
@@ -146,7 +146,6 @@ static int *computeFaceTopology(int *a1,
 
 #if DEBUG>1
     /* Check for repeated nodes:*/
-  {
     int i;
     fprintf(stderr, "face: ");
     for (i=0; i<8; ++i){
@@ -158,7 +157,6 @@ static int *computeFaceTopology(int *a1,
         }
     }
     fprintf(stderr, "\n");
-  }
 #endif
 
 
@@ -171,27 +169,37 @@ static int *computeFaceTopology(int *a1,
 
 
 
-/* a) If we assume that the index increase when z increase for
+/* a) If we assume that the index increase when z increase for each
-      each pillar (but only separately), we can use only the point indices.
+   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.
+   b) We assume input is preprocessed such that no intersections occur
-      This is convenient in the identification of (unique) intersections.
+   on the first and last lines, for instance by padding the grid with
+   extra cells.  This is convenient in the identification of (unique)
+   intersections.
 
 */
 
-#define lineintersection(a1,a2,b1,b2)(((a1>b1)&&(a2<b2))||((a1<b1)&&(a2>b2)))
+#define LINE_INTERSECTION(a1, a2, b1, b2)       \
-static int faceintersection(int *a1, int *a2, int *b1, int *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(a1[0], b1[0]) < MIN(a1[1], b1[1]) ||
-    max(a2[0],b2[0]) < min(a2[1],b2[1]) ||
+        MAX(a2[0], b2[0]) < MIN(a2[1], b2[1]) ||
-    lineintersection(a1[0], a2[0], b1[0], b2[0]) ||
+        LINE_INTERSECTION(a1[0], a2[0], b1[0], b2[0]) ||
-    lineintersection(a1[1], a2[1], b1[1], b2[1]);
+        LINE_INTERSECTION(a1[1], a2[1], b1[1], b2[1]);
 }
 
 
-#define meaningful_face !((a1[i]  ==INT_MIN) && (b1[j]  ==INT_MIN)) && \
+#define MEANINGFUL_FACE(i, j)                                   \
-                        !((a1[i+1]==INT_MAX) && (b1[j+1]==INT_MAX))
+    (! ((a1[i]   == INT_MIN) && (b1[j]   == INT_MIN)) &&        \
+     ! ((a1[i+1] == INT_MAX) && (b1[j+1] == INT_MAX)))
+
 
 /* work should be pointer to 2n ints initialised to zero . */
 void findconnections(int n, int *pts[4],
@@ -210,24 +218,34 @@ void findconnections(int n, int *pts[4],
     int *ibottom = work + n;
     int *f       = out->face_nodes + out->face_ptr[out->number_of_faces];
     int *c       = out->face_neighbors + 2*out->number_of_faces;
-  int *tmp;
 
     int k1  = 0;
     int k2  = 0;
 
     int i,j=0;
-  int intersect[4]= {-1};
+    int intersect[4];
+    int *tmp;
     /* for (i=0; i<2*n; work[i++]=-1); */
 
+    for (i = 0; i < 4; i++) { intersect[i] = -1; }
+
     for (i = 0; i < n - 1; ++i) {
-    if (a1[i] == a1[i+1] && a2[i] == a2[i+1]) continue;
+
+        /* pinched a-cell */
+        if ((a1[i] == a1[i + 1]) &&
+            (a2[i] == a2[i + 1])) {
+            continue;
+        }
 
 
+        while ((j < n-1) &&
+               ((b1[j] < a1[i + 1]) ||
+                (b2[j] < a2[i + 1])))
+        {
+            /* pinched b-cell */
+            if ((b1[j] == b1[j + 1]) &&
+                (b2[j] == b2[j + 1])) {
 
-
-    while(j<n-1 && (b1[j] < a1[i+1] || b2[j] < a2[i+1])){
-
-      if (b1[j] == b1[j+1] && b2[j] == b2[j+1]){
                 itop[j+1] = itop[j];
                 ++j;
                 continue;
@@ -241,14 +259,13 @@ void findconnections(int n, int *pts[4],
 
 
                 /* Completely matching faces */
-	if (a1[i]==b1[j] && a1[i+1]==b1[j+1] &&
+                if (((a1[i] == b1[j]) && (a1[i+1] == b1[j+1])) &&
-	    a2[i]==b2[j] && a2[i+1]==b2[j+1]){
+                    ((a2[i] == b2[j]) && (a2[i+1] == b2[j+1]))) {
 
-	  /* Add face to list of faces if not any first or last points are involved. */
+                    if (MEANINGFUL_FACE(i, j)) {
-	  if (meaningful_face){
 
-	    int cell_a = i%2 ? (i-1)/2 : -1;
+                        int cell_a = i%2 != 0 ? (i-1)/2 : -1;
-	    int cell_b = j%2 ? (j-1)/2 : -1;
+                        int cell_b = j%2 != 0 ? (j-1)/2 : -1;
 
                         if (cell_a != -1 || cell_b != -1){
                             *c++ = cell_a;
@@ -257,21 +274,17 @@ void findconnections(int n, int *pts[4],
                             /* face */
                             *f++ = a1[i];
                             *f++ = a2[i];
+
                             /* avoid duplicating nodes in pinched faces  */
-	      if (a2[i+1] != a2[i]) *f++ = a2[i+1]; 
+                            if (a2[i+1] != a2[i]) { *f++ = a2[i+1]; }
-	      if (a1[i+1] != a1[i]) *f++ = a1[i+1];
+                            if (a1[i+1] != a1[i]) { *f++ = a1[i+1]; }
 
                             out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
 
                         }
                         else{
                             ;
-	      /*
-		fprintf(stderr, 
-		      "Warning. For some reason we get face connecting void spaces\n");
-	      */
                         }
-
                     }
                 }
 
@@ -279,7 +292,8 @@ void findconnections(int n, int *pts[4],
                 else{
 
                     /* Find new intersection */
-	  if (lineintersection(a1[i+1],a2[i+1],b1[j+1],b2[j+1])) {
+                    if (LINE_INTERSECTION(a1[i+1], a2[i+1],
+                                          b1[j+1], b2[j+1])) {
                         itop[j+1] = out->number_of_nodes++;
 
                         /* store point numbers of intersecting lines */
@@ -300,15 +314,16 @@ void findconnections(int n, int *pts[4],
                     intersect[3] = itop[j+1];     /* i+1 x j+1 */
 
 
-	  /* Add face to list of faces if no INT_MIN or INT_MAX appear in a or b. */
+                    /* Add face to list of faces if no INT_MIN or
-	  if (meaningful_face){
+                     * INT_MAX appear in a or b. */
+                    if (MEANINGFUL_FACE(i,j)) {
 
                         /*
                            Even indices refer to space between cells,
                            odd indices refer to cells
                         */
-	    int cell_a = i%2 ? (i-1)/2 : -1;
+                        int cell_a = i%2 != 0 ? (i-1)/2 : -1;
-	    int cell_b = j%2 ? (j-1)/2 : -1;
+                        int cell_b = j%2 != 0 ? (j-1)/2 : -1;
 
 
 
@@ -322,32 +337,29 @@ void findconnections(int n, int *pts[4],
                         }
                         else{
                             ;
-	      /*
-		fprintf(stderr, 
-		      "Warning. For some reason we get face connecting void spaces\n");
-	      */
                         }
                     }
                 }
             }
 
             /* Update candidates for restart of j for in next i-iteration */
-      if (b1[j] < a1[i+1]) k1 = j;
+            if (b1[j] < a1[i+1]) { k1 = j; }
-      if (b2[j] < a2[i+1]) k2 = j;
+            if (b2[j] < a2[i+1]) { k2 = j; }
 
             j = j+1;
         }
 
 
 
-    /* Swap intersection records: top line of a[i,i+1] is bottom line of a[i+1,i+2] */
+        /* Swap intersection records: top line of a[i,i+1] is bottom
+         * line of a[i+1,i+2] */
         tmp = itop; itop = ibottom; ibottom = tmp;
 
         /* Zero out the "new" itop */
-    for(j=0;j<n; ++j) itop[j]=-1;
+        for(j=0;j<n; ++j) { itop[j]=-1; }
 
         /* Set j to appropriate start position for next i */
-    j = min(k1, k2);
+        j = MIN(k1, k2);
     }
 }
 

opm/core/grid/cpgpreprocess/facetopology.h

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

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

opm/core/grid/cpgpreprocess/grdecl.h

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

opm/core/grid/cpgpreprocess/mxgrdecl.c

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

opm/core/grid/cpgpreprocess/mxgrdecl.h

@@ -1,4 +1,4 @@
-//===========================================================================
+/*=========================================================================
 //
 // File: mxgrdecl.h
 //
@@ -10,7 +10,7 @@
 //
 // $Revision$
 //
-//===========================================================================
+//=======================================================================*/
 
 /*
   Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
@@ -32,11 +32,13 @@ You should have received a copy of the GNU General Public License
   along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPENRS_MXGRDECL_HEADER
+#ifndef OPM_MXGRDECL_HEADER
-#define OPENRS_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:                */

opm/core/grid/cpgpreprocess/preprocess.c

@@ -32,35 +32,65 @@ You should have received a copy of the GNU General Public License
   along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <stdlib.h>
-#include <math.h>
-#include <string.h>
 #include <assert.h>
-#include <stdio.h>
 #include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include "preprocess.h"
 #include "uniquepoints.h"
 #include "facetopology.h"
 
 
-#define min(i,j) ((i)<(j) ? (i) : (j))
+#define MIN(i,j) ((i)<(j) ? (i) : (j))
-#define max(i,j) ((i)>(j) ? (i) : (j))
+#define MAX(i,j) ((i)>(j) ? (i) : (j))
 
+static void
+compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len);
+
+static int
+checkmemory(int nz, struct processed_grid *out, int **intersections);
+
+static void
+process_vertical_faces(int direction,
+                       int **intersections,
+                       int *plist, int *work,
+                       struct processed_grid *out);
+
+static void
+process_horizontal_faces(int **intersections,
+                         int *plist,
+                         struct processed_grid *out);
+
+static int
+linearindex(const int dims[3], int i, int j, int k)
+{
+    assert (0 <= i);
+    assert (0 <= j);
+    assert (0 <= k);
+
+    assert (i < dims[0]);
+    assert (j < dims[1]);
+    assert (k < dims[2]);
+
+    return i + dims[0]*(j + dims[1]*k);
+}
 
 
 /*-----------------------------------------------------------------
   Given a vector <field> with k index running faster than i running
   faster than j, and Cartesian dimensions <dims>, find pointers to the
   (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of
-  field.
+  field.  */
- */
+static void
-static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
+igetvectors(int dims[3], int i, int j, int *field, int *v[])
 {
-
+    int im = MAX(1,       i  ) - 1;
-  int im = max(1,       i  ) - 1;
+    int ip = MIN(dims[0], i+1) - 1;
-  int ip = min(dims[0], i+1) - 1;
+    int jm = MAX(1,       j  ) - 1;
-  int jm = max(1,       j  ) - 1;
+    int jp = MIN(dims[1], j+1) - 1;
-  int jp = min(dims[1], j+1) - 1;
 
     v[0] = field + dims[2]*(im + dims[0]* jm);
     v[1] = field + dims[2]*(im + dims[0]* jp);
@@ -76,23 +106,27 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
 /*-----------------------------------------------------------------
   Special purpose
 
-  Convert from k-index to Cartesian index i+nx*(j+ny*k) for every other
+  Convert from k-index to Cartesian index i+nx*(j+ny*k) for every
-  element in neighbors.
+  other element in neighbors.
 
 */
-static void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len)
+static void
+compute_cell_index(const int dims[3], int i, int j,
+                   int *neighbors, int len)
 {
-
     int k;
-  if (i<0 || i>=dims[0] || j<0 || j >= dims[1]){
+
+    if (((i < 0) || (i >= dims[0])) || /* 'i' outside [0, dims[0]) */
+        ((j < 0) || (j >= dims[1]))) { /* 'j' outside [0, dims[1]) */
+
         for (k = 0; k < len; k += 2) {
-      neighbors[k] = -1;
+            neighbors[k] = -1;  /* Neighbour is outside domain */
         }
-  }else{
+    }
+    else {
         for (k = 0; k < len; k += 2) {
             if (neighbors[k] != -1) {
-        int tmp = i + dims[0]*(j + dims[1]*neighbors[k]);
+                neighbors[k] = linearindex(dims, i, j, neighbors[k]);
-        neighbors[k] = tmp;
             }
         }
     }
@@ -100,52 +134,60 @@ static void compute_cell_index(const int dims[3], int i, int j, int *neighbors,
 
 
 /*-----------------------------------------------------------------
-  Ensure there's sufficient memory
+  Ensure there's sufficient memory */
- */
+static int
-static int checkmemeory(int nz, struct processed_grid *out, int **intersections)
+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
-  int r = (2*nz+2)*(2*nz+2);
+     * of every single cell on one side of a fault connecting to all
-  int m = out->m;
+     * cells on the other side of the fault (i.e., an all-to-all cell
-  int n = out->n;
+     * connectivity pairing). */
+    r = (2*nz + 2) * (2*nz + 2);
+    m = out->m;
+    n = out->n;
 
     if (out->number_of_faces +  r > m) {
-    m += max(m*0.5, 2*r);
+        m += MAX(m / 2,  2 * r);
     }
     if (out->face_ptr[out->number_of_faces] + 6*r > n) {
-    n += max(n*0.5, 12*r);
+        n += MAX(n / 2, 12 * r);
     }
 
-  if (m != out->m){
+    ok = m == out->m;
-    void *p1 = realloc(out->face_neighbors, 2*m * sizeof(*out->face_neighbors));
+    if (! ok) {
-    void *p2 = realloc(*intersections, 4*m * sizeof(**intersections));
+        void *p1, *p2, *p3, *p4;
-    if (p1 && p2) {
+
-      out->face_neighbors = p1;
+        p1 = realloc(*intersections     , 4*m   * sizeof **intersections);
-      *intersections = p2;
+        p2 = realloc(out->face_neighbors, 2*m   * sizeof *out->face_neighbors);
-    } else {
+        p3 = realloc(out->face_ptr      , (m+1) * sizeof *out->face_ptr);
-      return 0;
+        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) {
+    if (ok && (n != out->n)) {
-    void *p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes);
+        void *p1;
-    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) {
+        p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes);
+
+        ok = p1 != NULL;
+
+        if (ok) {
             out->face_nodes = p1;
-      out->face_ptr = p2;
-      out->face_tag = p3;
-    } else {
-      return 0;
-    }
-
-    out->m = m;
             out->n          = n;
         }
+    }
 
-  return 1;
+    return ok;
 }
 
 /*-----------------------------------------------------------------
@@ -157,41 +199,52 @@ static int checkmemeory(int nz, struct processed_grid *out, int **intersections)
   direction == 0 : constant-i faces.
   direction == 1 : constant-j faces.
 */
-static void process_vertical_faces(int direction,
+static void
+process_vertical_faces(int direction,
                        int **intersections,
                        int *plist, int *work,
                        struct processed_grid *out)
 {
     int i,j;
-  int d[3];
     int *cornerpts[4];
-
+    int d[3];
+    int f;
     enum face_tag tag[] = { LEFT, BACK };
-
+    int *tmp;
     int nx = out->dimensions[0];
     int ny = out->dimensions[1];
     int nz = out->dimensions[2];
+    int startface;
+    int num_intersections;
+    int *ptr;
+    int len;
 
-  int startface, num_intersections, f, len, *ptr;
+    assert ((direction == 0) || (direction == 1));
+
+    d[0] = 2 * (nx + 0);
+    d[1] = 2 * (ny + 0);
+    d[2] = 2 * (nz + 1);
 
     for (j = 0; j < ny + direction; ++j) {
-    for (i=0; i<nx+1-direction; ++i){
+        for (i = 0; i < nx + (1 - direction); ++i) {
 
-      if (!checkmemeory(nz, out, intersections)){
+            if (! checkmemory(nz, out, intersections)) {
-        fprintf(stderr, "Could not allocat enough space in process_vertical_faces\n");
+                fprintf(stderr,
+                        "Could not allocate enough space in "
+                        "process_vertical_faces()\n");
                 exit(1);
             }
 
             /* Vectors of point numbers */
-      d[0] = 2*nx; d[1] = 2*ny; d[2] = 2+2*nz;
+            igetvectors(d, 2*i + direction, 2*j + (1 - direction),
-      igetvectors(d, 2*i+direction, 2*j+1-direction, plist, cornerpts);
+                        plist, cornerpts);
 
             if (direction == 1) {
                 /* 1   3       0   1    */
                 /*       --->           */
                 /* 0   2       2   3    */
                 /* rotate clockwise     */
-        int *tmp     = cornerpts[1];
+                tmp          = cornerpts[1];
                 cornerpts[1] = cornerpts[0];
                 cornerpts[0] = cornerpts[2];
                 cornerpts[2] = cornerpts[3];
@@ -204,13 +257,21 @@ static void process_vertical_faces(int direction,
             num_intersections = out->number_of_nodes -
                 out->number_of_nodes_on_pillars;
 
+            /* Establish new connections (faces) along pillar pair. */
             findconnections(2*nz + 2, cornerpts,
                             *intersections + 4*num_intersections,
                             work, out);
 
+            /* Start of ->face_neighbors[] for this set of connections. */
             ptr = out->face_neighbors + 2*startface;
+
+            /* Total number of cells (both sides) connected by this
+             * set of connections (faces). */
             len = 2*out->number_of_faces - 2*startface;
 
+            /* Derive inter-cell connectivity (i.e. ->face_neighbors)
+             * of global (uncompressed) cells for this set of
+             * connections (faces). */
             compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr    , len);
             compute_cell_index(out->dimensions, i            , j          , ptr + 1, len);
 
@@ -223,11 +284,6 @@ static void process_vertical_faces(int direction,
     }
 }
 
-static int linearindex(const int dims[3], int i, int j, int k)
-{
-  return i+dims[0]*(j+dims[1]*k);
-}
-
 
 /*-----------------------------------------------------------------
   For each horizontal face (i.e. k constant),
@@ -240,7 +296,8 @@ static int linearindex(const int dims[3], int i, int j, int k)
   ACTNUM==0)
 
 */
-static void process_horizontal_faces(int **intersections,
+static void
+process_horizontal_faces(int **intersections,
                          int *plist,
                          struct processed_grid *out)
 {
@@ -252,20 +309,25 @@ static void process_horizontal_faces(int **intersections,
 
     int *cell  = out->local_cell_index;
     int cellno = 0;
-
+    int *f, *n, *c[4];
-  int *f, *n, *c[4], prevcell, thiscell;
+    int prevcell, thiscell;
+    int idx;
 
     /* dimensions of plist */
     int  d[3];
-  d[0] = 2*nx; d[1] = 2*ny; d[2] = 2+2*nz;
+    d[0] = 2*nx;
+    d[1] = 2*ny;
+    d[2] = 2+2*nz;
 
 
     for(j=0; j<ny; ++j) {
         for (i=0; i<nx; ++i) {
 
 
-      if (!checkmemeory(nz, out, intersections)){
+            if (! checkmemory(nz, out, intersections)) {
-        fprintf(stderr, "Could not allocat enough space in process_horizontal_faces\n");
+                fprintf(stderr,
+                        "Could not allocate enough space in "
+                        "process_horizontal_faces()\n");
                 exit(1);
             }
 
@@ -290,7 +352,7 @@ static void process_horizontal_faces(int **intersections,
 
                     /* If the pinch is a cell: */
                     if (k%2){
-            int idx = linearindex(out->dimensions, i,j,(k-1)/2);
+                        idx = linearindex(out->dimensions, i,j,(k-1)/2);
                         cell[idx] = -1;
                     }
                 }
@@ -348,39 +410,72 @@ static void process_horizontal_faces(int **intersections,
 */
 static void approximate_intersection_pt(int *L, double *c, double *pt)
 {
+    double a;
+    double z0, z1, z2, z3;
+    double b1, b2;
+    double x1, y1;
+    double x2, y2;
+    double z;
 
-  double z0 = c[3*L[0]+2];
+    /* no intersection on pillars expected here! */
-  double z1 = c[3*L[1]+2];
+    assert (L[0] != L[2]);
-  double z2 = c[3*L[2]+2];
+    assert (L[1] != L[3]);
-  double 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 */
+    if (fabs((z1 - z0) - (z3 - z2)) > 0.0) {
+
+        a = (z2 - z0) / ((z1 - z0) - (z3 - z2));
+
+    } else {
 
-  double a = (z2-z0)/(z1-z0 - (z3-z2));
-  if (isinf(a) || isnan(a)){
         a = 0;
+
     }
 
+    /* the corresponding z-coordinate is */
+    z =  z0*(1.0 - a) + z1*a;
 
-  pt[0] = c[3*L[0]+0]* (1.0-a) + c[3*L[1]+0]* a;
-  pt[1] = c[3*L[0]+1]* (1.0-a) + c[3*L[1]+1]* a;
-  pt[2] = c[3*L[0]+2]* (1.0-a) + c[3*L[1]+2]* a;
 
+    /* find point (x1, y1, z) on pillar 1 */
+    b1 = (z2 - z) / (z2 - z0);
+    b2 = (z - z0) / (z2 - z0);
+    x1 = c[3*L[0] + 0]*b1 + c[3*L[2] + 0]*b2;
+    y1 = c[3*L[0] + 1]*b1 + c[3*L[2] + 1]*b2;
+
+    /* find point (x2, y2, z) on pillar 2 */
+    b1 = (z - z3) / (z1 - z3);
+    b2 = (z1 - z) / (z1 - z3);
+    x2 = c[3*L[1] + 0]*b1 + c[3*L[3] + 0]*b2;
+    y2 = c[3*L[1] + 1]*b1 + c[3*L[3] + 1]*b2;
+
+    /* horizontal lines are by definition ON the bilinear surface
+       spanned by L0, L1, L2 and L3.  find point (x, y, z) on
+       horizontal line between point (x1, y1, z) and (x2, y2, z).*/
+    pt[0] = x1*(1.0 - a) + x2*a;
+    pt[1] = y1*(1.0 - a) + y2*a;
+    pt[2] = z;
 }
 
 /*-----------------------------------------------------------------
-  Compute x,y and z coordinates for points on each pillar.
+  Compute x,y and z coordinates for points on each pillar.  Then,
-  Then, append x,y and z coordinates for extra points on faults.
+  append x,y and z coordinates for extra points on faults.  */
- */
 static void
 compute_intersection_coordinates(int                   *intersections,
                                  struct processed_grid *out)
 {
     int n  = out->number_of_nodes;
     int np = out->number_of_nodes_on_pillars;
-
+    int    k;
-  int k, *itsct;
     double *pt;
-
+    int    *itsct = intersections;
-  /* Make sure the space allocated for nodes match the number of node. */
+    /* Make sure the space allocated for nodes match the number of
+     * node. */
     void *p = realloc (out->node_coordinates, 3*n*sizeof(double));
     if (p) {
         out->node_coordinates = p;
@@ -392,7 +487,6 @@ compute_intersection_coordinates(int                   *intersections,
 
     /* Append intersections */
     pt    = out->node_coordinates + 3*np;
-  itsct = intersections;
 
     for (k=np; k<n; ++k){
         approximate_intersection_pt(itsct, out->node_coordinates, pt);
@@ -403,6 +497,257 @@ compute_intersection_coordinates(int                   *intersections,
 }
 
 
+/* ------------------------------------------------------------------ */
+static int*
+copy_and_permute_actnum(int nx, int ny, int nz, const int *in, int *out)
+/* ------------------------------------------------------------------ */
+{
+    int i,j,k;
+    int *ptr = out;
+
+    /* Permute actnum such that values of each vertical stack of cells
+     * are adjacent in memory, i.e.,
+     *
+     *    out = [in(0,0,:), in(1,0,:),..., in(nx-1, ny-1,:)]
+     *
+     * in MATLAB pseudo-code.
+     */
+    if (in != NULL) {
+        for (j = 0; j < ny; ++j) {
+            for (i = 0; i < nx; ++i) {
+                for (k = 0; k < nz; ++k) {
+                    *ptr++ = in[i + nx*(j + ny*k)];
+                }
+            }
+        }
+    }
+    else {
+        /* No explicit ACTNUM.  Assume all cells active. */
+        for (i = 0; i < nx * ny * nz; i++) {
+            out[ i ] = 1;
+        }
+    }
+
+    return out;
+}
+
+/* ------------------------------------------------------------------ */
+static double*
+copy_and_permute_zcorn(int nx, int ny, int nz, const double *in,
+                       double sign, double *out)
+/* ------------------------------------------------------------------ */
+{
+    int i,j,k;
+    double *ptr = out;
+    /* Permute zcorn such that values of each vertical stack of cells
+     * are adjacent in memory, i.e.,
+
+     out = [in(0,0,:), in(1,0,:),..., in(2*nx-1, 2*ny-1,:)]
+
+     in Matlab pseudo-code.
+    */
+    for (j=0; j<2*ny; ++j){
+        for (i=0; i<2*nx; ++i){
+            for (k=0; k<2*nz; ++k){
+                *ptr++ = sign * in[i+2*nx*(j+2*ny*k)];
+            }
+        }
+    }
+    return out;
+}
+
+/* ------------------------------------------------------------------ */
+static int
+get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
+               const double *zcorn, int *error)
+/* ------------------------------------------------------------------ */
+{
+    /* Ensure that zcorn (i.e., depth) is strictly nondecreasing in
+       the k-direction.  This is required by the processign algorithm.
+
+       1) if  z(i,j,k) <= z(i,j,k+1) for all (i,j,k), return 1.0
+
+       2) if -z(i,j,k) <=-z(i,j,k+1) for all (i,j,k), return -1.0
+
+       3) if (1) and (2) fails, return -1.0, and set *error = 1.
+
+    */
+    int    sign;
+    int    i, j, k;
+    int    c1, c2;
+    double z1, z2;
+
+    for (sign = 1; sign>-2; sign = sign - 2)
+    {
+        *error = 0;
+
+        for (j=0; j<2*ny; ++j){
+            for (i=0; i<2*nx; ++i){
+                for (k=0; k<2*nz-1; ++k){
+                    z1 = sign*zcorn[i+2*nx*(j+2*ny*(k))];
+                    z2 = sign*zcorn[i+2*nx*(j+2*ny*(k+1))];
+
+                    c1 = i/2 + nx*(j/2 + ny*k/2);
+                    c2 = i/2 + nx*(j/2 + ny*(k+1)/2);
+
+                    if (((actnum == NULL) ||
+                         (actnum[c1] && actnum[c2]))
+                        && (z2 < z1)) {
+
+                        fprintf(stderr, "\nZCORN should be strictly "
+                                "nondecreasing along pillars!\n");
+                        *error = 1;
+                        goto end;
+                    }
+                }
+            }
+        }
+
+    end:
+        if (!*error){
+            break;
+        }
+    }
+
+    if (*error){
+        fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n"
+                "Grid definition may be broken\n");
+    }
+
+    return sign;
+}
+
+
+static void
+ind2sub(const size_t nx,
+        const size_t ny,
+        const size_t nz,
+        size_t       c ,
+        size_t *i, size_t *j, size_t *k)
+{
+    assert (c < (nx * ny * nz));
+
+#if !defined(NDEBUG)
+    (void) nz;
+#endif
+
+    *i = c % nx;  c /= nx;
+    *j = c % ny;
+    *k = c / ny;
+}
+
+
+/* ---------------------------------------------------------------------- */
+static double
+vert_size(const struct grdecl *in,
+          const size_t         c ,
+          const size_t         off[8])
+/* ---------------------------------------------------------------------- */
+{
+    size_t        i, j, k, start;
+    double        dz;
+    const double *zcorn;
+
+    ind2sub(in->dims[0], in->dims[1], in->dims[2], c, &i, &j, &k);
+
+    zcorn = in->zcorn;
+    start = (k * off[4]) + (j * off[2]) + (i * 2);
+
+    for (k = 0, dz = 0.0; (! (fabs(dz) > 0)) && (k < 4); k++) {
+        dz = zcorn[start + off[k + 4]] - zcorn[start + off[k]];
+    }
+
+    return dz;
+}
+
+
+/* ---------------------------------------------------------------------- */
+static int
+is_lefthanded(const struct grdecl *in)
+/* ---------------------------------------------------------------------- */
+{
+    int           active, searching;
+    size_t        nx, ny, nz, c;
+    size_t        origin, imax, jmax;
+    size_t        off[8];
+    double        dx[2], dy[2], dz, triple;
+    const double *pt_coord;
+
+    nx = in->dims[0];
+    ny = in->dims[1];
+    nz = in->dims[2];
+
+    off[0] = 0;
+    off[1] = off[0] + 1;
+    off[2] = off[0] + (2 * nx);
+    off[3] = off[2] + 1;
+    off[4] = off[0] + ((2 * nx) * (2 * ny));
+    off[5] = off[4] + 1;
+    off[6] = off[4] + (2 * nx);
+    off[7] = off[6] + 1;
+
+    pt_coord = in->coord;
+
+    origin = 0;
+    imax   = (nx + 0) * 1        * (2 * 3);
+    jmax   = (nx + 1) * (ny + 0) * (2 * 3);
+
+    dx[0] = pt_coord[imax + 0] - pt_coord[origin + 0];
+    dy[0] = pt_coord[imax + 1] - pt_coord[origin + 1];
+
+    dx[1] = pt_coord[jmax + 0] - pt_coord[origin + 0];
+    dy[1] = pt_coord[jmax + 1] - pt_coord[origin + 1];
+
+    c = 0;  dz = 0.0;
+    do {
+        active = (in->actnum == NULL) || (in->actnum[c] != 0);
+
+        if (active) {
+            dz = vert_size(in, c, off);
+        }
+
+        searching = ! (active && (fabs(dz) > 0.0));
+
+        c += 1;
+    } while (searching && (c < (nx * ny * nz)));
+
+    assert (! searching);       /* active && (fabs(dz) > 0) */
+
+    /* Compute vector triple product to distinguish left-handed (<0)
+     * from right-handed (>0) coordinate systems. */
+    triple = dz * (dx[0]*dy[1] - dx[1]*dy[0]);
+
+    assert (fabs(triple) > 0.0);
+
+    return triple < 0.0;
+}
+
+
+/* ---------------------------------------------------------------------- */
+static void
+reverse_face_nodes(struct processed_grid *out)
+/* ---------------------------------------------------------------------- */
+{
+    int f, t, *i, *j;
+
+    for (f = 0; f < out->number_of_faces; f++) {
+        i = out->face_nodes + (out->face_ptr[f + 0] + 0);
+        j = out->face_nodes + (out->face_ptr[f + 1] - 1);
+
+        assert (i <= j);
+
+        while (i < j) {
+            t  = *i;
+            *i = *j;
+            *j =  t;
+
+            i += 1;
+            j -= 1;
+        }
+    }
+}
+
+
 /*-----------------------------------------------------------------
   Public interface
 */
@@ -410,105 +755,39 @@ void process_grdecl(const struct grdecl   *in,
                     double                tolerance,
                     struct processed_grid *out)
 {
-  int i,j,k;
-
-  int nx = in->dims[0];
-  int ny = in->dims[1];
-  int nz = in->dims[2];
-
-  int nc = nx*ny*nz;
     struct grdecl g;
 
-  int    *actnum, *iptr, sign, error;
+    size_t i;
-  double *zcorn , *dptr;
+    int    sign, error, left_handed;
+    int    cellnum;
 
-  int *intersections, *work, *plist, *ptr, *global_cell_index;
+    int    *actnum, *iptr;
+    int    *global_cell_index;
 
-  const int BIGNUM = 64;
+    double *zcorn;
 
-  g.dims[0] = nx;
+    const size_t BIGNUM = 64;
-  g.dims[1] = ny;
+    const int    nx = in->dims[0];
-  g.dims[2] = nz;
+    const int    ny = in->dims[1];
-  actnum    = malloc (nc *     sizeof *actnum);
+    const int    nz = in->dims[2];
-  zcorn     = malloc (nc * 8 * sizeof *zcorn);
+    const size_t nc = ((size_t) nx) * ((size_t) ny) * ((size_t) nz);
 
-  /* Permute actnum */
+    /* internal work arrays */
-  iptr = actnum;
+    int    *work;
-  if (in->actnum != NULL) {
+    int    *plist;
-    for (j=0; j<ny; ++j){
+    int    *intersections;
-      for (i=0; i<nx; ++i){
-        for (k=0; k<nz; ++k){
-          *iptr++ = in->actnum[i+nx*(j+ny*k)];
-        }
-      }
-    }
-  } else {
-    for (i=0; i<nx*ny*nz; ++i){
-      *iptr++ = 1;
-    }
-  }
-  g.actnum = actnum;
-
-
-  /* HACK */
-  /* Check that ZCORN is strictly nodecreaing along pillars.  If */
-  /* not, check if -ZCORN is strictly nondecreasing.             */
-  for (sign = 1; sign>-2; sign = sign - 2){
-    error = 0;
-
-    /* Ensure that zcorn is strictly nondecreasing in k-direction */
-    for (j=0; j<2*ny; ++j){
-      for (i=0; i<2*nx; ++i){
-        for (k=0; k<2*nz-1; ++k){
-          double z1 = sign*in->zcorn[i+2*nx*(j+2*ny*(k))];
-          double z2 = sign*in->zcorn[i+2*nx*(j+2*ny*(k+1))];
-
-          int c1 = i/2 + nx*(j/2 + ny*(k/2));
-          int c2 = i/2 + nx*(j/2 + ny*((k+1)/2));
-
-          if (((in->actnum == NULL) || (in->actnum[c1] && in->actnum[c2])) && (z2 < z1)){
-            fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n");
-
-            error = 1;
-            goto end;
-          }
-        }
-      }
-    }
-
-  end:
-    if (!error){
-      break;
-    }
-  }
-
-  if (error){
-    fprintf(stderr, "Attempt to reverse sign in ZCORN failed.\n"
-                    "Grid definition may be broken\n");
-  }
-
-  /* Permute zcorn */
-  dptr = zcorn;
-  for (j=0; j<2*ny; ++j){
-    for (i=0; i<2*nx; ++i){
-      for (k=0; k<2*nz; ++k){
-        *dptr++ = sign*in->zcorn[i+2*nx*(j+2*ny*k)];
-      }
-    }
-  }
 
 
 
 
-  g.zcorn = zcorn;
+    /* -----------------------------------------------------------------*/
-  g.coord = in->coord;
+    /* Initialize output structure:
-
+       1) allocate space for grid topology (which may need to be
-
+          increased)
-  /* Code below assumes k index runs fastests, ie. that dimensions of
+       2) set Cartesian imensions
-     table is permuted to (dims[2], dims[0], dims[1]) */
+    */
-
+    out->m                = (int) (BIGNUM / 3);
-  out->m                = BIGNUM/3;
+    out->n                = (int) BIGNUM;
-  out->n                = BIGNUM;
 
     out->face_neighbors   = malloc( BIGNUM      * sizeof *out->face_neighbors);
     out->face_nodes       = malloc( out->n      * sizeof *out->face_nodes);
@@ -516,35 +795,74 @@ void process_grdecl(const struct grdecl   *in,
     out->face_tag         = malloc( out->m      * sizeof *out->face_tag);
     out->face_ptr[0]      = 0;
 
-  out->dimensions[0]    = g.dims[0];
+    out->dimensions[0]    = in->dims[0];
-  out->dimensions[1]    = g.dims[1];
+    out->dimensions[1]    = in->dims[1];
-  out->dimensions[2]    = g.dims[2];
+    out->dimensions[2]    = in->dims[2];
     out->number_of_faces  = 0;
     out->number_of_nodes  = 0;
     out->number_of_cells  = 0;
 
     out->node_coordinates = NULL;
-  out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index);
+    out->local_cell_index = malloc(nc * sizeof *out->local_cell_index);
 
 
-  /* internal */
-  intersections = malloc(BIGNUM* sizeof(*intersections));
-
-  /* internal */
-  work          = malloc(2* (2*nz+2)*   sizeof(*work));
-  for(i=0; i<4*(nz+1); ++i)   work[i] = -1;
-
-  /* Allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */
-  plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist);
-
 
     /* Do actual work here:*/
 
+    /* -----------------------------------------------------------------*/
+    /* For each pillar, compare zcorn values for adjacent cells to
+     * find the unique node z-coordinates specified by the input.
+     * While here, enumerate unique points and assign point numbers
+     * (in plist) for each cornerpoint cell. In other words, plist has
+     * 8 node numbers for each cornerpoint cell.*/
+
+    /* initialize grdecl structure "g" that will be processd by
+     * "finduniquepoints" */
+    g.dims[0] = in->dims[0];
+    g.dims[1] = in->dims[1];
+    g.dims[2] = in->dims[2];
+
+    actnum    = malloc (nc *     sizeof *actnum);
+    g.actnum  = copy_and_permute_actnum(nx, ny, nz, in->actnum, actnum);
+
+    zcorn     = malloc (nc * 8 * sizeof *zcorn);
+    sign      = get_zcorn_sign(nx, ny, nz, in->actnum, in->zcorn, &error);
+    g.zcorn   = copy_and_permute_zcorn(nx, ny, nz, in->zcorn, sign, zcorn);
+
+    g.coord   = in->coord;
+
+
+    /* allocate space for cornerpoint numbers plus INT_MIN (INT_MAX)
+     * padding */
+    plist = malloc(8 * (nc + ((size_t)nx)*((size_t)ny)) * sizeof *plist);
+
     finduniquepoints(&g, plist, tolerance, out);
 
     free (zcorn);
     free (actnum);
 
+    /* Determine if coordinate system is left handed or not. */
+    left_handed = is_lefthanded(in);
+    if (left_handed) {
+        /* Reflect Y coordinates about XZ plane to create right-handed
+         * coordinate system whilst processing intersections. */
+        for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) {
+            out->node_coordinates[i] = -out->node_coordinates[i];
+        }
+    }
+
+    /* -----------------------------------------------------------------*/
+    /* Find face topology and face-to-cell connections */
+
+    /* internal */
+    work = malloc(2 * ((size_t) (2*nz + 2)) * sizeof *work);
+    for(i = 0; i < ((size_t)4) * (nz + 1); ++i) { work[i] = -1; }
+
+    /* internal array to store intersections */
+    intersections = malloc(BIGNUM* sizeof(*intersections));
+
+
+
     process_vertical_faces   (0, &intersections, plist, work, out);
     process_vertical_faces   (1, &intersections, plist, work, out);
     process_horizontal_faces (   &intersections, plist,       out);
@@ -552,34 +870,65 @@ void process_grdecl(const struct grdecl   *in,
     free (plist);
     free (work);
 
+    /* -----------------------------------------------------------------*/
+    /* (re)allocate space for and compute coordinates of nodes that
+     * arise from intesecting cells (faults) */
     compute_intersection_coordinates(intersections, out);
 
     free (intersections);
 
-
+    /* -----------------------------------------------------------------*/
-  /* Convert to local cell numbers in face_neighbors */
+    /* Enumerate compressed cells:
-  ptr=out->face_neighbors;;
+       -make array [0...#cells-1] of global cell numbers
-  for (i=0; i<out->number_of_faces*2; ++i, ++ptr){
+       -make [0...nx*ny*nz-1] array of local cell numbers,
-    if (*ptr != -1){
+       lexicographically ordered, used to remap out->face_neighbors
-      *ptr = out->local_cell_index[*ptr];
+    */
-    }
+    global_cell_index = malloc(nc * sizeof *global_cell_index);
-  }
+    cellnum = 0;
-
+    for (i = 0; i < nc; ++i) {
-  /* Invert global-to-local map */
-  global_cell_index = malloc(out->number_of_cells *
-                                  sizeof (*global_cell_index));
-  for (i=0; i<nx*ny*nz; ++i){
         if (out->local_cell_index[i] != -1) {
-      global_cell_index[out->local_cell_index[i]] = i;
+            global_cell_index[cellnum] = (int) i;
+            out->local_cell_index[i]   = cellnum;
+            cellnum++;
         }
     }
+
+    /* Remap out->face_neighbors */
+    iptr = out->face_neighbors;
+    for (i = 0; i < ((size_t) 2) * out->number_of_faces; ++i, ++iptr) {
+        if (*iptr != -1){
+            *iptr = out->local_cell_index[*iptr];
+        }
+    }
+
+
     free(out->local_cell_index);
     out->local_cell_index = global_cell_index;
 
-  /* HACK: undo sign reversal in ZCORN prepprocessing */
+    /* Reflect Y coordinate back to original position if left-handed
-  for (i=2; i<3*out->number_of_nodes; i = i+3)
+     * coordinate system was detected and handled earlier. */
-    out->node_coordinates[i]=sign*out->node_coordinates[i];
+    if (left_handed) {
+        for (i = 1; i < ((size_t) 3) * out->number_of_nodes; i += 3) {
+            out->node_coordinates[i] = -out->node_coordinates[i];
+        }
+    }
 
+    /* if sign==-1 in ZCORN preprocessing, the sign of the
+     * z-coordinate need to change before we finish */
+    if (sign == -1)
+    {
+        for (i = 2; i < ((size_t) 3) * out->number_of_nodes; i += 3)
+            out->node_coordinates[i] *= sign;
+    }
+
+    /* If an odd number of coordinate reflections were applied, the
+     * processing routines--especially facetopology()--will produce
+     * node orderings that lead to normals pointing from 2 to 1.
+     * Reverse nodes to reestablish expected normal direction (and
+     * positive cell volumes). */
+    if (left_handed ^ (sign == -1)) {
+        reverse_face_nodes(out);
+    }
 }
 
 /*-------------------------------------------------------*/

opm/core/grid/cpgpreprocess/preprocess.h

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

opm/core/grid/cpgpreprocess/processgrid.c

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

opm/core/grid/cpgpreprocess/processgrid.m

@@ -1,12 +1,14 @@
-function varargout = processgrid(varargin)
+function G = processgrid(varargin)
 %Compute grid topology and geometry from pillar grid description.
 %
 % SYNOPSIS:
-%   G = processGRDECL(grdecl)
+%   G = processgrid(grdecl)
+%   G = processgrid(grdecl,ztol)
 %
 % PARAMETERS:
 %   grdecl - Raw pillar grid structure, as defined by function
 %            'readGRDECL', with fields COORDS, ZCORN and, possibly, ACTNUM.
+%   ztol   - tolerance for unique points 
 %
 % RETURNS:
 %   G      - Valid grid definition containing connectivity, cell
@@ -32,10 +34,33 @@ function varargout = processgrid(varargin)
 % $Date$
 % $Revision$
 
-% Build MEX edition of same.
+   G = processgrid_mex(varargin{:});
-%
+   G.griddim = 3;
-buildmex CFLAGS='$CFLAGS -Wall -fPIC' processgrid.c preprocess.c ...
+   G = splitDisconnectedGrid(G, false);   
-         uniquepoints.c facetopology.c sparsetable.c mxgrdecl.c
+end
 
-% Call MEX edition.
+function G = splitDisconnectedGrid(G, verbose)
-[varargout{1:nargout}] = processgrid(varargin{:});
+   % 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

opm/core/grid/cpgpreprocess/readvector.cpp

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

opm/core/grid/cpgpreprocess/readvector.hpp

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

opm/core/grid/cpgpreprocess/sparsetable.c

@@ -1,100 +0,0 @@
-/*===========================================================================
-//
-// File: sparsetable.c
-//
-// Created: Fri Jun 19 08:48:04 2009
-//
-// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
-//
-// $Date$
-//
-// $Revision$
-//
-//==========================================================================*/
-
-/*
-Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
-Copyright 2009, 2010 Statoil ASA.
-
-This file is part of The Open Reservoir Simulator Project (OpenRS).
-
-OpenRS is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or 
-(at your option) any later version.
-
-OpenRS is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include <stdlib.h>
-#include <math.h>
-#include <string.h>
-#include <assert.h>
-#include <stdio.h>
-
-
-#include "sparsetable.h"
-
-void free_sparse_table (sparse_table_t *tab)
-{
-  if(tab->ptr)  free(tab->ptr);
-  if(tab->data) free(tab->data);
-
-  free(tab);
-}
-
-sparse_table_t *malloc_sparse_table(int m, int n, int datasz)
-{
-  size_t alloc_sz;
-  sparse_table_t *tab = malloc(sizeof *tab);
-  tab->m = m; 
-  tab->n = n;
-  tab->position = 0;
-  if (!(tab->ptr  = malloc((m+1) * sizeof (*tab->ptr)))){
-    fprintf(stderr, "Could not allocate space for sparse ptr\n");
-    free_sparse_table(tab);
-    return NULL;
-  }
-  
-  alloc_sz  = datasz;
-  alloc_sz *= n;
-  if(!(tab->data = malloc(alloc_sz))){
-    fprintf(stderr, "Could not allocate space for sparse data(%d)\n", n);
-    free_sparse_table(tab);
-    return NULL;
-  }    
-  
-  return tab;
-}
-
-sparse_table_t *realloc_sparse_table(sparse_table_t *tab, int m, int n, int datasz)
-{  
-  void *p = realloc(tab->ptr, (m+1) * sizeof (*tab->ptr));
-  if (p){
-    tab->ptr = p;
-  }else{
-    fprintf(stderr, "Could not reallocate space for sparse ptr\n");
-    free_sparse_table(tab);
-    return NULL;
-  }
-  
-  p = realloc(tab->data, n * datasz);
-  if(p){
-    tab->data = p;
-  }else{
-    fprintf(stderr, "Could not reallocate space for sparse data(%d)\n", n);
-    free_sparse_table(tab);
-    return NULL;
-  }
-  
-  tab->m = m; 
-  tab->n = n;
-  
-  return tab;
-}

opm/core/grid/cpgpreprocess/sparsetable.h

@@ -1,53 +0,0 @@
-/*===========================================================================
-//
-// File: sparsetable.h
-//
-// Created: Fri Jun 19 08:47:45 2009
-//
-// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
-//
-// $Date$
-//
-// $Revision$
-//
-//===========================================================================*/
-
-/*
-Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
-Copyright 2009, 2010 Statoil ASA.
-
-This file is part of The Open Reservoir Simulator Project (OpenRS).
-
-OpenRS is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or 
-(at your option) any later version.
-
-OpenRS is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#ifndef OPENRS_SPARSETABLE_HEADER
-#define OPENRS_SPARSETABLE_HEADER
-
-typedef struct{
-  int m;        /* number of rows          */
-  int *ptr;     /* row pointer of size m+1 */
-  int position; /* first position in ptr that is not filled. */
-
-  int n;        /* size of data            */
-  void *data;   /* sparse table data       */
-  
-
-} sparse_table_t;
-
-void             free_sparse_table     (sparse_table_t *tab);
-sparse_table_t   *malloc_sparse_table  (int m, int n, int datasz);
-sparse_table_t   *realloc_sparse_table (sparse_table_t *tab, int m, int n, int datasz);
-
-#endif /* OPENRS_SPARSETABLE_HEADER */

opm/core/grid/cpgpreprocess/uniquepoints.c

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

opm/core/grid/cpgpreprocess/uniquepoints.h

@@ -32,15 +32,15 @@ You should have received a copy of the GNU General Public License
   along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPENRS_UNIQUEPOINTS_HEADER
+#ifndef OPM_UNIQUEPOINTS_HEADER
-#define OPENRS_UNIQUEPOINTS_HEADER
+#define OPM_UNIQUEPOINTS_HEADER
 
 int finduniquepoints(const struct grdecl *g,  /* input */
                      int                 *p,  /* for each z0 in zcorn, z0 = z[p0] */
                      double               t,  /* tolerance*/
                      struct processed_grid *out);
 
-#endif /* OPENRS_UNIQUEPOINTS_HEADER */
+#endif /* OPM_UNIQUEPOINTS_HEADER */
 
 /* Local Variables:    */
 /* c-basic-offset:4    */

opm/core/linalg/sparse_sys.h

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

tests/Makefile.am

@@ -19,7 +19,6 @@ test_column_extract				\
 test_lapack					\
 test_read_vag					\
 test_readpolymer				\
-test_readvector					\
 test_sf2p					\
 test_writeVtkData				\
 unit_test
@@ -49,8 +48,6 @@ test_lapack_LDADD         = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
 test_readpolymer_SOURCES  = test_readpolymer.cpp
 test_read_vag_SOURCES     = test_read_vag.cpp
 
-test_readvector_SOURCES   = test_readvector.cpp
-
 test_sf2p_SOURCES         = test_sf2p.cpp
 
 test_writeVtkData_SOURCES = test_writeVtkData.cpp

tests/test_readvector.cpp

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