346 lines
9.2 KiB
C
346 lines
9.2 KiB
C
/*===========================================================================
|
|
//
|
|
// File: uniquepoints.c
|
|
//
|
|
// Created: Fri Jun 19 08:46:05 2009
|
|
//
|
|
// Author: Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
|
|
//
|
|
// $Date$
|
|
//
|
|
// $Revision$
|
|
//
|
|
//==========================================================================*/
|
|
|
|
/*
|
|
Copyright 2009 SINTEF ICT, Applied Mathematics.
|
|
Copyright 2009 Statoil ASA.
|
|
|
|
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 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/>.
|
|
*/
|
|
|
|
#include <assert.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include <string.h>
|
|
#include <limits.h>
|
|
#include <float.h>
|
|
#include <stdio.h>
|
|
|
|
|
|
#include "sparsetable.h"
|
|
#include "preprocess.h"
|
|
#include "uniquepoints.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)
|
|
|
|
|
|
/*-----------------------------------------------------------------
|
|
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
|
|
*/
|
|
static int createSortedList(double *list, int n, int m,
|
|
const double *z[], const int *a[])
|
|
{
|
|
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 less than <tolerance> apart in <list> of increasing
|
|
doubles.
|
|
*/
|
|
static int uniquify(int n, double *list, double tolerance)
|
|
{
|
|
assert (!(tolerance < 0.0));
|
|
|
|
if (n<1) return 0;
|
|
int i;
|
|
int pos = 0;
|
|
double val = list[pos++];/* Keep first value */
|
|
|
|
for (i=1; i<n; ++i){
|
|
if (val + tolerance < list [i]){
|
|
val = list[i];
|
|
list[pos++] = val;
|
|
}
|
|
}
|
|
|
|
/*
|
|
Preserve outer z-boundary.
|
|
|
|
This operation is a no-op in the case
|
|
|
|
list[pos-2] + tolerance < list[n-1].
|
|
|
|
If, however, the second to last point is less than <tolerance>
|
|
away from the last point (list[n-1]), we remove this
|
|
second-to-last point as it cannot be distinguished from "final"
|
|
point.
|
|
*/
|
|
list[pos-1] = list[n-1];
|
|
|
|
return pos;
|
|
}
|
|
|
|
|
|
/*-----------------------------------------------------------------
|
|
Along single pillar:
|
|
*/
|
|
static int assignPointNumbers(int begin,
|
|
int end,
|
|
const double *zlist,
|
|
int n,
|
|
const double *zcorn,
|
|
const 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. */
|
|
|
|
|
|
const double *z = zcorn;
|
|
const 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) || ( zlist[k] + tolerance < z[i])){
|
|
fprintf(stderr, "Cannot associate zcorn values with given list\n");
|
|
fprintf(stderr, "of z-coordinates to given tolerance\n");
|
|
return 0;
|
|
}
|
|
|
|
*p++ = k;
|
|
}
|
|
*p++ = INT_MAX;/* Padding to ease processing of faults */
|
|
|
|
|
|
return 1;
|
|
}
|
|
|
|
|
|
|
|
/*-----------------------------------------------------------------
|
|
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[])
|
|
{
|
|
|
|
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);
|
|
}
|
|
|
|
|
|
/*-----------------------------------------------------------------
|
|
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;
|
|
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);
|
|
}
|
|
|
|
|
|
/*-----------------------------------------------------------------
|
|
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].
|
|
*/
|
|
static void interpolate_pillar(const double *coord, double *pt)
|
|
{
|
|
double a = (pt[2]-coord[2])/(coord[5]-coord[2]);
|
|
if (isinf(a) || isnan(a)){
|
|
a = 0;
|
|
}
|
|
|
|
pt[0] = coord[0] + a*(coord[3]-coord[0]);
|
|
pt[1] = coord[1] + a*(coord[4]-coord[1]);
|
|
}
|
|
|
|
/*-----------------------------------------------------------------
|
|
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*/
|
|
double tolerance,
|
|
struct processed_grid *out)
|
|
|
|
{
|
|
int nx = out->dimensions[0];
|
|
int ny = out->dimensions[1];
|
|
int nz = out->dimensions[2];
|
|
|
|
|
|
/* ztab->data may need extra space temporarily due to simple boundary treatement */
|
|
int npillarpoints = 8*(nx+1)*(ny+1)*nz;
|
|
int npillars = (nx+1)*(ny+1);
|
|
sparse_table_t *ztab = malloc_sparse_table(npillars,
|
|
npillarpoints,
|
|
sizeof(double));
|
|
|
|
|
|
|
|
int nc = g->dims[0]*g->dims[1]*g->dims[2];
|
|
out->node_coordinates = malloc (3*8*nc*sizeof(*out->node_coordinates));
|
|
|
|
double *zlist = ztab->data; /* casting void* to double* */
|
|
int *zptr = ztab->ptr;
|
|
|
|
int i,j,k;
|
|
|
|
int d1[3] = {2*g->dims[0], 2*g->dims[1], 2*g->dims[2]};
|
|
int len = 0;
|
|
double *zout = zlist;
|
|
int pos = 0;
|
|
|
|
zptr[pos++] = zout - zlist;
|
|
|
|
double *coord = (double*)g->coord;
|
|
double *pt = out->node_coordinates;
|
|
|
|
/* 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){
|
|
|
|
const int *a[4];
|
|
const 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, tolerance);
|
|
|
|
/* Assign unique points */
|
|
for (k=0; k<len; ++k){
|
|
pt[2] = zout[k];
|
|
interpolate_pillar(coord, pt);
|
|
pt += 3;
|
|
}
|
|
|
|
/* Increment pointer to sparse table of unique zcorn values */
|
|
zout = zout + len;
|
|
zptr[pos++] = zout - zlist;
|
|
|
|
coord += 6;
|
|
}
|
|
}
|
|
out->number_of_nodes_on_pillars = zptr[pos-1];
|
|
out->number_of_nodes = zptr[pos-1];
|
|
|
|
/* Loop over all vertical sets of zcorn values, assign point numbers */
|
|
int *p = plist;
|
|
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);
|
|
|
|
const int *a = g->actnum + cix;
|
|
const double *z = g->zcorn + zix;
|
|
|
|
if (!assignPointNumbers(zptr[pix], zptr[pix+1], zlist,
|
|
2*g->dims[2], z, a, p, tolerance)){
|
|
fprintf(stderr, "Something went wrong in assignPointNumbers");
|
|
return 0;
|
|
}
|
|
|
|
p += 2 + 2*g->dims[2];
|
|
}
|
|
}
|
|
|
|
|
|
free_sparse_table(ztab);
|
|
|
|
|
|
return 1;
|
|
}
|
|
|
|
|