OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Moved permutation of ACTNUM and ZCORN to process_grdecl. The input is

Browse Source 
now unaltered on exit.
This commit is contained in:
Jostein R. Natvig
2009-06-19 11:31:16 +00:00
parent 12c6e7170c
commit 7cc7596fc8
5 changed files with 89 additions and 71 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
mxgrdecl.c
View File

@@ -47,7 +47,7 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
/*-------------------------------------------------------*/
void mx_init_grdecl(struct grdecl *g, const mxArray *s)
{
int i,j,k,n;
int i,n;
mxArray *field;
int numel;
@@ -71,20 +71,9 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s)
numel != g->dims[0]*g->dims[1]*g->dims[2] ){
mexErrMsgTxt("ACTNUM field must be nx*ny*nz numbers int32");
}
/* grdecl.actnum = permute(actnum, [3,1,2]); */
int *actnum = mxGetData(field);
int *a = malloc(n* sizeof(*g->actnum));
int *iptr = a;
for (j=0; j<g->dims[1]; ++j){
for (i=0; i<g->dims[0]; ++i){
for (k=0; k<g->dims[2]; ++k){
*iptr++ = actnum[i+g->dims[0]*(j+g->dims[1]*k)];
}
}
}
g->actnum = a;
g->actnum = mxGetData(field);
field = mxGetField(s, 0, "COORD");
numel = mxGetNumberOfElements(field);
@@ -101,28 +90,5 @@ void mx_init_grdecl(struct grdecl *g, const mxArray *s)
numel != 8*g->dims[0]*g->dims[1]*g->dims[2]){
mexErrMsgTxt("ZCORN field must have 8*nx*ny*nz doubles.");
}
double *zcorn = mxGetPr(field);
/* grdecl.zcorn = permute(zcorn, [3,1,2]); */
double *z = malloc(n*8*sizeof(*g->zcorn));
double *dptr = z;
for (j=0; j<2*g->dims[1]; ++j){
for (i=0; i<2*g->dims[0]; ++i){
for (k=0; k<2*g->dims[2]; ++k){
*dptr++ = zcorn[i+2*g->dims[0]*(j+2*g->dims[1]*k)];
}
}
}
g->zcorn = z;
g->zcorn = mxGetPr(field);
}
/* Free stuff that was allocated in initgrdecl. */
/*-------------------------------------------------------*/
void free_grdecl(struct grdecl *g)
{
free((double*)g->zcorn); g->zcorn = NULL;
free((double*)g->actnum); g->actnum = NULL;
}

1
mxgrdecl.h
View File

@@ -36,7 +36,6 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
#define OPENRS_MXGRDECL_HEADER
void mx_init_grdecl (struct grdecl *g, const mxArray *s);
void free_grdecl (struct grdecl *g);
#endif // OPENRS_MXGRDECL_HEADER

63
preprocess.c
View File

