From 45e3b4cede75b92b9b2d99d8c76677ce11ff6835 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Wed, 17 Jun 2009 13:08:59 +0000 Subject: [PATCH] First seemingly working version. Norner is processed. There are some differences between C and Matlab versions of the processed grid, though. --- facetopology.c | 92 ++++++++++++------- preprocess.c | 244 +++++++++++++++++++++++++------------------------ preprocess.h | 1 + processgrid.c | 132 +++++++++++++++++++++++++- sparsetable.c | 4 +- test.m | 11 ++- uniquepoints.c | 21 +++-- 7 files changed, 335 insertions(+), 170 deletions(-) diff --git a/facetopology.c b/facetopology.c index a894dd29..7c098b3e 100644 --- a/facetopology.c +++ b/facetopology.c @@ -3,6 +3,7 @@ #include #include #include +#include #include "sparsetable.h" #include "facetopology.h" @@ -32,7 +33,7 @@ static int *computeFaceTopology(int *a1, int intersect[4], int *faces) { - int mask[8]; + int mask[8] = {-1}; /* Which pillar points should we use? */ if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; } @@ -45,47 +46,54 @@ static int *computeFaceTopology(int *a1, /* but not all pillar points. This is encoded in mask. */ - mask[1] = intersect[0]; /* top-top */ - mask[3] = 0; - mask[5] = intersect[3]; /* bottom-bottom*/ - mask[7] = 0; + mask[1] = intersect[3]; /* top-top */ + mask[3] = -1; + mask[5] = intersect[0]; /* bottom-bottom*/ + mask[7] = -1; /* bottom-top */ - if (intersect[1]){ + if (intersect[1] != -1){ if(a1[0] > b1[1]){ /* intersection[1] left of (any) intersection[0] */ - mask[0] = 0; - mask[6] = 0; + mask[0] = -1; + mask[6] = -1; mask[7] = intersect[1]; } else{ - mask[2] = 0; - mask[4] = 0; + mask[2] = -1; + mask[4] = -1; mask[3] = intersect[1]; } } /* top-bottom */ - if (intersect[2]){ + if (intersect[2] != -1){ if(a1[1] < b1[0]){ /* intersection[2] left of (any) intersection[3] */ - mask[0] = 0; - mask[6] = 0; + mask[0] = -1; + mask[6] = -1; mask[7] = intersect[2]; } else{ - mask[2] = 0; - mask[4] = 0; + mask[2] = -1; + mask[4] = -1; mask[3] = intersect[2]; } } +/* fprintf(stderr, "intersect: %d %d %d %d\n", */ +/* intersect[0], */ +/* intersect[1], */ +/* intersect[2], */ +/* intersect[3]); */ +/* fprintf(stderr, "face: "); */ int k; int *f = faces; for (k=0; k<8; ++k){ - if(mask[k]){ + if(mask[k] != -1){ *f++ = mask[k]; +/* fprintf(stderr, "%d ", f[-1]); */ } } - +/* fprintf(stderr, "\n"); */ return f; } @@ -109,6 +117,10 @@ static int faceintersection(int *a1, int *a2, int *b1, int *b2) lineintersection(a1[0], a2[0], b1[0], b2[0]); } + +#define meaningful_face !((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], int *ptnumber, @@ -126,7 +138,6 @@ 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[*fpos]; */ int *f = (int*)ftab->data + ftab->ptr[ftab->position]; int *c = neighbors + 2*ftab->position; @@ -134,7 +145,8 @@ void findconnections(int n, int *pts[4], int k2 = 0; int i,j=0; - int intersect[4]; + int intersect[4]= {-1}; + for (i = 0; iptr[++ftab->position] = f - (int*)ftab->data; } @@ -220,14 +237,21 @@ void findconnections(int n, int *pts[4], j = j+1; } - /* Swap intersection records: top line of a is next bottom line of a */ + + + /* Swap intersection records: top line of a[i,i+1] is bottom line of a[i+1,i+2] */ int *tmp; tmp = itop; itop = ibottom; ibottom = tmp; + /* Zero out the "new" itop, set j to appropriate start position for next i */ + for(j=0;jmin(k1,k2);--j) itop[j-1]=0; + /* int k;for (k=0;k=dims[0] || j<0 || j >= dims[1]){ for(k=0; kposition; - int num_intersections = *npoints - npillarpoints; - + /* Vectors of point numbers */ int d[] = {2*dims[0], 2*dims[1], 2+2*dims[2]}; igetvectors(d, 2*i+direction, 2*j+1-direction, plist, cornerpts); - + if(direction==1){ /* swap */ - int *tmp = cornerpts[2]; + int *tmp = cornerpts[2]; cornerpts[2] = cornerpts[1]; cornerpts[1] = tmp; } + + int startface = ftab->position; + int num_intersections = *npoints - npillarpoints; findconnections(2*dims[2]+2, cornerpts, npoints, *intersections+4*num_intersections, *neighbors, work, ftab); - + int *ptr = *neighbors + 2*startface; int len = 2*ftab->position - 2*startface; + compute_cell_index(dims, i-1+direction, j-direction, ptr, len); compute_cell_index(dims, i, j, ptr+1, len); - } } } - - -/*-------------------------------------------------------*/ -void process_horizontal_faces(const struct grdecl *g, int *cell_index, - sparse_table_t *faces, int **neighbors, - int **intersections, int *plist) +int global_cell(const int dims[3], int i, int j, int k) +{ + return i+dims[0]*(j+dims[1]*k); +} +#include +/*-------------------------------------------------------*/ +void process_horizontal_faces(const struct grdecl *g, + int *cell_index, + int *ncells, + sparse_table_t *faces, + int **neighbors, + int **intersections, + int *plist) { - int nx = g->dims[0]; - int ny = g->dims[1]; - int i,j,k; - int *cell = cell_index; + int *cell = cell_index; int cellno = 0; /* compute constant-k faces */ @@ -187,54 +195,67 @@ void process_horizontal_faces(const struct grdecl *g, int *cell_index, int *f = (int*)faces->data; int *newface = (int*)faces->data + faces->ptr[faces->position]; - int *newneighbors = *neighbors + 2*faces->position; + int *newneighbors = *neighbors + 2*faces->position; - - + /* Vectors of point numbers */ int *c[4]; - int d[] = {2*g->dims[0], 2*g->dims[1], 2+2*g->dims[2]}; + int d[] = {2*g->dims[0], 2*g->dims[1], 2+2*g->dims[2]}; igetvectors(d, 2*i+1, 2*j+1, plist, c); - int kprev = -1; + int previous_global_cell = -1; for (k = 1; kdims[2]*2+1; ++k){ + /* Skip if space between face k and face k+1 is collapsed. */ - /* Note that inactive cells (with ACNUM==0) have all been */ - /* collased in finduniquepoints. */ - if (c[0][k] == c[0][k+1] && - c[1][k] == c[1][k+1] && - c[2][k] == c[2][k+1] && - c[3][k] == c[3][k+1] && - k%2){ - *cell++ = -1; - kprev = -1; + /* Note that inactive cells (with ACTNUM==0) have all been */ + /* collapsed in finduniquepoints. */ + if (c[0][k] == c[0][k+1] && c[1][k] == c[1][k+1] && + c[2][k] == c[2][k+1] && c[3][k] == c[3][k+1]){ - }else{ - - /* Add face ### - room for some refinement here #### */ - *newface++ = c[0][k]; - *newface++ = c[1][k]; - *newface++ = c[2][k]; - *newface++ = c[3][k]; + if (k%2) cell[global_cell(g->dims, i,j,(k-1)/2)] = -1; - faces->ptr[++faces->position] = newface - f; - *newneighbors++ = kprev != -1 ? i+nx*(j+ny*(kprev-1)/2) : -1; - *newneighbors++ = k != g->dims[2]*2 ? i+nx*(j+ny*(k -1)/2) : -1; - kprev = k; + } + else{ - if (k%2)*cell++ = cellno++; + if (k%2){ + /* Add face */ + *newface++ = c[0][k]; + *newface++ = c[1][k]; + *newface++ = c[3][k]; + *newface++ = c[2][k]; + + faces->ptr[++faces->position] = newface - f; + *newneighbors++ = previous_global_cell; + *newneighbors++ = previous_global_cell = global_cell(g->dims, i,j,(k-1)/2); + + cell[global_cell(g->dims, i,j,(k-1)/2)] = cellno++; + + } + else{ + if (previous_global_cell != -1){ + /* Add face */ + *newface++ = c[0][k]; + *newface++ = c[1][k]; + *newface++ = c[3][k]; + *newface++ = c[2][k]; + + faces->ptr[++faces->position] = newface - f; + *newneighbors++ = previous_global_cell; + *newneighbors++ = previous_global_cell = -1; + } + } } } } } - + *ncells = cellno; } -static void approximate_intersection_pt(int *intersection, - double *coordinates, double *pt) +static void approximate_intersection_pt(int *intersection, + double *coordinates, double *pt) { double z0 = coordinates[3*intersection[0]+2]; double z1 = coordinates[3*intersection[1]+2]; @@ -253,11 +274,11 @@ static void approximate_intersection_pt(int *intersection, } static void compute_node_coordinates(const struct grdecl *g, - double *coordinates, - int *intersections, - sparse_table_t *pillarz, - int npillarpoints, - int npoints) + double *coordinates, + int *intersections, + sparse_table_t *pillarz, + int npillarpoints, + int npoints) { int i,k; int nx = g->dims[0]; @@ -269,9 +290,6 @@ static void compute_node_coordinates(const struct grdecl *g, /* Loop over pillars */ int pillar = 0; for (i=0; i< (nx+1)*(ny+1); ++i){ - - fprintf(stderr, "pillar: %f %f %f - %f %f %f\n", - c[0], c[1], c[2], c[3],c[4],c[5]); /* Loop over unique zcorn values - may be none */ @@ -283,7 +301,6 @@ static void compute_node_coordinates(const struct grdecl *g, /* Compute x- and y- coordinate */ interpolate_pillar(c, pt); - fprintf(stderr, "pt : %f %f %f\n", pt[0], pt[1], pt[2]); pt += 3; } @@ -295,10 +312,10 @@ static void compute_node_coordinates(const struct grdecl *g, /* Append intersections */ int *itsct = intersections; for (k=npillarpoints; kptr[npillars]; + + void *p = realloc_sparse_table (pillarz, npillars, npillarpoints, sizeof(double)); if (p) { pillarz = p; @@ -341,7 +360,7 @@ void processGrdecl(const struct grdecl *g, double tol, struct processed_grid *ou } - fprintf(stderr, "process face geometry\n"); + /* Process face geometry and cell-face topology on constant-i pillar pairs */ @@ -355,84 +374,70 @@ void processGrdecl(const struct grdecl *g, double tol, struct processed_grid *ou int *intersections = malloc(BIGNUM* sizeof(*intersections)); - int *work = calloc(2* (2*nz+2), sizeof(*work)); + int *work = malloc(2* (2*nz+2)* sizeof(*work)); + for(i=0; i<2* (2*nz+2); ++i){ + work[i] = -1; + } int npoints = npillarpoints; faces->position = 0; faces->ptr[0] = 0; /* faces with constant i */ - process_vertical_faces(g->dims, 0, faces, - &neighbors, &intersections, &npoints, + 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); - fprintf(stderr, "processing k-constant faces and cell numbering \n"); - fprintf(stderr, "Cell numbering is warped due to k running fastest\n"); - int *cell_index = malloc(nx*ny*nz*sizeof(*cell_index)); - process_horizontal_faces(g, cell_index, faces, &neighbors, &intersections, plist); + int *cell_index = calloc(nx*ny*nz,sizeof(*cell_index)); + int ncells; + + process_horizontal_faces(g, cell_index, &ncells, faces, &neighbors, + &intersections, plist); - - /* */ - /* */ - /* */ - /* D E B U G CODE */ - /* */ - /* */ - /* */ - -#if 1 - fprintf(stderr, "\nfaces\nnumfaces %d\n", faces->position); - for (i=0; iposition; ++i){ - for (k=faces->ptr[i]; kptr[i+1]; ++k){ - fprintf(stderr, "%d ", ((int*)faces->data)[k]); + /* Convert to local cell numbers in face_neighbors */ + int *ptr=neighbors;; + for (i=0; iposition*2; ++i, ++ptr){ + if (*ptr != -1){ + *ptr = cell_index[*ptr]; } - fprintf(stderr, "\n"); } - - - fprintf(stderr, "\nneighbors\n"); - int *iptr = neighbors; - for(k=0; kposition; ++k){ - fprintf(stderr, " (%d %d)\n", iptr[0], iptr[1]); - ++iptr; - ++iptr; - } - - fprintf(stderr, "\nline intersections:\n"); - iptr = intersections; - int numintersections = npoints - npillarpoints; - for(k=0; kface_nodes = faces->data; out->face_ptr = faces->ptr; out->face_neighbors = neighbors; - out->number_of_nodes = 0; + out->number_of_nodes = npoints; out->node_coordinates = coordinates; - out->number_of_cells = 0; + out->number_of_cells = ncells; out->local_cell_index = cell_index; - } diff --git a/preprocess.h b/preprocess.h index 18aa8430..3f0f7ee3 100644 --- a/preprocess.h +++ b/preprocess.h @@ -22,6 +22,7 @@ struct processed_grid{ double *node_coordinates; /* 3 doubles per node, sequentially */ int number_of_cells; /* number of active cells */ + int number_of_active_cells; int *local_cell_index; /* Global to local map */ }; diff --git a/processgrid.c b/processgrid.c index 2f4199f3..55a62674 100644 --- a/processgrid.c +++ b/processgrid.c @@ -7,6 +7,133 @@ #include "preprocess.h" #include "mxgrdecl.h" + +void fill_grid(mxArray **out, struct processed_grid *grid) +{ + const char *names[] = {"nodes", "faces", "cells", "cellFaces", "faceNodes"}; + mxArray *G = mxCreateStructMatrix(1,1,5,names); + + int i,j; + double *ptr; + + + /* nodes */ + const char *n2[] = {"num", "coords"}; + mxArray *nodes = mxCreateStructMatrix(1,1,2,n2); + mxSetField(nodes, 0, "num", mxCreateDoubleScalar(grid->number_of_nodes)); + + mxArray *nodecoords = mxCreateDoubleMatrix(grid->number_of_nodes, 3, mxREAL); + ptr = mxGetPr(nodecoords); + + for (j=0;j<3;++j){ + for(i=0; inumber_of_nodes; ++i){ + ptr[i+grid->number_of_nodes*j] = grid->node_coordinates[3*i+j]; + } + } + mxSetField(nodes, 0, "coords", nodecoords); + mxSetField(G, 0, "nodes", nodes); + + + /* faces */ + const char *n3[] = {"num", "neighbors", "numNodes"}; + mxArray *faces = mxCreateStructMatrix(1,1,3,n3); + + + mxSetField(faces, 0, "num", mxCreateDoubleScalar(grid->number_of_faces)); + + mxArray *faceneighbors = mxCreateDoubleMatrix(grid->number_of_faces, 2, mxREAL); + /* int i,j; */ + ptr = mxGetPr(faceneighbors); + for(j=0; j<2; ++j){ + for (i=0; inumber_of_faces; ++i){ + int ix = grid->face_neighbors[2*i+j]; + if (ix == -1){ + ptr[i+grid->number_of_faces*j] = 0; + }else{ + ptr[i+grid->number_of_faces*j] = ix+1;/* grid->local_cell_index[ix]+1; */ + } + } + } + mxSetField(faces, 0, "neighbors", faceneighbors); + + mxArray *numnodes = mxCreateDoubleMatrix(grid->number_of_faces, 1, mxREAL); + ptr = mxGetPr(numnodes); + for (i=0; inumber_of_faces; ++i){ + ptr[i] = grid->face_ptr[i+1]-grid->face_ptr[i]; + } + mxSetField(faces, 0, "numNodes", numnodes); + + mxSetField(G, 0, "faces", faces); + + + const char *n4[] = {"num", "numFaces", "map"}; + mxArray *cells = mxCreateStructMatrix(1,1,3,n4); + + mxSetField(cells, 0, "num", mxCreateDoubleScalar(grid->number_of_cells)); + mxArray *map = mxCreateDoubleMatrix(grid->number_of_cells, 1, mxREAL); + ptr = mxGetPr(map); + for(i=0; inumber_of_cells; ++i){ + ptr[i] = grid->local_cell_index[i]+1; + } + mxSetField(cells, 0, "map", map); + + + mxArray *numfaces = mxCreateDoubleMatrix(grid->number_of_cells, 1, mxREAL); + + ptr = mxGetPr(numfaces); + for(i=0; inumber_of_cells; ++i){ + ptr[i] = 0.0; + } + for (i=0; i<2*grid->number_of_faces; ++i){ + int c=grid->face_neighbors[i]; + if(c != -1) { + ptr[c]++; + } + } + mxSetField(cells, 0, "numFaces", numfaces); + + mxSetField(G, 0, "cells", cells); + + + + int *counter = calloc(grid->number_of_cells, sizeof(*counter)); + int num_half_faces = 0; + for(i=0; inumber_of_cells; ++i){ + counter[i] = num_half_faces; + num_half_faces += ptr[i]; + } + + mxArray *cellfaces = mxCreateDoubleMatrix(num_half_faces, 1, mxREAL); + ptr = mxGetPr(cellfaces); + + + for (i=0; inumber_of_faces; ++i){ + int c1 = grid->face_neighbors[2*i]; + int c2 = grid->face_neighbors[2*i+1]; + if(c1 != -1) ptr[counter[c1]++] = i+1; + if(c2 != -1) ptr[counter[c2]++] = i+1; + } + + mxSetField(G, 0, "cellFaces", cellfaces); + + + int n = grid->face_ptr[grid->number_of_faces]; + + mxArray *facenodes = mxCreateDoubleMatrix(n, 1, mxREAL); + ptr = mxGetPr(facenodes); + for (i=0; iface_nodes[i]+1; + } + + mxSetField(G, 0, "faceNodes", facenodes); + free(counter); + + + out[0] = G; + +} + + /* Gateway routine for Matlab mex function. */ /*-------------------------------------------------------*/ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) @@ -19,8 +146,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) processGrdecl(&g, 0.0, &out); - if (plhs == 0){ - ;/* write to matlab struct */ + if (plhs >0){ + /* write to matlab struct */ + fill_grid(plhs, &out); } /* Free whatever was allocated in initGrdecl. */ diff --git a/sparsetable.c b/sparsetable.c index 1fd1260b..e88f00c1 100644 --- a/sparsetable.c +++ b/sparsetable.c @@ -28,7 +28,7 @@ sparse_table_t *malloc_sparse_table(int m, int n, int datasz) if(!(tab->data = malloc(n * datasz))){ - fprintf(stderr, "Could not allocate space for sparse data\n"); + fprintf(stderr, "Could not allocate space for sparse data(%d)\n", n); free_sparse_table(tab); return NULL; } @@ -51,7 +51,7 @@ sparse_table_t *realloc_sparse_table(sparse_table_t *tab, int m, int n, int data if(p){ tab->data = p; }else{ - fprintf(stderr, "Could not reallocate space for sparse data\n"); + fprintf(stderr, "Could not reallocate space for sparse data(%d)\n", n); free_sparse_table(tab); return NULL; } diff --git a/test.m b/test.m index c4922761..afa12e00 100644 --- a/test.m +++ b/test.m @@ -1,6 +1,11 @@ - -g=simpleGrdecl([4, 2, 4], @(x) -0.055+0.11*x+0.011 ); -G=processGRDECL(g); +nx = 5; +ny = 7; +nz = 11; +%g=simpleGrdecl([nx, ny, nz], @(x) -0.055+0.11*x+0.011 ); +g = makeModel3([200,220,30]); +%G=processGRDECL(g); %clf,plotGrid(G);view(3); g.ACTNUM=int32(g.ACTNUM); + +grdecl = readGRDECL(fullfile(ROOTDIR, 'examples','grids','GSmodel.grdecl')); diff --git a/uniquepoints.c b/uniquepoints.c index 394b98db..e4fec259 100644 --- a/uniquepoints.c +++ b/uniquepoints.c @@ -47,9 +47,11 @@ static int createSortedList(double *list, int n, int m, /*-------------------------------------------------------*/ static int uniquify(int n, double *list, double tolerance) { + if (n<1) return 0; int i; - int pos = 1; /* Keep first value */ - double val = list[pos]; + int pos = 0; + double val = list[pos++];/* Keep first value */ + for (i=1; i tolerance){ val = list[i]; @@ -61,7 +63,7 @@ static int uniquify(int n, double *list, double tolerance) if (list[n-1] - val > tolerance){ list[pos-1] = list[n-1]; } - + return pos; } @@ -112,6 +114,7 @@ static int* assignPointNumbers(int begin, } *p++ = INT_MAX;/* Padding to ease processing of faults */ + return p; } @@ -169,8 +172,8 @@ void finduniquepoints(const struct grdecl *g, double *zout = zlist; int pos = 0; - zptr[0] = 0; - + zptr[pos++] = zout - zlist; + /* 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){ @@ -184,10 +187,10 @@ void finduniquepoints(const struct grdecl *g, len = createSortedList( zout, d1[2], 4, z, a); len = uniquify (len, zout, 0.0); - + /* Increment pointer to sparse table of unique zcorn values */ zout = zout + len; - zptr[++pos] = zout - zlist; + zptr[pos++] = zout - zlist; } } @@ -197,7 +200,7 @@ void finduniquepoints(const struct grdecl *g, int *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); @@ -212,8 +215,8 @@ void finduniquepoints(const struct grdecl *g, assignPointNumbers(zptr[pix], zptr[pix+1], zlist, 2*g->dims[2], z, a, p, 0.0); - p += 2 + 2*g->dims[2]; + p += 2 + 2*g->dims[2]; } } }