Initial checkin of cornerpoint processing code.
Current status -------------- Given vectors ZCORN, COORD and ACTNUM as well as the Cartesian dimensions these vectors implicitly refer to, the code is currently capable of * Identify unique points along each pillar * Assign point numbers for each point specified in ZCORN * Compute face topology, i.e., the corners that define the geometry of the faces as well as the cells that are connected through the face. * Identify and compute intesections that occur in the processing of face topology. What remains is * Handle the face geometry of boundary faces. (simple) * Compute point coordinates of the final point list. * Put all pieces together in a tidy manner.
This commit is contained in:
Jostein R. Natvig
commit
e04eb2f728
13
grdecl.h
Normal file
13
grdecl.h
Normal file
|
|
@ -0,0 +1,13 @@
|
|||
#ifndef GRDECL_H
|
||||
#define GRDECL_H
|
||||
|
||||
struct Grdecl{
|
||||
int n;
|
||||
int dims[3];
|
||||
double *coord;
|
||||
double *zcorn;
|
||||
int *actnum;
|
||||
};
|
||||
|
||||
|
||||
#endif
|
||||
24
matalloc.c
Normal file
24
matalloc.c
Normal file
|
|
@ -0,0 +1,24 @@
|
|||
#include <stdlib.h>
|
||||
#include "matalloc.h"
|
||||
|
||||
double **dmatalloc(int m, int n)
|
||||
{
|
||||
double **ret = malloc(m*sizeof *ret);
|
||||
ret[0] = malloc(m*n*sizeof *ret[0]);
|
||||
int i;
|
||||
for (i=1; i<m; ++i) {
|
||||
ret[i] = ret[i-1]+n*sizeof *ret[0];
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
|
||||
int **imatalloc(int m, int n)
|
||||
{
|
||||
int **ret = malloc(m*sizeof *ret);
|
||||
ret[0] = malloc(m*n*sizeof *ret[0]);
|
||||
int i;
|
||||
for (i=1; i<m; ++i) {
|
||||
ret[i] = ret[i-1]+n*sizeof *ret[0];
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
7
matalloc.h
Normal file
7
matalloc.h
Normal file
|
|
@ -0,0 +1,7 @@
|
|||
#ifndef MATALLOC_H_
|
||||
#define MATALLOC_H
|
||||
|
||||
double **dmatalloc(int m, int n);
|
||||
int **imatalloc(int m, int n);
|
||||
|
||||
#endif
|
||||
68
mxgrdecl.c
Normal file
68
mxgrdecl.c
Normal file
|
|
@ -0,0 +1,68 @@
|
|||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
#include <assert.h>
|
||||
#include <mex.h>
|
||||
|
||||
|
||||
#include "grdecl.h"
|
||||
|
||||
|
||||
|
||||
/* Get COORD, ZCORN, ACTNUM and DIMS from mxArray. */
|
||||
/*-------------------------------------------------------*/
|
||||
void mxInitGrdecl(struct Grdecl *g, const mxArray *prhs[])
|
||||
{
|
||||
int i,j,k;
|
||||
|
||||
|
||||
g->coord = mxGetPr(mxGetField(prhs[0], 0, "COORD"));
|
||||
double *tmp = mxGetPr(mxGetField(prhs[0], 0, "cartDims"));
|
||||
g->n = 1;
|
||||
for (i=0; i<3; ++i){
|
||||
g->dims[i] = tmp[i];
|
||||
g->n *= tmp[i];
|
||||
}
|
||||
mexPrintf("dimensions: %d %d %d\n",
|
||||
g->dims[0],
|
||||
g->dims[1],
|
||||
g->dims[2]);
|
||||
|
||||
|
||||
|
||||
/* grdecl.actnum = permute(actnum, [3,1,2]); */
|
||||
int *actnum = mxGetData(mxGetField(prhs[0], 0, "ACTNUM"));
|
||||
g->actnum = malloc(g->n* sizeof(*g->actnum));
|
||||
int *iptr = g->actnum;
|
||||
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)];
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* grdecl.zcorn = permute(zcorn, [3,1,2]); */
|
||||
double *zcorn = mxGetPr(mxGetField(prhs[0], 0, "ZCORN"));
|
||||
g->zcorn = malloc(g->n*8*sizeof(*g->zcorn));
|
||||
double *dptr = g->zcorn;
|
||||
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)];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* Free stuff that was allocated in initgrdecl. */
|
||||
/*-------------------------------------------------------*/
|
||||
void freeGrdecl(struct Grdecl *g)
|
||||
{
|
||||
free(g->zcorn); g->zcorn = NULL;
|
||||
free(g->actnum); g->actnum = NULL;
|
||||
}
|
||||
|
||||
7
mxgrdecl.h
Normal file
7
mxgrdecl.h
Normal file
|
|
@ -0,0 +1,7 @@
|
|||
#ifndef MXGRDECL_H
|
||||
#define MXGRDECL_H
|
||||
void mxInitGrdecl (struct Grdecl *g, const mxArray *prhs[]);
|
||||
void freeGrdecl (struct Grdecl *g);
|
||||
|
||||
#endif
|
||||
|
||||
437
processgrid.c
Normal file
437
processgrid.c
Normal file
|
|
@ -0,0 +1,437 @@
|
|||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
#include <assert.h>
|
||||
#include <mex.h>
|
||||
|
||||
|
||||
#include "grdecl.h"
|
||||
#include "uniquepoints.h"
|
||||
#include "mxgrdecl.h"
|
||||
#include "matalloc.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))
|
||||
|
||||
|
||||
|
||||
/*------------------------------------------------------*/
|
||||
/* */
|
||||
/* */
|
||||
/* Find connections for each pair of pillars */
|
||||
/* */
|
||||
/* */
|
||||
/* */
|
||||
/*------------------------------------------------------*/
|
||||
|
||||
|
||||
/* Determine face geometry first, then compute intersections. */
|
||||
/* All intersections that occur are present in the final face geometry.*/
|
||||
int *computeFaceTopology(int *a1,
|
||||
int *a2,
|
||||
int *b1,
|
||||
int *b2,
|
||||
int intersect[4],
|
||||
int *faces)
|
||||
{
|
||||
int mask[8];
|
||||
|
||||
/* Which pillar points should we use? */
|
||||
if (a1[1] > b1[1]){ mask[0] = b1[1]; } else { mask[0] = a1[1]; }
|
||||
if (a2[1] > b2[1]){ mask[2] = b2[1]; } else { mask[2] = a2[1]; }
|
||||
if (a2[0] > b2[0]){ mask[4] = a2[0]; } else { mask[4] = b2[0]; }
|
||||
if (a1[0] > b1[0]){ mask[6] = a1[0]; } else { mask[6] = b1[0]; }
|
||||
|
||||
/* Get shape of face: */
|
||||
/* each new intersection will be part of the new face, */
|
||||
/* 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;
|
||||
|
||||
/* bottom-top */
|
||||
if (intersect[1]){
|
||||
if(a1[0] > b1[1]){ /* intersection[1] left of (any) intersection[0] */
|
||||
mask[0] = 0;
|
||||
mask[6] = 0;
|
||||
mask[7] = intersect[1];
|
||||
}
|
||||
else{
|
||||
mask[2] = 0;
|
||||
mask[4] = 0;
|
||||
mask[3] = intersect[1];
|
||||
}
|
||||
}
|
||||
|
||||
/* top-bottom */
|
||||
if (intersect[2]){
|
||||
if(a1[1] < b1[0]){ /* intersection[2] left of (any) intersection[3] */
|
||||
mask[0] = 0;
|
||||
mask[6] = 0;
|
||||
mask[7] = intersect[2];
|
||||
}
|
||||
else{
|
||||
mask[2] = 0;
|
||||
mask[4] = 0;
|
||||
mask[3] = intersect[2];
|
||||
}
|
||||
}
|
||||
|
||||
int k;
|
||||
int *f = faces;
|
||||
for (k=0; k<8; ++k){
|
||||
if(mask[k]){
|
||||
*f++ = mask[k];
|
||||
}
|
||||
}
|
||||
|
||||
return f;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/* values in za and zb must be increasing. */
|
||||
|
||||
|
||||
/* a) If we assume that the index increase when z increase for
|
||||
each pillar (but only separately), we can skip the z
|
||||
|
||||
b) We assume no intersections occur on the first and last lines.
|
||||
|
||||
|
||||
|
||||
*/
|
||||
/* #define overlap(a1,a2,b1,b2) max(a1,b1) < min(a2,b2) */
|
||||
#define intersection(a1,a2,b1,b2)(((a1>b1)&&(a2<b2))||((a1<b1)&&(a2>b2)))
|
||||
|
||||
int faceintersection(int *a1, int *a2, int *b1, int *b2)
|
||||
{
|
||||
/* overlap(a1[0], a1[1], b1[0], b1[1]) || */
|
||||
/* overlap(a2[0], a2[1], b2[0], b2[1]) || */
|
||||
|
||||
return
|
||||
max(a1[0],b1[0]) < min(a1[1], b1[1]) ||
|
||||
max(a2[0],b2[0]) < min(a2[1],b2[1]) ||
|
||||
intersection(a1[0], a2[0], b1[0], b2[0]);
|
||||
}
|
||||
|
||||
/* work should be pointer to 2n ints initialised to zero . */
|
||||
void findconnections(int n, int *pts[4], int *ptnumber, int *intersections,
|
||||
int *neighbors,
|
||||
int *faces, int *fptr, int *fpos, int *work)
|
||||
{
|
||||
/* vectors of point numbers for faces a(b) on pillar 1(2) */
|
||||
int *a1 = pts[0];
|
||||
int *a2 = pts[1];
|
||||
int *b1 = pts[2];
|
||||
int *b2 = pts[3];
|
||||
|
||||
/* Intersection record for top line and bottomline of a */
|
||||
int *itop = work; /* calloc(n, sizeof(*itop)); */
|
||||
int *ibottom = work + n; /* calloc(n, sizeof(*ibottom)); */
|
||||
int *f = faces + fptr[*fpos];
|
||||
int *c = neighbors;
|
||||
|
||||
int k1 = 0;
|
||||
int k2 = 0;
|
||||
|
||||
int i,j=0;
|
||||
int intersect[4];
|
||||
|
||||
for (i = 0; i<n-1; ++i){
|
||||
if (a1[i] == a1[i+1] && a2[i] == a2[i+1]) continue;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
while(j<n-1 && (b1[j] < a1[i+1] || b2[j] < a2[i+1])){
|
||||
|
||||
if (b1[j] == b1[j+1] && b2[j] == b2[j+1]){
|
||||
itop[j+1] = itop[j];
|
||||
++j;
|
||||
continue;
|
||||
}
|
||||
|
||||
|
||||
/* --------------------------------------------------------- */
|
||||
/* face a(i,i+1) and face b(j,j+1) have nonzero intersection */
|
||||
/* --------------------------------------------------------- */
|
||||
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))){
|
||||
fprintf(stderr, "here: %d %d\n", i%2 ? (i-1)/2 : -1,
|
||||
j%2 ? (j-1)/2 : -1);
|
||||
*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))){
|
||||
*f++ = a1[i];
|
||||
*f++ = a1[i+1];
|
||||
*f++ = a2[i+1];
|
||||
*f++ = a2[i];
|
||||
fptr[++(*fpos)] = f-faces;
|
||||
}
|
||||
}
|
||||
|
||||
/* Non-matching faces */
|
||||
else{
|
||||
|
||||
/* Find new intersection */
|
||||
if (intersection(a1[i+1],a2[i+1],b1[j+1],b2[j+1])) {
|
||||
itop[j+1] = (*ptnumber)++;
|
||||
|
||||
/* store point numbers of intersecting lines */
|
||||
*intersections++ = a1[i+1];
|
||||
*intersections++ = a2[i+1];
|
||||
*intersections++ = b1[j+1];
|
||||
*intersections++ = b2[j+1];
|
||||
|
||||
|
||||
}else{
|
||||
itop[j+1] = 0;
|
||||
}
|
||||
|
||||
/* Update intersection record */
|
||||
intersect[0] = ibottom[j ]; /* i x j */
|
||||
intersect[1] = ibottom[j+1]; /* i x j+1 */
|
||||
intersect[2] = itop[j ]; /* i+1 x j */
|
||||
intersect[3] = itop[j+1]; /* i+1 x 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))){
|
||||
f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f);
|
||||
fptr[++(*fpos)] = f-faces;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Update candidates for restart of j for in next i-iteration */
|
||||
if (b1[j] < a1[i+1]) k1 = j;
|
||||
if (b2[j] < a2[i+1]) k2 = j;
|
||||
|
||||
j = j+1;
|
||||
}
|
||||
|
||||
/* Swap intersection records: top line of a is next bottom line of a */
|
||||
int *tmp;
|
||||
tmp = itop; itop = ibottom; ibottom = tmp;
|
||||
|
||||
for(;j>min(k1,k2);--j) itop[j-1]=0;
|
||||
/* Now, j = min(k1, k2) */
|
||||
|
||||
}
|
||||
|
||||
/* free(itop); */
|
||||
/* free(ibottom); */
|
||||
/* return (cells-begin)/2; */
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
void interpolate_pillar(double *coord, double *x, double *y, double *z)
|
||||
{
|
||||
double a = (*z-coord[2])/(coord[5]-coord[2]);
|
||||
*x = coord[0] + a*(coord[3]-coord[0]);
|
||||
*y = coord[1] + a*(coord[4]-coord[1]);
|
||||
}
|
||||
|
||||
|
||||
/*-------------------------------------------------------*/
|
||||
static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
|
||||
{
|
||||
|
||||
int im = max(1, i ) - 1;
|
||||
int ip = min(dims[0], i+1) - 1;
|
||||
int jm = max(1, j ) - 1;
|
||||
int jp = min(dims[1], j+1) - 1;
|
||||
fprintf(stderr, "%d %d %d %d\n", im,ip, jm,jp);
|
||||
v[0] = field + dims[2]*(im + dims[0]* jm);
|
||||
v[1] = field + dims[2]*(im + dims[0]* jp);
|
||||
v[2] = field + dims[2]*(ip + dims[0]* jm);
|
||||
v[3] = field + dims[2]*(ip + dims[0]* jp);
|
||||
}
|
||||
|
||||
|
||||
/* Gateway routine for Matlab mex function. */
|
||||
/*-------------------------------------------------------*/
|
||||
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
||||
{
|
||||
int i,j;
|
||||
|
||||
/* Set up data passed from Matlab */
|
||||
struct Grdecl g;
|
||||
mxInitGrdecl(&g, prhs);
|
||||
|
||||
|
||||
/* Code below assumes k index runs fastests, ie. that dimensions of
|
||||
table is permuted to (dims[2], dims[0], dims[1]) */
|
||||
/* ---------------------------------------------------------------*/
|
||||
|
||||
|
||||
/* Set up space for return values */
|
||||
|
||||
/* zlist may need extra space temporarily due to simple boundary treatement */
|
||||
|
||||
int sz = 8*(g.dims[0]+1)*(g.dims[1]+1)*g.dims[2];
|
||||
int numpillars = (g.dims[0]+1)*(g.dims[1]+1);
|
||||
double *zlist = malloc(sz*sizeof(*zlist));
|
||||
int *zptr = malloc((numpillars+1)*sizeof(*zptr));
|
||||
int *plist = malloc( 2*g.dims[0]*2*g.dims[1]*(2*g.dims[2]+2)*sizeof(int));
|
||||
|
||||
|
||||
fprintf(stderr, "Allocate %d ints for plist\n", 2*g.dims[0]*2*g.dims[1]*(2*g.dims[2]+2));
|
||||
|
||||
finduniquepoints(&g, zlist, zptr, plist);
|
||||
|
||||
|
||||
#if 0
|
||||
for (i=0; i<numpillars ; ++i){
|
||||
fprintf(stderr, "pillar %d\n", i);
|
||||
for (k=zptr[i]; k<zptr[i+1]; ++k){
|
||||
fprintf(stderr, "%f\n", zlist[k]);
|
||||
}
|
||||
|
||||
mexPrintf("\n");
|
||||
}
|
||||
|
||||
for (i=0; i<8*g.n; ++i) mexPrintf("%d\n", plist[i]);
|
||||
#endif
|
||||
|
||||
|
||||
/*======================================================*/
|
||||
|
||||
fprintf(stderr, "process face geomtery\n");
|
||||
/* Process face geometry and cell-face topology */
|
||||
/* internal constant i pillar pairs */
|
||||
int n = 2 + 2*g.dims[2];
|
||||
int k;
|
||||
|
||||
const int BIGNUM = 100;
|
||||
int *faces = malloc(BIGNUM*sz*sizeof(int));
|
||||
int *fptr = malloc(BIGNUM*sz*sizeof(int));
|
||||
|
||||
|
||||
|
||||
int *pts[4];
|
||||
int *neighbors = malloc(2* n*n* sizeof(*neighbors));
|
||||
int *intersections = malloc(4* n*n* sizeof(*intersections));
|
||||
int *work = calloc(2* n, sizeof(*work));
|
||||
int d[3] = {2*g.dims[0], 2*g.dims[1], 2+2*g.dims[2]};
|
||||
|
||||
int numfaces = 0;
|
||||
int numpts = zptr[numpillars];
|
||||
|
||||
fptr[0]=0;
|
||||
|
||||
for (j=0; j<g.dims[1]; ++j) {
|
||||
/* Add boundary faces */
|
||||
igetvectors(d, 0, 2*j+1, plist, pts);
|
||||
findconnections(n, pts, &numpts, intersections,
|
||||
neighbors+2*numfaces, faces, fptr, &numfaces,
|
||||
work);
|
||||
|
||||
/* internal faces */
|
||||
for (i=1; i<g.dims[0]; ++i){
|
||||
|
||||
/* Vectors of point numbers */
|
||||
igetvectors(d, 2*i, 2*j+1, plist, pts);
|
||||
|
||||
int start = numpts;
|
||||
findconnections(n, pts, &numpts, intersections,
|
||||
neighbors+2*numfaces, faces, fptr, &numfaces,
|
||||
work);
|
||||
|
||||
|
||||
|
||||
/* Compute cell numbers from cell k-index (stored in neighbors) */
|
||||
for (k=0; k<numfaces; ++k){
|
||||
if (neighbors[2*k] != -1){
|
||||
neighbors[2*k] = i-1 + g.dims[0]*(j + g.dims[1]*neighbors[2*k]);
|
||||
}
|
||||
if (neighbors[2*k+1] != -1){
|
||||
|
||||
neighbors[2*k+1] = i + g.dims[0]*(j + g.dims[1]*neighbors[2*k+1]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
#if 1
|
||||
fprintf(stderr, "\nfaces\nnumfaces %d\n", numfaces);
|
||||
for (i=0; i<numfaces; ++i){
|
||||
for (k=fptr[i]; k<fptr[i+1]; ++k){
|
||||
fprintf(stderr, "%d ", faces[k]);
|
||||
}
|
||||
fprintf(stderr, "\n");
|
||||
}
|
||||
|
||||
|
||||
fprintf(stderr, "\nneighbors\n");
|
||||
int *iptr = neighbors;
|
||||
for(k=0; k<numfaces; ++k){
|
||||
fprintf(stderr, " (%d %d)\n", iptr[0], iptr[1]);
|
||||
iptr+=2;
|
||||
}
|
||||
|
||||
fprintf(stderr, "\nline intersections:\n");
|
||||
iptr = intersections;
|
||||
int numintersections = numpts - start;
|
||||
for(k=0; k<numintersections; ++k){
|
||||
fprintf(stderr, " (%d %d %d %d)\n", iptr[0], iptr[1], iptr[2], iptr[3]);
|
||||
iptr+=4;
|
||||
}
|
||||
#endif
|
||||
|
||||
/* compute intersections */
|
||||
|
||||
}
|
||||
|
||||
/* Add boundary face */
|
||||
|
||||
|
||||
}
|
||||
#if 0
|
||||
int xstride = 2*g.dims[2];
|
||||
int ystride = 4*g.dims[2]*g.dims[0];
|
||||
int zstride = 1;
|
||||
|
||||
/* internal constant j pillar pairs */
|
||||
for (i=0; i<g.dims[0]; ++i){
|
||||
for (j=0; j<g.dims[1]-1; ++j){
|
||||
int p1 = 2*i + ystride* (2*j+1);
|
||||
int p2 = p1 + xstride;
|
||||
process_pair(p1, p2, g.dims, g.zcorn, ystride, zstride);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
free (intersections);
|
||||
free (faces);
|
||||
free (fptr);
|
||||
free (neighbors);
|
||||
free (zptr);
|
||||
free (plist);
|
||||
free (zlist);
|
||||
/* Free whatever was allocated in initGrdecl. */
|
||||
freeGrdecl(&g);
|
||||
|
||||
}
|
||||
|
||||
6
test.m
Normal file
6
test.m
Normal file
|
|
@ -0,0 +1,6 @@
|
|||
|
||||
g=simpleGrdecl([2, 1, 4], @(x) -0.05+0.1*x+0.01 );
|
||||
G=processGRDECL(g);
|
||||
%clf,plotGrid(G);view(3);
|
||||
|
||||
g.ACTNUM=int32(g.ACTNUM);
|
||||
219
uniquepoints.c
Normal file
219
uniquepoints.c
Normal file
|
|
@ -0,0 +1,219 @@
|
|||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
#include <limits.h>
|
||||
#include <stdio.h>
|
||||
|
||||
|
||||
#include "grdecl.h"
|
||||
#include "uniquepoints.h"
|
||||
#include "matalloc.h"
|
||||
|
||||
|
||||
#define min(i,j) ((i)<(j) ? (i) : (j))
|
||||
#define max(i,j) ((i)>(j) ? (i) : (j))
|
||||
#define overlap(a1,a2,b1,b2) max(a1,b1) < min(a2,b2)
|
||||
|
||||
|
||||
/* Compar 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 */
|
||||
/*-------------------------------------------------------*/
|
||||
static int createSortedList(double *list, int n, int m,
|
||||
double *z[], int *a[])
|
||||
{
|
||||
fprintf(stderr, "\n");
|
||||
int i,j;
|
||||
double *ptr = list;
|
||||
for (i=0; i<n; ++i){
|
||||
for (j=0; j<m; ++j){
|
||||
if (a[j][i/2]) *ptr++ = z[j][i];
|
||||
/* else fprintf(stderr, "skipping point in inactive cell\n"); */
|
||||
}
|
||||
}
|
||||
|
||||
qsort(list, ptr-list, sizeof(double), compare);
|
||||
return ptr-list;
|
||||
}
|
||||
|
||||
|
||||
/* Remove points that are closer than tolerance in list */
|
||||
/* of increasing doubles */
|
||||
/*-------------------------------------------------------*/
|
||||
static int uniquify(int n, double *list, double tolerance)
|
||||
{
|
||||
int i;
|
||||
int pos = 1; /* Keep first value */
|
||||
double val = list[pos];
|
||||
for (i=1; i<n; ++i){
|
||||
if (list[i] - val > tolerance){
|
||||
val = list[i];
|
||||
list[pos++] = val;
|
||||
}
|
||||
}
|
||||
|
||||
/* Keep last value (one way or the other...) */
|
||||
if (list[n-1] - val > tolerance){
|
||||
list[pos-1] = list[n-1];
|
||||
}
|
||||
|
||||
return pos;
|
||||
}
|
||||
|
||||
|
||||
/* Along single pillar: */
|
||||
static int* assignPointNumbers(int begin,
|
||||
int end,
|
||||
double *zlist,
|
||||
int n,
|
||||
double *zcorn,
|
||||
int *actnum,
|
||||
int *plist,
|
||||
double tolerance)
|
||||
{
|
||||
/* n - number of cells */
|
||||
/* zlist - list of len unique z-values */
|
||||
/* start - number of unique z-values processed before. */
|
||||
|
||||
int i, k;
|
||||
/* All points should now be within tolerance of a listed point. */
|
||||
|
||||
|
||||
double *z = zcorn;
|
||||
int *a = actnum;
|
||||
int *p = plist;
|
||||
|
||||
k = begin;
|
||||
*p++ = INT_MIN; /* Padding to ease processing of faults */
|
||||
for (i=0; i<n; ++i){
|
||||
|
||||
/* Skip inactive cells */
|
||||
if (!a[i/2]) {
|
||||
p[0] = p[-1]; /* Inactive cells are collapsed leaving void space.*/
|
||||
++p;
|
||||
continue;
|
||||
}
|
||||
|
||||
/* Find next k such that zlist[k] < z[i] < zlist[k+1] */
|
||||
while (k < end && zlist[k] + tolerance < z[i]){
|
||||
k++;
|
||||
}
|
||||
|
||||
/* assert (k < len && z[i] - zlist[k] <= tolerance) */
|
||||
if (k == end || z[i] - zlist[k] > tolerance){
|
||||
fprintf(stderr, "What!?\n");
|
||||
}
|
||||
|
||||
*p++ = k;
|
||||
}
|
||||
*p++ = INT_MAX;/* Padding to ease processing of faults */
|
||||
|
||||
return p;
|
||||
}
|
||||
|
||||
|
||||
/*-------------------------------------------------------*/
|
||||
static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
|
||||
{
|
||||
|
||||
int im = max(1, i ) - 1;
|
||||
int ip = min(dims[0], i+1) - 1;
|
||||
int jm = max(1, j ) - 1;
|
||||
int jp = min(dims[1], j+1) - 1;
|
||||
|
||||
v[0] = field + dims[2]*(im + dims[0]* jm);
|
||||
v[1] = field + dims[2]*(im + dims[0]* jp);
|
||||
v[2] = field + dims[2]*(ip + dims[0]* jm);
|
||||
v[3] = field + dims[2]*(ip + dims[0]* jp);
|
||||
}
|
||||
|
||||
/*-------------------------------------------------------*/
|
||||
static void dgetvectors(int dims[3], int i, int j, double *field, double *v[])
|
||||
{
|
||||
|
||||
int im = max(1, i ) - 1;
|
||||
int ip = min(dims[0], i+1) - 1;
|
||||
int jm = max(1, j ) - 1;
|
||||
int jp = min(dims[1], j+1) - 1;
|
||||
|
||||
v[0] = field + dims[2]*(im + dims[0]* jm);
|
||||
v[1] = field + dims[2]*(im + dims[0]* jp);
|
||||
v[2] = field + dims[2]*(ip + dims[0]* jm);
|
||||
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) */
|
||||
/*-------------------------------------------------------*/
|
||||
void finduniquepoints(struct Grdecl *g,
|
||||
/* return values: */
|
||||
double *zlist, /* list of z-values on each pillar*/
|
||||
int *zptr, /* pointer to start of each pillar*/
|
||||
int *plist) /* list of point numbers on each pillar*/
|
||||
|
||||
{
|
||||
int i,j;
|
||||
|
||||
int d1[3] = {2*g->dims[0], 2*g->dims[1], 2*g->dims[2]};
|
||||
int len = 0;
|
||||
double *zout = zlist;
|
||||
int pos = 1;
|
||||
|
||||
zptr[0] = 0;
|
||||
|
||||
/* 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){
|
||||
|
||||
int *a[4];
|
||||
double *z[4];
|
||||
|
||||
/* Get positioned pointers for actnum and zcorn data */
|
||||
igetvectors(g->dims, i, j, g->actnum, a);
|
||||
dgetvectors(d1, 2*i, 2*j, g->zcorn, z);
|
||||
|
||||
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] = zptr[pos-1] + len;
|
||||
++pos;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* Loop over all vertical sets of zcorn values, assign point numbers */
|
||||
int *p = plist;
|
||||
int nz = 2*g->dims[2];
|
||||
|
||||
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);
|
||||
|
||||
/* cell column position */
|
||||
int cix = g->dims[2]*((i/2) + (j/2)*g->dims[0]);
|
||||
|
||||
/* zcorn column position */
|
||||
int zix = 2*g->dims[2]*(i+2*g->dims[0]*j);
|
||||
|
||||
int *a = g->actnum + cix;
|
||||
double *z = g->zcorn + zix;
|
||||
|
||||
p = assignPointNumbers(zptr[pix], zptr[pix+1], zlist,
|
||||
nz, z, a, p, 0.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
8
uniquepoints.h
Normal file
8
uniquepoints.h
Normal file
|
|
@ -0,0 +1,8 @@
|
|||
#ifndef UNIQUEPOINTS_H
|
||||
#define UNIQUEPOINTS_H
|
||||
void finduniquepoints(struct Grdecl *g,
|
||||
double *zout,
|
||||
int *zptr,
|
||||
int *plist);
|
||||
|
||||
#endif
|
||||
Loading…
Reference in New Issue
Block a user