From 2024d80518404bd85d886a7cfdac8d69d9b76609 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Fri, 19 Jun 2009 09:17:44 +0000 Subject: [PATCH] Some comments. --- preprocess.c | 121 ++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 85 insertions(+), 36 deletions(-) diff --git a/preprocess.c b/preprocess.c index e820d85f..d7c77924 100644 --- a/preprocess.c +++ b/preprocess.c @@ -43,27 +43,34 @@ along with OpenRS. If not, see . #include "uniquepoints.h" #include "facetopology.h" -/* No checking of input arguments in this code! */ + #define min(i,j) ((i)<(j) ? (i) : (j)) #define max(i,j) ((i)>(j) ? (i) : (j)) -static void interpolate_pillar(const double *coord, double *pt) +/*----------------------------------------------------------------- + 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]. + */ +istatic void interpolate_pillar(const double *coord, double *pt) { - double a; - if (fabs(coord[5]-coord[2]) < DBL_EPSILON){ + double a = (pt[2]-coord[2])/(coord[5]-coord[2]); + if (isinf(a) || isnan(a)){ a = 0; } - else{ - a = (pt[2]-coord[2])/(coord[5]-coord[2]); - } + pt[0] = coord[0] + a*(coord[3]-coord[0]); pt[1] = coord[1] + a*(coord[4]-coord[1]); } -/*-------------------------------------------------------*/ +/*----------------------------------------------------------------- + Given a vector with k index running faster than i running + faster than j, and Cartesian dimensions , find pointers to the + (i-1, j-1, 0), (i-1, j, 0), (i, j-1, 0) and (i, j, 0) elements of + field. + */ static void igetvectors(int dims[3], int i, int j, int *field, int *v[]) { @@ -79,23 +86,17 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[]) } -/*-------------------------------------------------------*/ -void free_processed_grid(struct processed_grid *g) -{ - if( g ){ - free ( g->face_nodes ); - free ( g->face_ptr ); - free ( g->face_neighbors ); - free ( g->node_coordinates ); - free ( g->local_cell_index ); - } -} +/*----------------------------------------------------------------- + Special purpose -/*-------------------------------------------------------*/ + Convert from k-index to Cartesian index i+nx*(j+ny*k) for every other + element in neighbors. + + */ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len) { @@ -115,8 +116,9 @@ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len } -/* Ensure there's sufficient memory */ -/*-------------------------------------------------------*/ +/*----------------------------------------------------------------- + Ensure there's sufficient memory + */ int checkmemeory(int nz, sparse_table_t *ftab, int **neighbors, int **intersections) { @@ -154,18 +156,22 @@ int checkmemeory(int nz, sparse_table_t *ftab, int **neighbors, int **intersecti return 1; } -/*-------------------------------------------------------*/ +/*----------------------------------------------------------------- + For each vertical face (i.e. i or j constant), + -find point numbers for the corners and + -cell neighbors. + -new points on faults defined by two intgersecting lines. + + direction == 0 : constant-i faces. + direction == 1 : constant-j faces. + */ void process_vertical_faces(const int dims[3], int direction, sparse_table_t *ftab, int **neighbors, int **intersections, int *npoints, int npillarpoints, int *plist, int *work) { - /* direction == 0 : constant i-faces - direction == 1 : constant j-faces */ - int i,j; int *cornerpts[4]; - /* constant i- or j-faces */ for (j=0; j. Exclude + cells that are have collapsed coordinates. (This includes cells with + ACTNUM==0) + + */ void process_horizontal_faces(const struct grdecl *g, int *cell_index, int *ncells, @@ -252,10 +269,11 @@ void process_horizontal_faces(const struct grdecl *g, /* 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]){ - - if (k%2) cell[linearindex(g->dims, i,j,(k-1)/2)] = -1; - + 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); + cell[idx] = -1; } else{ @@ -298,15 +316,29 @@ void process_horizontal_faces(const struct grdecl *g, } +/*----------------------------------------------------------------- + On input, + L points to 4 ints that indirectly refers to points in c. + c points to array of coordinates [x0,y0,z0,x1,y1,z1,...,xn,yn,zn]. + pt points to array of 3 doubles. + + On output, + pt holds coordinates to intersection between lines given by point + numbers L[0]-L[1] and L[2]-L[3]. +*/ static void approximate_intersection_pt(int *L, double *c, double *pt) { + double z0 = c[3*L[0]+2]; double z1 = c[3*L[1]+2]; double z2 = c[3*L[2]+2]; double z3 = c[3*L[3]+2]; double a = (z2-z0)/(z1-z0 - (z3-z2)); - assert(z1-z0 - (z3-z2)!=0); + if (isinf(a) || isnan(a)){ + a = 0; + } + pt[0] = c[3*L[0]+0]* (1.0-a) + c[3*L[1]+0]* a; pt[1] = c[3*L[0]+1]* (1.0-a) + c[3*L[1]+1]* a; @@ -314,7 +346,10 @@ 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, double *coordinates, int *intersections, @@ -363,8 +398,9 @@ static void compute_node_coordinates(const struct grdecl *g, } -/* Gateway routine. */ -/*-------------------------------*/ +/*----------------------------------------------------------------- + Public interface + */ void process_grdecl(const struct grdecl *g, double tolerance, struct processed_grid *out) @@ -514,3 +550,16 @@ void process_grdecl(const struct grdecl *g, out->number_of_cells = ncells; out->local_cell_index = cell_index; } + +/*-------------------------------------------------------*/ +void free_processed_grid(struct processed_grid *g) +{ + if( g ){ + free ( g->face_nodes ); + free ( g->face_ptr ); + free ( g->face_neighbors ); + free ( g->node_coordinates ); + free ( g->local_cell_index ); + } +} +