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

Tag faces (possible values: {LEFT,BACK,TOP}) during processing.

Browse Source 
This is needed to support non-trivial boundaryId() return values.
This commit is contained in:
Bård Skaflestad committed 2009-07-03 18:18:59 +00:00
1 parent 8bae310891
commit 987d59c92c
3 files changed
+140 -114
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

preprocess.c
+91 -75
View File
@@ -20,7 +20,7 @@ This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
@@ -50,9 +50,9 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
/*-----------------------------------------------------------------
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
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(int dims[3], int i, int j, int *field, int *v[])
@@ -77,12 +77,12 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
/*-----------------------------------------------------------------
Special purpose
Convert from k-index to Cartesian index i+nx*(j+ny*k) for every other
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)
{
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]){
@@ -96,7 +96,7 @@ void compute_cell_index(const int dims[3], int i, int j, int *neighbors, int len
neighbors[k] = tmp;
}
}
}
}
}
@@ -110,31 +110,35 @@ int checkmemeory(int nz, struct processed_grid *out, int **intersections)
int r = (2*nz+2)*(2*nz+2);
int m = out->m;
int n = out->n;
if(out->number_of_faces + r> m){
if(out->number_of_faces + r > m){
m += max(m*0.5, 2*r);
}
if (out->face_ptr[out->number_of_faces] + 6*r > n){
n += max(n*0.5, 12*r);
}
if (m != out->m){
void *p1 = realloc(out->face_neighbors, 2*m * sizeof(*out->face_neighbors));
void *p1 = realloc(out->face_neighbors, 2*m * sizeof(*out->face_neighbors));
void *p2 = realloc(*intersections, 4*m * sizeof(**intersections));
if (p1 && p2){
out->face_neighbors = p1;
if (p1 && p2) {
out->face_neighbors = p1;
*intersections = p2;
}else{
} else {
return 0;
}
}
if (m != out->m || n != out->n){
void *p1 = realloc(out->face_nodes, n*sizeof(double));
void *p2 = realloc(out->face_ptr, m*sizeof(double));
if (p1 && p2){
if (m != out->m || n != out->n) {
void *p1 = realloc(out->face_nodes, n * sizeof *out->face_nodes);
void *p2 = realloc(out->face_ptr, (m+1) * sizeof *out->face_ptr);
void *p3 = realloc(out->face_tag, m * sizeof *out->face_tag);
if (p1 && p2 && p3) {
out->face_nodes = p1;
out->face_ptr = p2;
}else{
out->face_tag = p3;
} else {
return 0;
}
@@ -146,22 +150,24 @@ int checkmemeory(int nz, struct processed_grid *out, int **intersections)
}
/*-----------------------------------------------------------------
For each vertical face (i.e. i or j constant),
For each vertical face (i.e. i or j constant),
-find point numbers for the corners and
-cell neighbors.
-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(int direction,
int **intersections,
int *plist, int *work,
void process_vertical_faces(int direction,
int **intersections,
int *plist, int *work,
struct processed_grid *out)
{
int i,j;
int *cornerpts[4];
enum face_tag tag[] = { LEFT, BACK };
int nx = out->dimensions[0];
int ny = out->dimensions[1];
int nz = out->dimensions[2];
@@ -171,13 +177,13 @@ void process_vertical_faces(int direction,
if (!checkmemeory(nz, out, intersections)){
fprintf(stderr, "Could not allocat enough space in process_vertical_faces\n");
exit(1);
exit(1);
}
/* Vectors of point numbers */
int d[] = {2*nx, 2*ny, 2+2*nz};
igetvectors(d, 2*i+direction, 2*j+1-direction, plist, cornerpts);
if(direction==1){
/* 1 3 0 1 */
/* ---> */
@@ -195,17 +201,22 @@ void process_vertical_faces(int direction,
/* int num_intersections = *npoints - npillarpoints; */
int num_intersections = out->number_of_nodes -
out->number_of_nodes_on_pillars;
findconnections(2*nz+2, cornerpts,
*intersections+4*num_intersections,
findconnections(2*nz+2, cornerpts,
*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(out->dimensions, i-1+direction, j-direction, ptr, len);
compute_cell_index(out->dimensions, i, j, ptr+1, len);
/* Tag the new faces */
int f = startface;
for (; f < out->number_of_faces; ++f) {
out->face_tag[f] = tag[direction];
}
}
}
}
@@ -217,17 +228,17 @@ static int linearindex(const int dims[3], int i, int j, int k)
/*-----------------------------------------------------------------
For each horizontal face (i.e. k constant),
For each horizontal face (i.e. k constant),
-find point numbers for the corners and
-cell neighbors.
-cell neighbors.
Also define map from logically Cartesian
cell index to local cell index 0, ..., #<active cells>. Exclude
cells that are have collapsed coordinates. (This includes cells with
Also define map from logically Cartesian
cell index to local cell index 0, ..., #<active cells>. Exclude
cells that are have collapsed coordinates. (This includes cells with
ACTNUM==0)
*/
void process_horizontal_faces(int **intersections,
void process_horizontal_faces(int **intersections,
int *plist,
struct processed_grid *out)
{
@@ -239,7 +250,7 @@ void process_horizontal_faces(int **intersections,
int *cell = out->local_cell_index;
int cellno = 0;
/* dimensions of plist */
int d[] = {2*nx, 2*ny, 2+2*nz};
@@ -247,21 +258,21 @@ void process_horizontal_faces(int **intersections,
for(j=0; j<ny; ++j) {
for (i=0; i<nx; ++i) {
if (!checkmemeory(nz, out, intersections)){
fprintf(stderr, "Could not allocat enough space in process_horizontal_faces\n");
exit(1);
}
int *f = out->face_nodes + out->face_ptr[out->number_of_faces];
int *n = out->face_neighbors + 2*out->number_of_faces;
int *f = out->face_nodes + out->face_ptr[out->number_of_faces];
int *n = out->face_neighbors + 2*out->number_of_faces;
/* Vectors of point numbers */
int *c[4];
igetvectors(d, 2*i+1, 2*j+1, plist, c);
int prevcell = -1;
@@ -272,11 +283,11 @@ void process_horizontal_faces(int **intersections,
/* 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 the pinch is a cell: */
if (k%2){
int idx = linearindex(out->dimensions, i,j,(k-1)/2);
cell[idx] = -1;
cell[idx] = -1;
}
}
else{
@@ -288,14 +299,15 @@ void process_horizontal_faces(int **intersections,
*f++ = c[3][k];
*f++ = c[1][k];
out->face_tag[ out->number_of_faces] = TOP;
out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
int thiscell = linearindex(out->dimensions, i,j,(k-1)/2);
*n++ = prevcell;
*n++ = prevcell = thiscell;
cell[thiscell] = cellno++;
}
else{
if (prevcell != -1){
@@ -305,7 +317,9 @@ void process_horizontal_faces(int **intersections,
*f++ = c[3][k];
*f++ = c[1][k];
out->face_tag[ out->number_of_faces] = TOP;
out->face_ptr[++out->number_of_faces] = f - out->face_nodes;
*n++ = prevcell;
*n++ = prevcell = -1;
}
@@ -319,14 +333,14 @@ void process_horizontal_faces(int **intersections,
/*-----------------------------------------------------------------
On input,
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].
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)
{
@@ -352,8 +366,8 @@ 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_intersection_coordinates(int *intersections,
static void
compute_intersection_coordinates(int *intersections,
struct processed_grid *out)
{
int n = out->number_of_nodes;
@@ -369,7 +383,7 @@ compute_intersection_coordinates(int *intersections,
fprintf(stderr, "Could not allocate extra space for intersections\n");
}
/* Append intersections */
int k;
double *pt = out->node_coordinates + 3*np;
@@ -387,24 +401,24 @@ compute_intersection_coordinates(int *intersections,
/*-----------------------------------------------------------------
Public interface
*/
void process_grdecl(const struct grdecl *in,
double tolerance,
void process_grdecl(const struct grdecl *in,
double tolerance,
struct processed_grid *out)
{
int i,j,k;
int nx = in->dims[0];
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));
int *actnum = malloc (nc * sizeof *actnum);
double *zcorn = malloc (nc * 8 * sizeof *zcorn);
/* Permute actnum */
int *iptr = actnum;
@@ -438,22 +452,23 @@ void process_grdecl(const struct grdecl *in,
out->m = BIGNUM/3;
out->n = BIGNUM;
out->face_neighbors = malloc(BIGNUM* sizeof(*out->face_neighbors));
out->face_nodes = malloc(out->n * sizeof(int));
out->face_ptr = malloc(out->m * sizeof(int));
out->face_neighbors = malloc( BIGNUM * sizeof *out->face_neighbors);
out->face_nodes = malloc( out->n * sizeof *out->face_nodes);
out->face_ptr = malloc((out->m + 1) * sizeof *out->face_ptr);
out->face_tag = malloc( out->m * sizeof *out->face_tag);
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_nodes = 0;
out->number_of_cells = 0;
out->node_coordinates = NULL;
out->local_cell_index = malloc(nx*ny*nz *sizeof(int));
out->local_cell_index = malloc(nx*ny*nz * sizeof *out->local_cell_index);
/* internal */
int *intersections = malloc(BIGNUM* sizeof(*intersections));
@@ -462,24 +477,24 @@ void process_grdecl(const struct grdecl *in,
for(i=0; i<4*(nz+1); ++i) work[i] = -1;
/* Allocate space for cornerpoint numbers plus INT_MIN (INT_MAX) padding */
int *plist = malloc( 4*nx*ny*(2*nz+2) * sizeof(int));
int *plist = malloc( 4*nx*ny*(2*nz+2) * sizeof *plist);
/* Do actual work here:*/
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_vertical_faces (1, &intersections, plist, work, out);
process_horizontal_faces (&intersections, plist, out);
free (plist);
free (work);
compute_intersection_coordinates(intersections, out);
compute_intersection_coordinates(intersections, out);
free (intersections);
@@ -491,9 +506,9 @@ void process_grdecl(const struct grdecl *in,
*ptr = out->local_cell_index[*ptr];
}
}
/* Invert global-to-local map */
int *global_cell_index = malloc(out->number_of_cells *
int *global_cell_index = malloc(out->number_of_cells *
sizeof (*global_cell_index));
for (i=0; i<nx*ny*nz; ++i){
if(out->local_cell_index[i]!=-1){
@@ -508,8 +523,9 @@ void process_grdecl(const struct grdecl *in,
void free_processed_grid(struct processed_grid *g)
{
if( g ){
free ( g->face_nodes );
free ( g->face_nodes );
free ( g->face_ptr );
free ( g->face_tag );
free ( g->face_neighbors );
free ( g->node_coordinates );
free ( g->local_cell_index );
preprocess.h
+40 -37
View File
@@ -13,23 +13,23 @@
//===========================================================================
/*
Copyright 2009 SINTEF ICT, Applied Mathematics.
Copyright 2009 Statoil ASA.
Copyright 2009 SINTEF ICT, Applied Mathematics.
Copyright 2009 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_PREPROCESS_HEADER
@@ -41,37 +41,40 @@ extern "C" {
#endif
/* Input structure holding raw cornerpoint spec. */
struct grdecl{
int dims[3];
const double *coord;
const double *zcorn;
const int *actnum;
};
struct grdecl{
int dims[3];
const double *coord;
const double *zcorn;
const int *actnum;
};
/* Constant: I J K */
enum face_tag { LEFT, BACK, TOP };
/* 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 */
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 */
enum face_tag *face_tag;
int number_of_nodes;
int number_of_nodes_on_pillars;
double *node_coordinates; /* 3 doubles per node, sequentially */
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 */
int *local_cell_index; /* Global to local map */
};
int number_of_cells; /* number of active cells */
int *local_cell_index; /* Global to local map */
};
void process_grdecl (const struct grdecl *g,
double tol,
struct processed_grid *out);
void free_processed_grid(struct processed_grid *g);
void process_grdecl (const struct grdecl *g,
double tol,
struct processed_grid *out);
void free_processed_grid(struct processed_grid *g);
#ifdef __cplusplus
processgrid.c
+9 -2
View File
@@ -35,8 +35,8 @@ void fill_grid(mxArray **out, struct processed_grid *grid)
/* faces */
const char *n3[] = {"num", "neighbors", "numNodes"};
mxArray *faces = mxCreateStructMatrix(1,1,3,n3);
const char *n3[] = {"num", "neighbors", "numNodes", "tag"};
mxArray *faces = mxCreateStructMatrix(1,1,4,n3);
mxSetField(faces, 0, "num", mxCreateDoubleScalar(grid->number_of_faces));
@@ -63,6 +63,13 @@ void fill_grid(mxArray **out, struct processed_grid *grid)
}
mxSetField(faces, 0, "numNodes", numnodes);
mxArray *tags = mxCreateDoubleMatrix(grid->number_of_faces, 1, mxREAL);
ptr = mxGetPr(tags);
for (i = 0; i < grid->number_of_faces; ++i) {
ptr[i] = grid->face_tag[i] + 1;
}
mxSetField(faces, 0, "tag", tags);
mxSetField(G, 0, "faces", faces);