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

Refactoring code, adding correct handling of faces on boundary.

Browse Source 
Packing unstructured data in csr-like strucuture.
This commit is contained in:
Jostein R. Natvig 2009-06-11 23:35:19 +00:00
parent 82ff742b37
commit f02fbe0b9d
7 changed files with 204 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
facetopology.c
View File

@ -3,6 +3,7 @@
#include <string.h>
#include <assert.h>
#include "sparsetable.h"
#include "facetopology.h"
/* No checking of input arguments in this code! */
@ -94,7 +95,7 @@ static int *computeFaceTopology(int *a1,
each pillar (but only separately), we can use only the point indices.
b) We assume no intersections occur on the first and last lines.
This is convenient in the identification of (unique) intersections.
This is convenient in the identification of (unique) intersections.
*/
@ -108,13 +109,12 @@ static int faceintersection(int *a1, int *a2, int *b1, int *b2)
}
/* work should be pointer to 2n ints initialised to zero . */
void findconnections(int n, int *pts[4],
int *ptnumber,
void findconnections(int n, int *pts[4],
int *ptnumber,
int *intersectionlist,
int *neighbors,
int *faces,
int *fptr,
int *fpos, int *work)
int *work,
sparse_table_t *ftab)
{
/* vectors of point numbers for faces a(b) on pillar 1(2) */
int *a1 = pts[0];
@ -125,8 +125,9 @@ 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 = faces + fptr[*fpos];
int *c = neighbors;
/* int *f = (int*)ftab->data + ftab->ptr[*fpos]; */
int *f = (int*)ftab->data + ftab->ptr[ftab->position];
int *c = neighbors + 2*ftab->position;
int k1 = 0;
int k2 = 0;
@ -174,7 +175,7 @@ void findconnections(int n, int *pts[4],
*f++ = a1[i+1];
*f++ = a2[i+1];
*f++ = a2[i];
fptr[++(*fpos)] = f-faces;
ftab->ptr[++ftab->position] = f - (int*)ftab->data;
}
}
@ -206,7 +207,7 @@ 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))){
f = computeFaceTopology(a1+i, a2+i, b1+j, b2+j, intersect, f);
fptr[++(*fpos)] = f-faces;
ftab->ptr[++ftab->position] = f - (int*)ftab->data;
}
}
}
@ -229,4 +230,3 @@ void findconnections(int n, int *pts[4],
/* Now, j = min(k1, k2) and itop is all zeros */
}
}

9
facetopology.h
View File

@ -2,12 +2,11 @@
#define FACETOPOLOGY_H
void findconnections(int n, int *pts[4],
int *ptnumber,
void findconnections(int n, int *pts[4],
int *ptnumber,
int *intersectionlist,
int *neighbors,
int *faces,
int *fptr,
int *fpos, int *work);
int *work,
sparse_table_t *ftab);
#endif

143
processgrid.c
View File