@@ -401,20 +401,57 @@ static void compute_node_coordinates(const struct grdecl *g,
/*-----------------------------------------------------------------
Public interface
*/
void process_grdecl(const struct grdecl *g,
void process_grdecl(const struct grdecl *in,
double tolerance,
struct processed_grid *out)
{
int i;
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;
g.dims[0] = nx;
g.dims[1] = ny;
g.dims[2] = nz;
int *actnum = malloc (nc * sizeof( *actnum));
double *zcorn = malloc (nc * 8 * sizeof (*zcorn));
/* Permute actnum */
int *iptr = actnum;
for (j=0; j<ny; ++j){
for (i=0; i<nx; ++i){
for (k=0; k<nz; ++k){
*iptr++ = in->actnum[i+nx*(j+ny*k)];
}
}
}
g.actnum = actnum;
/* Permute zcorn */
double *dptr = zcorn;
for (j=0; j<2*ny; ++j){
for (i=0; i<2*nx; ++i){
for (k=0; k<2*nz; ++k){
*dptr++ = in->zcorn[i+2*nx*(j+2*ny*k)];
}
}
}
g.zcorn = zcorn;
g.coord = in->coord;
/* Code below assumes k index runs fastests, ie. that dimensions of
table is permuted to (dims[2], dims[0], dims[1]) */
int nx = g->dims[0];
int ny = g->dims[1];
int nz = g->dims[2];
int *cell_index = calloc(nx*ny*nz,sizeof(*cell_index));
@@ -436,11 +473,10 @@ void process_grdecl(const struct grdecl *g,
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);
finduniquepoints(&g, plist, pillarz, tolerance);
npillarpoints = pillarz->ptr[npillars];
void *p = realloc_sparse_table (pillarz, npillars, npillarpoints, sizeof(double));
if (p) {
pillarz = p;
@@ -483,11 +519,11 @@ void process_grdecl(const struct grdecl *g,
faces->ptr[0] = 0;
/* faces with constant i */
process_vertical_faces(g->dims, 0, faces,
process_vertical_faces(g.dims, 0, faces,
&neighbors, &intersections, &npoints,
npillarpoints, plist, work);
/* faces with constant j */
process_vertical_faces(g->dims, 1, faces,
process_vertical_faces(g.dims, 1, faces,
&neighbors, &intersections, &npoints,
npillarpoints, plist, work);
@@ -495,7 +531,7 @@ void process_grdecl(const struct grdecl *g,
process_horizontal_faces(g, cell_index, &ncells, faces, &neighbors,
process_horizontal_faces(&g, cell_index, &ncells, faces, &neighbors,
&intersections, plist);
@@ -527,18 +563,19 @@ void process_grdecl(const struct grdecl *g,
}
}
/* 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,
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;

1
processgrid.c
View File

@@ -159,6 +159,5 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
/* Free whatever was allocated in initGrdecl. */
free_grdecl(&g);
free_processed_grid(&out);
}

51
uniquepoints.c
View File

@@ -50,16 +50,18 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
#define overlap(a1,a2,b1,b2) max(a1,b1) < min(a2,b2)
/* Compar function passed to qsort */
/*-------------------------------------------------------*/
/*-----------------------------------------------------------------
Compare function passed to qsort
*/
static int compare(const void *a, const void *b)
{
if (*(double*)a < *(double*) b) return -1;
else return 1;
}
/* Creat sorted list of z-values in zcorn with actnum==1 */
/*-------------------------------------------------------*/
/*-----------------------------------------------------------------
Creat sorted list of z-values in zcorn with actnum==1
*/
static int createSortedList(double *list, int n, int m,
const double *z[], const int *a[])
{
@@ -77,9 +79,10 @@ static int createSortedList(double *list, int n, int m,
}
/* Remove points that are closer than tolerance in list */
/* of increasing doubles */
/*-------------------------------------------------------*/
/*-----------------------------------------------------------------
Remove points that are closer than tolerance in list
of increasing doubles
*/
static int uniquify(int n, double *list, double tolerance)
{
if (n<1) return 0;
@@ -103,7 +106,9 @@ 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,
@@ -157,7 +162,13 @@ static int assignPointNumbers(int begin,
}
/*-------------------------------------------------------*/
/*-----------------------------------------------------------------
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 igetvectors(const int dims[3], int i, int j,
const int *field, const int *v[])
{
@@ -173,8 +184,15 @@ static void igetvectors(const int dims[3], int i, int j,
v[3] = field + dims[2]*(ip + dims[0]* jp);
}
/*-------------------------------------------------------*/
static void dgetvectors(const int dims[3], int i, int j, const double *field, const double *v[])
/*-----------------------------------------------------------------
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[])
{
int im = max(1, i ) - 1;
@@ -188,10 +206,11 @@ static void dgetvectors(const int dims[3], int i, int j, const double *field, co
v[3] = field + dims[2]*(ip + dims[0]* jp);
}
/* Assign point numbers p such that "zlist(p)==zcorn". */
/* Assume that coordinate number is arranged in a */
/* sequence such that the natural index is (k,i,j) */
/*-------------------------------------------------------*/
/*-----------------------------------------------------------------
Assign point numbers p such that "zlist(p)==zcorn".
Assume that coordinate number is arranged in a
sequence such that the natural index is (k,i,j)
*/
int finduniquepoints(const struct grdecl *g,
/* return values: */
int *plist, /* list of point numbers on each pillar*/
@@ -200,8 +219,6 @@ int finduniquepoints(const struct grdecl *g,
{
tolerance = tolerance < DBL_EPSILON ? DBL_EPSILON : tolerance;
double *zlist = ztab->data; /* casting void* to double* */
int *zptr = ztab->ptr;