Refactoring. Moving return struct around instead of a bunch of

auxillary variables.

facetopology.c

@@ -39,6 +39,7 @@ along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 #include <stdio.h>
 #include <limits.h>
 
+#include "preprocess.h"
 #include "sparsetable.h"
 #include "facetopology.h"
 
@@ -148,11 +149,9 @@ static int faceintersection(int *a1, int *a2, int *b1, int *b2)
 
 /* work should be pointer to 2n ints initialised to zero . */
 void findconnections(int n, int *pts[4],
-		     int *ptnumber,
 		     int *intersectionlist,
-		     int *neighbors,
 		     int *work,
-		     sparse_table_t *ftab)
+		     struct processed_grid *out)
 {
   /* vectors of point numbers for faces a(b) on pillar 1(2) */
   int *a1 = pts[0];
@@ -163,8 +162,8 @@ void findconnections(int n, int *pts[4],
   /* Intersection record for top line and bottomline of a */
   int *itop    = work;
   int *ibottom = work + n;
-  int *f       = (int*)ftab->data     + ftab->ptr[ftab->position];
+  int *f       = out->face_nodes + out->face_ptr[out->number_of_faces];
-  int *c       = neighbors + 2*ftab->position;
+  int *c       = out->face_neighbors + 2*out->number_of_faces;
 
   int k1  = 0;
   int k2  = 0;
@@ -213,7 +212,9 @@ void findconnections(int n, int *pts[4],
 	      *f++ = a1[i+1];
 	      *f++ = a2[i+1];
 	      *f++ = a2[i];
-	      ftab->ptr[++ftab->position] = f - (int*)ftab->data;
+
+	      out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
+
 	    }
 	    else{
 	      ;
@@ -231,7 +232,7 @@ void findconnections(int n, int *pts[4],
 
 	  /* Find new intersection */
 	  if (lineintersection(a1[i+1],a2[i+1],b1[j+1],b2[j+1])) {
-	    itop[j+1] = (*ptnumber)++;
+	    itop[j+1] = out->number_of_nodes++;
 
 	    /* store point numbers of intersecting lines */
 	    *intersectionlist++ = a1[i+1];
@@ -268,7 +269,8 @@ void findconnections(int n, int *pts[4],
 	      *c++ = cell_b;
 	      
 	      f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f);
-	      ftab->ptr[++ftab->position] = f - (int*)ftab->data;
+
+	      out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
 	    }
 	    else{
 	      ;

facetopology.h

@@ -37,11 +37,9 @@ along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 
 
 void findconnections(int n, int *pts[4],
-		     int *ptnumber,
 		     int *intersectionlist,
-		     int *neighbors,
                      int *work,
-                     sparse_table_t *ftab);
+		     struct processed_grid *out);
 
 #endif // OPENRS_FACETOPOLOGY_HEADER
 

preprocess.c

@@ -49,22 +49,6 @@ along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 
 
 
-/*-----------------------------------------------------------------
-  Given a z coordinate, find x and y coordinates on line defined by
-  coord.  Coord points to a vector of 6 doubles [x0,y0,z0,x1,y1,z1].
- */
-static void interpolate_pillar(const double *coord, double *pt)
-{
-  double a = (pt[2]-coord[2])/(coord[5]-coord[2]);
-  if (isinf(a) || isnan(a)){
-    a = 0;
-  }
-
-  pt[0]       = coord[0] + a*(coord[3]-coord[0]);
-  pt[1]       = coord[1] + a*(coord[4]-coord[1]);
-}
-
-
 /*-----------------------------------------------------------------
   Given a vector <field> with k index running faster than i running 
   faster than j, and Cartesian dimensions <dims>, find pointers to the 
@@ -119,38 +103,43 @@ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len
 /*-----------------------------------------------------------------
   Ensure there's sufficient memory
  */
-int checkmemeory(int nz, sparse_table_t *ftab, int **neighbors, int **intersections)
+int checkmemeory(int nz, struct processed_grid *out, int **intersections)
 {
 
   /* Ensure there is enough space */
   int r = (2*nz+2)*(2*nz+2);
-  int m = ftab->m;
+  int m = out->m;
-  int n = ftab->n;
+  int n = out->n;
-  if(ftab->position +  r> m){
+  if(out->number_of_faces +  r> m){
     m += max(m*0.5, 2*r);
   }
-  if (ftab->ptr[ftab->position] + 6*r > n){
+  if (out->face_ptr[out->number_of_faces] + 6*r > n){
     n += max(n*0.5, 12*r);
   }
       
-  if (m != ftab->m){
+  if (m != out->m){
-    void *p1 = realloc(*neighbors,     2*m * sizeof(**neighbors));
+    void *p1 = realloc(out->face_neighbors,     2*m * sizeof(*out->face_neighbors));
-    void *p2 = realloc(*intersections, 4*m * sizeof(**neighbors));
+    void *p2 = realloc(*intersections, 4*m * sizeof(**intersections));
     if (p1 && p2){
-      *neighbors     = p1;
+      out->face_neighbors     = p1;
       *intersections = p2;
     }else{
       return 0;
     }
   }
 
-  if (m != ftab->m || n != ftab->n){
+  if (m != out->m || n != out->n){
-    void *p = realloc_sparse_table(ftab, m, n, sizeof(int));
+    void *p1 = realloc(out->face_nodes, n*sizeof(double));
-    if (p){
+    void *p2 = realloc(out->face_ptr, m*sizeof(double));
-      ftab = p;
+    if (p1 && p2){
+      out->face_nodes = p1;
+      out->face_ptr = p2;
     }else{
       return 0;
     }
+
+    out->m = m;
+    out->n = n;
   }
 
   return 1;
@@ -165,24 +154,28 @@ int checkmemeory(int nz, sparse_table_t *ftab, int **neighbors, int **intersecti
   direction == 0 : constant-i faces.
   direction == 1 : constant-j faces.
  */
-void process_vertical_faces(const int dims[3], int direction, sparse_table_t *ftab,  
+void process_vertical_faces(int direction,  
-			    int **neighbors, int **intersections, int *npoints, 
+			   int **intersections, 
-			    int npillarpoints, int *plist, int *work)
+			    int *plist, int *work, 
+			    struct processed_grid *out)
 {
   int i,j;
   int *cornerpts[4];
 
-  for (j=0; j<dims[1]+direction; ++j) {
+  int nx = out->dimensions[0];
-    for (i=0; i<dims[0]+1-direction; ++i){
+  int ny = out->dimensions[1];
+  int nz = out->dimensions[2];
 
-      if (!checkmemeory(dims[2], ftab, neighbors, intersections)){
+  for (j=0; j<ny+direction; ++j) {
-	fprintf(stderr, "Could not allocat enough space\n");
+    for (i=0; i<nx+1-direction; ++i){
+
+      if (!checkmemeory(nz, out, intersections)){
+	fprintf(stderr, "Could not allocat enough space in process_vertical_faces\n");
 	exit(1);      
       }
-
       
       /* Vectors of point numbers */
-      int d[] = {2*dims[0], 2*dims[1], 2+2*dims[2]};
+      int d[] = {2*nx, 2*ny, 2+2*nz};
       igetvectors(d, 2*i+direction, 2*j+1-direction, plist, cornerpts);
       
       if(direction==1){
@@ -192,17 +185,21 @@ void process_vertical_faces(const int dims[3], int direction, sparse_table_t *ft
 	cornerpts[1] = tmp;
       }
 
-      int startface = ftab->position;
+      /* int startface = ftab->position; */
-      int num_intersections = *npoints - npillarpoints;
+      int startface = out->number_of_faces;
-      findconnections(2*dims[2]+2, cornerpts, npoints, 
+      /* int num_intersections = *npoints - npillarpoints; */
-		      *intersections+4*num_intersections, 
+      int num_intersections = out->number_of_nodes -
-		      *neighbors, work, ftab);
+	                      out->number_of_nodes_on_pillars;
       
-      int *ptr = *neighbors + 2*startface;
+      findconnections(2*nz+2, cornerpts, 
-      int len  = 2*ftab->position - 2*startface;
+		      *intersections+4*num_intersections, 
+		      work, out);
+      
+      int *ptr = out->face_neighbors + 2*startface;
+      int len = 2*out->number_of_faces - 2*startface;
 
-      compute_cell_index(dims, i-1+direction, j-direction, ptr,   len);
+      compute_cell_index(out->dimensions, i-1+direction, j-direction, ptr,   len);
-      compute_cell_index(dims, i,   j, ptr+1, len);
+      compute_cell_index(out->dimensions, i,   j, ptr+1, len);
 
     }
   }
@@ -225,34 +222,35 @@ static int linearindex(const int dims[3], int i, int j, int k)
   ACTNUM==0)
 
  */
-void process_horizontal_faces(const struct grdecl *g, 
+void process_horizontal_faces(int **intersections, 
-			      int                 *cell_index,
+			      int *plist,
-			      int                 *ncells,
+			      struct processed_grid *out)
-			      sparse_table_t      *faces, 
-			      int                 **neighbors, 
-			      int                 **intersections, 
-			      int                 *plist)
 {
   int i,j,k;
-  int *cell  = cell_index;
+
+  int nx = out->dimensions[0];
+  int ny = out->dimensions[1];
+  int nz = out->dimensions[2];
+
+  int *cell  = out->local_cell_index;
   int cellno = 0;
   
   /* dimensions of plist */
-  int  d[] = {2*g->dims[0], 2*g->dims[1], 2+2*g->dims[2]};
+  int  d[] = {2*nx, 2*ny, 2+2*nz};
 
 
-  for(j=0; j<g->dims[1]; ++j) {
+  for(j=0; j<ny; ++j) {
-    for (i=0; i<g->dims[0]; ++i) {
+    for (i=0; i<nx; ++i) {
 
-
+      
-      if (!checkmemeory(g->dims[2], faces, neighbors, intersections)){
+      if (!checkmemeory(nz, out, intersections)){
-	fprintf(stderr, "Could not allocat enough space\n");
+	fprintf(stderr, "Could not allocat enough space in process_horizontal_faces\n");
 	exit(1);
       }
 
       
-      int *f = (int*)faces->data + faces->ptr[faces->position];
+      int *f = out->face_nodes + out->face_ptr[out->number_of_faces];
-      int *n =       *neighbors  + 2*faces->position;
+      int *n =       out->face_neighbors  + 2*out->number_of_faces;
 
  
       /* Vectors of point numbers */
@@ -262,7 +260,7 @@ void process_horizontal_faces(const struct grdecl *g,
       int prevcell = -1;
 
 
-      for (k = 1; k<g->dims[2]*2+1; ++k){
+      for (k = 1; k<nz*2+1; ++k){
 
       	
 	/* Skip if space between face k and face k+1 is collapsed. */
@@ -272,7 +270,7 @@ void process_horizontal_faces(const struct grdecl *g,
 	    c[2][k] == c[2][k+1] && c[3][k] == c[3][k+1] &&
 	    k%2){
 	
-	  int idx = linearindex(g->dims, i,j,(k-1)/2);
+	  int idx = linearindex(out->dimensions, i,j,(k-1)/2);
 	  cell[idx] = -1;  
 
 	}
@@ -285,9 +283,9 @@ void process_horizontal_faces(const struct grdecl *g,
 	    *f++ = c[3][k];
 	    *f++ = c[2][k];
 
-	    faces->ptr[++faces->position] = f - (int*)faces->data;
+	    out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
-	    
+
-	    int thiscell = linearindex(g->dims, i,j,(k-1)/2);
+	    int thiscell = linearindex(out->dimensions, i,j,(k-1)/2);
 
 	    *n++ = prevcell;
 	    *n++ = prevcell = thiscell;
@@ -303,7 +301,7 @@ void process_horizontal_faces(const struct grdecl *g,
 	      *f++ = c[3][k];
 	      *f++ = c[2][k];
 
-	      faces->ptr[++faces->position] = f - (int*)faces->data;
+	      out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
 	      *n++ = prevcell;
 	      *n++ = prevcell = -1;
 	    }
@@ -312,7 +310,7 @@ void process_horizontal_faces(const struct grdecl *g,
       }
     }
   }
-  *ncells = cellno;
+  out->number_of_cells = cellno;
 }
 
 
@@ -350,50 +348,34 @@ static void approximate_intersection_pt(int *L, double *c, double *pt)
   Compute x,y and z coordinates for points on each pillar.
   Then, append x,y and z coordinates for extra points on faults.
  */
-static void compute_node_coordinates(const struct grdecl *g, 
+static void 
-				     double              *coordinates, 
+compute_intersection_coordinates(int                   *intersections, 
-				     int                 *intersections, 
+				 struct processed_grid *out)
-				     sparse_table_t      *pillarz, 
-				     int                 npillarpoints,
-				     int                 npoints)
 {
-  int i,k;
+  int n  = out->number_of_nodes;
-  int nx = g->dims[0];
+  int np = out->number_of_nodes_on_pillars;
-  int ny = g->dims[1];
-
-  double       *pt          = coordinates;
-  const double *c           = g->coord;
-
-  /* Loop over pillars */
-  int pillar = 0;
-  for (i=0; i< (nx+1)*(ny+1); ++i){
 
 
-    /* Loop over unique zcorn values - may be none */
+  /* Make sure the space allocated for nodes match the number of node. */
-    for (k=pillarz->ptr[pillar]; k<pillarz->ptr[pillar+1]; ++k) {
+  void *p = realloc (out->node_coordinates, 3*n*sizeof(double));
-	
+  if (p) {
-      /* Assign z-coordinate */
+    out->node_coordinates = p;
-      pt[2] = ((double*)pillarz->data)[k];
+  }
-
+  else {
-      /* Compute x- and y- coordinate */
+    fprintf(stderr, "Could not allocate extra space for intersections\n");
-      interpolate_pillar(c, pt);
-	
-      pt += 3;
-    }
-
-    ++pillar;
-    c += 6;
   }
 
   
   /* Append intersections */
-  int *itsct = intersections;
+  int    k;
-  for (k=npillarpoints; k<npoints; ++k){
+  double *pt    = out->node_coordinates + 3*np;
-
+  int    *itsct = intersections;
-    approximate_intersection_pt(itsct, coordinates, pt);
 
+  for (k=np; k<n; ++k){
+    approximate_intersection_pt(itsct, out->node_coordinates, pt);
     pt    += 3;
     itsct += 4;
+
   }
 }
 
@@ -448,114 +430,80 @@ void process_grdecl(const struct grdecl   *in,
   /* Code below assumes k index runs fastests, ie. that dimensions of
      table is permuted to (dims[2], dims[0], dims[1]) */
 
-  
-
-
-  
-
-
-  int *cell_index = calloc(nx*ny*nz,sizeof(*cell_index));
-  int ncells;
-  out->local_cell_index = cell_index;
-  /* out->number_of_cells = ncells */
-
-
-
-
-  /* ztab->data may need extra space temporarily due to simple boundary treatement  */
-  int            npillarpoints = 8*(nx+1)*(ny+1)*nz; 
-  int            npillars      = (nx+1)*(ny+1);
-  sparse_table_t *pillarz      = malloc_sparse_table(npillars, 
-						     npillarpoints, 
-						     sizeof(double));
-
-  /* Allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */
-  int    *plist = malloc( 4*nx*ny*(2*nz+2) * sizeof(int));
-
-  /* Fill plist of cornerpoint numbers and ztab of unique zcorn values. */
-  finduniquepoints(&g, plist, pillarz, tolerance);
-
-  npillarpoints = pillarz->ptr[npillars];
-
-  void *p = realloc_sparse_table (pillarz, npillars, npillarpoints, sizeof(double));
-  if (p) {
-    pillarz = p;
-  }else{
-    fprintf(stderr, "Could not reallocate space\n");
-    exit(1);
-  }
-
-
-
-  /* Process face geometry and cell-face topology on constant-i pillar pairs */
-
-
   const int BIGNUM = 64;
-  /* Unstructured storage of face nodes */
+  out->m                = BIGNUM/3;
-  sparse_table_t *faces = malloc_sparse_table(BIGNUM/3,
+  out->n                = BIGNUM;
-					     BIGNUM,
-					     sizeof(int));
-  out->face_ptr   = faces->ptr;
-  out->face_nodes = faces->data;
-  /* out->number_of_faces  faces->position; */
 
-  int *neighbors     = malloc(BIGNUM* sizeof(*neighbors));
+  out->face_neighbors   = malloc(BIGNUM* sizeof(*out->face_neighbors));
-  out->face_neighbors = neighbors;
+  out->face_nodes       = malloc(out->n * sizeof(int));
+  out->face_ptr         = malloc(out->m * sizeof(int));
+  out->face_ptr[0]      = 0;
 
+  out->dimensions[0]    = g.dims[0];
+  out->dimensions[1]    = g.dims[1];
+  out->dimensions[2]    = g.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(int));
+
+  
   /* internal */
   int *intersections = malloc(BIGNUM* sizeof(*intersections));
 
   /* internal */
   int *work          = malloc(2* (2*nz+2)*   sizeof(*work));
-  for(i=0; i<2* (2*nz+2); ++i){
+  for(i=0; i<4*(nz+1); ++i)   work[i] = -1;
-    work[i] = -1;
-  }
 
-  /* out->number_of_points */
+  /* Allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */
-  int npoints        = npillarpoints;
+  int    *plist = malloc( 4*nx*ny*(2*nz+2) * sizeof(int));
-  
-  faces->position = 0;
-  faces->ptr[0]   = 0;
 
-  /* faces with constant i */
-  process_vertical_faces(g.dims, 0, faces,
-			 &neighbors, &intersections, &npoints,
-			 npillarpoints, plist, work);
-  /* faces with constant j */
-  process_vertical_faces(g.dims, 1, faces,  
-			 &neighbors, &intersections, &npoints, 
-			 npillarpoints, plist, work);
-  
-  free(work);
-  
-  
 
-  process_horizontal_faces(&g, cell_index, &ncells, faces, &neighbors, 
+  /* Do actual work here:*/
-			   &intersections, plist);
+
+  finduniquepoints(&g, plist, tolerance, out);
+
+  free (zcorn);
+  free (actnum);
+
+  process_vertical_faces   (0, &intersections, plist, work, out);
+  process_vertical_faces   (1, &intersections, plist, work, out);  
+  process_horizontal_faces (&intersections, plist, out);
+
+  free (plist);
+  free (work);
+
+  compute_intersection_coordinates(intersections, out);  
+
+  free (intersections);
 
 
   /* Convert to local cell numbers in face_neighbors */
-  int *ptr=neighbors;;
+  int *ptr=out->face_neighbors;;
-  for (i=0; i<faces->position*2; ++i, ++ptr){
+  for (i=0; i<out->number_of_faces*2; ++i, ++ptr){
     if (*ptr != -1){
-      *ptr = cell_index[*ptr];
+      *ptr = out->local_cell_index[*ptr];
     }
   }
  
   /* Invert global-to-local map */
-  ptr = cell_index;
+  ptr = out->local_cell_index;
   for (i=0; i<nx*ny*nz; ++i){
-    if(cell_index[i]!=-1){
+    if(out->local_cell_index[i]!=-1){
       *ptr++ = i;
     }
   }
 
+
+
   /* Shrink memory allocated for cell_index */
-  if (ptr != cell_index){   /* always !*/
+  if (ptr != out->local_cell_index){   /* always !*/
-    p = realloc(cell_index, (ptr-cell_index)*sizeof(*cell_index));
+    void *p = realloc(out->local_cell_index, 
+		      (ptr-out->local_cell_index)*sizeof(*out->local_cell_index));
     if (p){
-      cell_index = p;
+      out->local_cell_index = p;
     }
     else{
       fprintf(stderr, "Could not reallocate space for index map\n");
@@ -563,29 +511,8 @@ void process_grdecl(const struct grdecl   *in,
     }
   }
 
-  /* compute node coordinates on pillars and new intersections */  
-    double       *coordinates = malloc(3*npoints * sizeof(*coordinates));
-  out->node_coordinates = coordinates;
-  compute_node_coordinates(&g, coordinates, intersections, pillarz, 
-			   npillarpoints, npoints);  
   
-  free_sparse_table(pillarz);
-  free (intersections);
-  free (plist);
 
-  free (zcorn);
-  free (actnum);
-
-
-  out->number_of_faces  = faces->position;
-  out->face_nodes       = faces->data;
-  out->face_ptr         = faces->ptr;
-
-  out->face_neighbors   = neighbors;
-  out->number_of_nodes  = npoints;
-  out->node_coordinates = coordinates;
-  out->number_of_cells  = ncells;
-  out->local_cell_index = cell_index;
 }
 
 /*-------------------------------------------------------*/

preprocess.h

@@ -47,12 +47,16 @@ struct grdecl{
 
 /* Output structure holding grid topology */
 struct processed_grid{
+  int m;
+  int n;
+  int    dimensions[3];     /* Cartesian dimension                                */
   int    number_of_faces;
   int    *face_nodes;       /* Nodes numbers of each face sequentially.           */
   int    *face_ptr;         /* Start position for each face in face_nodes.        */
   int    *face_neighbors;   /* Global cell numbers.  2 ints per face sequentially */
   
-  int    number_of_nodes;      
+  int    number_of_nodes; 
+  int    number_of_nodes_on_pillars;
   double *node_coordinates; /* 3 doubles per node, sequentially                   */
   
   int    number_of_cells;   /* number of active cells                             */

test.m

@@ -9,4 +9,4 @@ g = makeModel3([5,5,300]);
 g.ACTNUM=int32(g.ACTNUM);
 
 grdecl = readGRDECL(fullfile(ROOTDIR, 'examples','grids','GSmodel.grdecl'));
-grdecl = readGRDECL(fullfile(ROOTDIR, 'projects','co2','Sector15.grdecl'));
+%grdecl = readGRDECL(fullfile(ROOTDIR, 'projects','co2','Sector15.grdecl'));

uniquepoints.c

@@ -41,7 +41,7 @@ along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 
 
 #include "sparsetable.h"
-#include "grdecl.h"
+#include "preprocess.h"
 #include "uniquepoints.h"
 #include "matalloc.h"
 
@@ -206,6 +206,22 @@ static void dgetvectors(const int dims[3], int i, int j,
   v[3] = field + dims[2]*(ip + dims[0]* jp);
 }
 
+
+/*-----------------------------------------------------------------
+  Given a z coordinate, find x and y coordinates on line defined by
+  coord.  Coord points to a vector of 6 doubles [x0,y0,z0,x1,y1,z1].
+ */
+static void interpolate_pillar(const double *coord, double *pt)
+{
+  double a = (pt[2]-coord[2])/(coord[5]-coord[2]);
+  if (isinf(a) || isnan(a)){
+    a = 0;
+  }
+
+  pt[0]       = coord[0] + a*(coord[3]-coord[0]);
+  pt[1]       = coord[1] + a*(coord[4]-coord[1]);
+}
+
 /*-----------------------------------------------------------------
  Assign point numbers p such that "zlist(p)==zcorn".
  Assume that coordinate number is arranged in a
@@ -214,15 +230,31 @@ static void dgetvectors(const int dims[3], int i, int j,
 int finduniquepoints(const struct grdecl *g,
 		                            /* return values: */
 		     int           *plist, /* list of point numbers on each pillar*/
-		     sparse_table_t *ztab, 
+		     double tolerance, 
-		     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  */
+  int            npillarpoints = 8*(nx+1)*(ny+1)*nz; 
+  int            npillars      = (nx+1)*(ny+1);
+  sparse_table_t *ztab      = malloc_sparse_table(npillars, 
+						     npillarpoints, 
+						     sizeof(double));
+
+
+
+  int nc = g->dims[0]*g->dims[1]*g->dims[2];
+  out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates));
 
   double *zlist = ztab->data; /* casting void* to double* */
   int     *zptr = ztab->ptr;
 
-  int     i,j;
+  int     i,j,k;
 
   int     d1[3]  = {2*g->dims[0], 2*g->dims[1], 2*g->dims[2]};
   int     len    = 0;
@@ -231,6 +263,9 @@ int finduniquepoints(const struct grdecl *g,
 
   zptr[pos++] = zout - zlist;
 
+  double *coord = (double*)g->coord;
+  double *pt    = out->node_coordinates;
+
   /* Loop over pillars, find unique points on each pillar */
   for (j=0; j < g->dims[1]+1; ++j){
     for (i=0; i < g->dims[0]+1; ++i){
@@ -245,13 +280,22 @@ int finduniquepoints(const struct grdecl *g,
       len = createSortedList(     zout, d1[2], 4, z, a);
       len = uniquify        (len, zout, tolerance);      
 
+      /* Assign unique points */
+      for (k=0; k<len; ++k){
+	pt[2] = zout[k]; 
+	interpolate_pillar(coord, pt);
+	pt += 3;
+      }
+
       /* Increment pointer to sparse table of unique zcorn values */
       zout        = zout + len;
       zptr[pos++] = zout - zlist;
+
+      coord += 6;
     }
   }
-
+  out->number_of_nodes_on_pillars = zptr[pos-1];
-
+  out->number_of_nodes            = zptr[pos-1];
 
   /* Loop over all vertical sets of zcorn values, assign point numbers */
   int *p = plist;
@@ -279,6 +323,11 @@ int finduniquepoints(const struct grdecl *g,
       p += 2 + 2*g->dims[2];
     }
   }
+
+
+  free_sparse_table(ztab);
+
+
   return 1;
 }
 

uniquepoints.h

@@ -37,8 +37,8 @@ along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 
 int finduniquepoints(const struct grdecl *g,  /* input */
 		     int                 *p,  /* for each z0 in zcorn, z0 = z[p0] */
-		     sparse_table_t      *z,  /* list of uniq zcorn valules for each pillar*/
+		     double               t,  /* tolerance*/
-		     double               t); /* tolerance*/
+		     struct processed_grid *out);
 
 #endif // OPENRS_UNIQUEPOINTS_HEADER