@ -4,7 +4,7 @@
#include <assert.h>
#include <mex.h>
#include "sparsetable.h"
#include "grdecl.h"
#include "uniquepoints.h"
#include "mxgrdecl.h"
@ -33,7 +33,7 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
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);
@ -41,6 +41,8 @@ static void igetvectors(int dims[3], int i, int j, int *field, int *v[])
}
/* Gateway routine for Matlab mex function. */
/*-------------------------------------------------------*/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
@ -59,32 +61,15 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
/* 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));
/* Unstructured storage of unique zcorn values for each pillar. */
/* ztab->data may need extra space temporarily due to simple boundary treatement */
int sz = 8*(g.dims[0]+1)*(g.dims[1]+1)*g.dims[2]; /* Big for convenient bc */
int numpillars = (g.dims[0]+1)*(g.dims[1]+1);
sparse_table_t *ztab = malloc_sparse_table(numpillars, sz, sizeof(double));
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
finduniquepoints(&g, plist, ztab);
/*======================================================*/
@ -96,69 +81,109 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int k;
const int BIGNUM = 100;
int *faces = malloc(BIGNUM*sz*sizeof(int));
int *fptr = malloc(BIGNUM*sz*sizeof(int));
/* Unstructured storage of face nodes */
sparse_table_t *ftab = malloc_sparse_table(BIGNUM*sz*sizeof(int),
BIGNUM*sz*sizeof(int),
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];
int startpts = numpts;
fptr[0]=0;
int startface;
int numpts = ztab->ptr[numpillars];
int startpts = numpts;
int isodd;
ftab->ptr[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);
startface = ftab->position;
findconnections(n, pts, &numpts, intersections, neighbors, work, ftab);
/* odd */
for(k=2*startface+1; k<2*ftab->position; k+=2){
if (neighbors[k] != -1){
neighbors[k] = g.dims[0]*(j + g.dims[1]*neighbors[k]);
}
}
/* even */
for(k=2*startface; k<2*ftab->position; k+=2){
neighbors[k] = -1;
}
/* -------------- */
/* internal faces */
/* -------------- */
for (i=1; i<g.dims[0]; ++i){
/* Vectors of point numbers */
igetvectors(d, 2*i, 2*j+1, plist, pts);
int startfaces = numfaces;
findconnections(n, pts, &numpts, intersections,
neighbors+2*numfaces, faces, fptr, &numfaces,
work);
startface = ftab->position;
findconnections(n, pts, &numpts, intersections, neighbors, work, ftab);
/* Compute cell numbers from cell k-index (stored in neighbors) */
for (k=startfaces; 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]);
/* For even k, neighbors[k] refers to stack (i-1,j),
for odd k, neighbors[k] refers to stack (i,j) of cells */
isodd = 1;
for(k=2*startface; k<2*ftab->position; ++k){
isodd = !isodd;
if (neighbors[k] != -1){
neighbors[k] = i-1+isodd + g.dims[0]*(j + g.dims[1]*neighbors[k]);
}
}
/* compute intersections */
}
/* ------------------ */
/* Add boundary face */
/* ------------------ */
startface = ftab->position;
igetvectors(d, 2*g.dims[0], 2*j+1, plist, pts);
findconnections(n, pts, &numpts, intersections,
neighbors+2*numfaces, faces, fptr, &numfaces,
work);
findconnections(n, pts, &numpts, intersections, neighbors, work, ftab);
/* even indices */
for(k=2*startface; k<2*ftab->position; k+=2){
if (neighbors[k] != -1){
neighbors[k] = g.dims[0]-1 + g.dims[0]*(j + g.dims[1]*neighbors[k]);
}
}
/* odd indices */
for(k=2*startface+1; k<2*ftab->position; k+=2){
neighbors[k] = -1;
}
}
/* */
/* */
/* */
/* D E B U G CODE */
/* */
/* */
/* */
#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, "\nfaces\nnumfaces %d\n", ftab->position);
for (i=0; i<ftab->position; ++i){
for (k=ftab->ptr[i]; k<ftab->ptr[i+1]; ++k){
fprintf(stderr, "%d ", ((int*)ftab->data)[k]);
}
fprintf(stderr, "\n");
}
@ -166,7 +191,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
fprintf(stderr, "\nneighbors\n");
int *iptr = neighbors;
for(k=0; k<numfaces; ++k){
for(k=0; k<ftab->position; ++k){
fprintf(stderr, " (%d %d)\n", iptr[0], iptr[1]);
++iptr;
++iptr;
@ -198,13 +223,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
}
#endif
free_sparse_table(ztab);
free_sparse_table(ftab);
free (intersections);
free (faces);
free (fptr);
free (neighbors);
free (zptr);
free (plist);
free (zlist);
/* Free whatever was allocated in initGrdecl. */
freeGrdecl(&g);

58
sparsetable.c Normal file
View File

@ -0,0 +1,58 @@
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <mex.h>
#include "sparsetable.h"
void free_sparse_table (sparse_table_t *tab)
{
if(tab->ptr) free(tab->ptr);
if(tab->data) free(tab->data);
free(tab);
}
sparse_table_t *malloc_sparse_table(int m, int n, int datasz)
{
sparse_table_t *tab = malloc(sizeof *tab);
tab->m = m;
tab->n = n;
tab->position = 0;
/* fprintf(stderr, "sizeof tab = %ld\n", sizeof (*tab->ptr)); */
if (!(tab->ptr = malloc((m+1) * sizeof (*tab->ptr)))){
fprintf(stderr, "Could not allocate space for sparse ptr\n");
free_sparse_table(tab);
return NULL;
}
if(!(tab->data = malloc(n * datasz))){
fprintf(stderr, "Could not allocate space for sparse data\n");
free_sparse_table(tab);
return NULL;
}
return tab;
}
sparse_table_t *realloc_sparse_table(sparse_table_t *tab, int m, int n, int datasz)
{
tab->m = m;
tab->n = n;
if (!(tab->ptr = realloc(tab->ptr, (m+1) * sizeof (*tab->ptr)))){
fprintf(stderr, "Could not reallocate space for sparse ptr\n");
free_sparse_table(tab);
return NULL;
}
if(!(tab->data = realloc(tab->data, n * datasz))){
fprintf(stderr, "Could not reallocate space for sparse data\n");
free_sparse_table(tab);
return NULL;
}
return tab;
}

22
sparsetable.h Normal file
View File

@ -0,0 +1,22 @@
#ifndef SPARSETABLE_H
#define SPARSETABLE_H
typedef struct{
int m; /* number of rows */
int *ptr; /* row pointer of size m+1 */
int position; /* position in ptr that is not filled. */
int n; /* size of data */
void *data; /* sparse table data */
} sparse_table_t;
void free_sparse_table (sparse_table_t *tab);
sparse_table_t *malloc_sparse_table (int m, int n, int datasz);
sparse_table_t *realloc_sparse_table (sparse_table_t *tab, int m, int n, int datasz);
#endif

44
uniquepoints.c
View File

@ -5,11 +5,11 @@
#include <stdio.h>
#include "sparsetable.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)
@ -69,13 +69,13 @@ static int uniquify(int n, double *list, double tolerance)
/* Along single pillar: */
static int* assignPointNumbers(int begin,
int end,
double *zlist,
int n,
double *zcorn,
int *actnum,
int *plist,
double tolerance)
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 */
@ -109,7 +109,6 @@ static int* assignPointNumbers(int begin,
if (k == end || z[i] - zlist[k] > tolerance){
fprintf(stderr, "What!?\n");
}
*p++ = k;
}
*p++ = INT_MAX;/* Padding to ease processing of faults */
@ -154,20 +153,24 @@ static void dgetvectors(int dims[3], int i, int j, double *field, double *v[])
/*-------------------------------------------------------*/
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 *plist, /* list of point numbers on each pillar*/
sparse_table_t *ztab)
{
double *zlist = ztab->data; /* casting void* to double* */
int *zptr = ztab->ptr;
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;
int pos = 0;
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){
@ -183,9 +186,8 @@ void finduniquepoints(struct Grdecl *g,
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;
zout = zout + len;
zptr[++pos] = zout - zlist;
}
}
@ -193,8 +195,6 @@ void finduniquepoints(struct Grdecl *g,
/* 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){
@ -210,8 +210,10 @@ void finduniquepoints(struct Grdecl *g,
int *a = g->actnum + cix;
double *z = g->zcorn + zix;
p = assignPointNumbers(zptr[pix], zptr[pix+1], zlist,
nz, z, a, p, 0.0);
assignPointNumbers(zptr[pix], zptr[pix+1], zlist,
2*g->dims[2], z, a, p, 0.0);
p += 2 + 2*g->dims[2];
}
}
}

5
uniquepoints.h
View File

@ -1,8 +1,7 @@
#ifndef UNIQUEPOINTS_H
#define UNIQUEPOINTS_H
void finduniquepoints(struct Grdecl *g,
double *zout,
int *zptr,
int *plist);
int *plist,
sparse_table_t *ztab);
#endif