First seemingly working version. Norner is processed. There are some
differences between C and Matlab versions of the processed grid, though.
This commit is contained in:
Jostein R. Natvig
parent
a4b48c494b
commit
45e3b4cede
@ -3,6 +3,7 @@
|
||||
#include <string.h>
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
#include <limits.h>
|
||||
|
||||
#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; i<n-1; ++i){
|
||||
if (a1[i] == a1[i+1] && a2[i] == a2[i+1]) continue;
|
||||
@ -142,9 +154,13 @@ void findconnections(int n, int *pts[4],
|
||||
|
||||
|
||||
|
||||
|
||||
while(j<n-1 && (b1[j] < a1[i+1] || b2[j] < a2[i+1])){
|
||||
|
||||
/*debug:*/
|
||||
/* int k; */
|
||||
/* for (k=0; k<n; ++k) itop[k]=-1; */
|
||||
/* for (k=0; k<n; ++k) ibottom[k]=-1; */
|
||||
|
||||
if (b1[j] == b1[j+1] && b2[j] == b2[j+1]){
|
||||
itop[j+1] = itop[j];
|
||||
++j;
|
||||
@ -158,20 +174,17 @@ void findconnections(int n, int *pts[4],
|
||||
if (faceintersection(a1+i, a2+i, b1+j, b2+j)){
|
||||
|
||||
|
||||
/* Add neighbors to list of neighbors if not any first or */
|
||||
/* last points are involved in face geometry. */
|
||||
if (!((i==0) && (j==0)) && !((i==n-2) && (j==n-2))){
|
||||
*c++ = i%2 ? (i-1)/2 : -1;
|
||||
*c++ = j%2 ? (j-1)/2 : -1;
|
||||
}
|
||||
|
||||
|
||||
/* Completely matching faces */
|
||||
if (a1[i]==b1[j] && a1[i+1]==b1[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 (!((i==0) && (j==0)) && !((i==n-2) && (j==n-2))){
|
||||
if (meaningful_face){
|
||||
/* Add neighbors cells */
|
||||
*c++ = i%2 ? (i-1)/2 : -1;
|
||||
*c++ = j%2 ? (j-1)/2 : -1;
|
||||
|
||||
/* face */
|
||||
*f++ = a1[i];
|
||||
*f++ = a1[i+1];
|
||||
*f++ = a2[i+1];
|
||||
@ -195,7 +208,7 @@ void findconnections(int n, int *pts[4],
|
||||
|
||||
|
||||
}else{
|
||||
itop[j+1] = 0;
|
||||
itop[j+1] = -1;
|
||||
}
|
||||
|
||||
/* Update intersection record */
|
||||
@ -206,7 +219,11 @@ void findconnections(int n, int *pts[4],
|
||||
|
||||
|
||||
/* Add face to list of faces if not any first or last points are involved. */
|
||||
if (!((i==0) && (j==0)) && !((i==n-2) && (j==n-2))){
|
||||
if (meaningful_face){
|
||||
/* neighbor cells */
|
||||
*c++ = i%2 ? (i-1)/2 : -1;
|
||||
*c++ = j%2 ? (j-1)/2 : -1;
|
||||
|
||||
f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f);
|
||||
ftab->ptr[++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;j<n; ++j) itop[j]=-1;
|
||||
|
||||
/* Zero out itop, set j to appropriate start position for next i */
|
||||
for(;j>min(k1,k2);--j) itop[j-1]=0;
|
||||
/* int k;for (k=0;k<n; ++k){ */
|
||||
/* if (itop [k] != -1) fprintf(stderr, "itop not -1\n"); */
|
||||
/* if (ibottom[k] != -1) fprintf(stderr, "ibottom not -1\n"); */
|
||||
/* } */
|
||||
|
||||
/* Now, j = min(k1, k2) and itop is all zeros */
|
||||
j = min(k1, k2);
|
||||
}
|
||||
}
|
||||
|
||||
244
preprocess.c
244
preprocess.c
@ -58,6 +58,7 @@ void free_processed_grid(struct processed_grid *g)
|
||||
/*-------------------------------------------------------*/
|
||||
void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len)
|
||||
{
|
||||
|
||||
int k;
|
||||
if (i<0 || i>=dims[0] || j<0 || j >= dims[1]){
|
||||
for(k=0; k<len; k+=2){
|
||||
@ -66,7 +67,8 @@ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len
|
||||
}else{
|
||||
for(k=0; k<len; k+=2){
|
||||
if (neighbors[k] != -1){
|
||||
neighbors[k] = i + dims[0]*(j + dims[1]*neighbors[k]);
|
||||
int tmp = i + dims[0]*(j + dims[1]*neighbors[k]);
|
||||
neighbors[k] = tmp;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -122,7 +124,8 @@ void process_vertical_faces(const int dims[3], int direction, sparse_table_t *ft
|
||||
|
||||
int i,j;
|
||||
int *cornerpts[4];
|
||||
/* constant i-faces */
|
||||
|
||||
/* constant i- or j-faces */
|
||||
for (j=0; j<dims[1]+direction; ++j) {
|
||||
for (i=0; i<dims[0]+1-direction; ++i){
|
||||
|
||||
@ -131,47 +134,52 @@ void process_vertical_faces(const int dims[3], int direction, sparse_table_t *ft
|
||||
exit(1);
|
||||
}
|
||||
|
||||
int startface = ftab->position;
|
||||
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 <limits.h>
|
||||
/*-------------------------------------------------------*/
|
||||
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; k<g->dims[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; k<npoints; ++k){
|
||||
/* Approximate intersection */
|
||||
|
||||
approximate_intersection_pt(itsct, coordinates, pt);
|
||||
fprintf(stderr, "%f %f %f\n", pt[0], pt[1], pt[2]);
|
||||
pt += 3;
|
||||
|
||||
pt += 3;
|
||||
itsct += 4;
|
||||
}
|
||||
}
|
||||
@ -308,7 +325,7 @@ static void compute_node_coordinates(const struct grdecl *g,
|
||||
/*-------------------------------*/
|
||||
void processGrdecl(const struct grdecl *g, double tol, struct processed_grid *out)
|
||||
{
|
||||
int i,k;
|
||||
int i;
|
||||
|
||||
/* Code below assumes k index runs fastests, ie. that dimensions of
|
||||
table is permuted to (dims[2], dims[0], dims[1]) */
|
||||
@ -332,6 +349,8 @@ void processGrdecl(const struct grdecl *g, double tol, struct processed_grid *ou
|
||||
finduniquepoints(g, plist, pillarz);
|
||||
|
||||
npillarpoints = pillarz->ptr[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; i<faces->position; ++i){
|
||||
for (k=faces->ptr[i]; k<faces->ptr[i+1]; ++k){
|
||||
fprintf(stderr, "%d ", ((int*)faces->data)[k]);
|
||||
/* Convert to local cell numbers in face_neighbors */
|
||||
int *ptr=neighbors;;
|
||||
for (i=0; i<faces->position*2; ++i, ++ptr){
|
||||
if (*ptr != -1){
|
||||
*ptr = cell_index[*ptr];
|
||||
}
|
||||
fprintf(stderr, "\n");
|
||||
}
|
||||
|
||||
|
||||
fprintf(stderr, "\nneighbors\n");
|
||||
int *iptr = neighbors;
|
||||
for(k=0; k<faces->position; ++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; k<numintersections; ++k){
|
||||
fprintf(stderr, " (%d %d %d %d)\n", iptr[0], iptr[1], iptr[2], iptr[3]);
|
||||
iptr+=4;
|
||||
}
|
||||
|
||||
|
||||
/* fprintf(stderr, "number of cells : %d\n", nx*ny*nz); */
|
||||
/* for (k=0; k<nx*ny*nz; ++k){ */
|
||||
/* fprintf(stderr, "cell index %d is %d\n", k, cell_index[k]); */
|
||||
/* } */
|
||||
#endif
|
||||
|
||||
/* compute node coordinates on pillars and new intersections */
|
||||
|
||||
/* Invert global-to-local map */
|
||||
ptr = cell_index;
|
||||
for (i=0; i<nx*ny*nz; ++i){
|
||||
if(cell_index[i]!=-1){
|
||||
*ptr++ = i;
|
||||
}
|
||||
}
|
||||
|
||||
/* Shrink memory allocated for cell_index */
|
||||
p = realloc(cell_index, (ptr-cell_index)*sizeof(*cell_index));
|
||||
if (p){
|
||||
cell_index = p;
|
||||
}
|
||||
else{
|
||||
fprintf(stderr, "Could not reallocate space for index map\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
|
||||
/* compute node coordinates on pillars and new intersections */
|
||||
double *coordinates = malloc(3*npoints * sizeof(*coordinates));
|
||||
|
||||
compute_node_coordinates(g, coordinates, intersections, pillarz,
|
||||
npillarpoints, npoints);
|
||||
compute_node_coordinates(g, coordinates, intersections, pillarz,
|
||||
npillarpoints, npoints);
|
||||
|
||||
|
||||
free_sparse_table(pillarz);
|
||||
|
||||
|
||||
free_sparse_table(pillarz);
|
||||
|
||||
/* compute intersections */
|
||||
free (intersections);
|
||||
|
||||
free (plist);
|
||||
@ -442,9 +447,8 @@ void processGrdecl(const struct grdecl *g, double tol, struct processed_grid *ou
|
||||
out->face_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;
|
||||
|
||||
}
|
||||
|
||||
@ -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 */
|
||||
};
|
||||
|
||||
|
||||
132
processgrid.c
132
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; i<grid->number_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; i<grid->number_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; i<grid->number_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; i<grid->number_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; i<grid->number_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; i<grid->number_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; i<grid->number_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; i<n; ++i) {
|
||||
ptr[i] = grid->face_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. */
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
11
test.m
11
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'));
|
||||
|
||||
@ -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<n; ++i){
|
||||
if (list[i] - val > 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];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user