Import geometry calculation from the MATLAB Reservoir Simulation Toolbox.
Author: Jostein R. Natvig Hook up to build.
This commit is contained in:
Bård Skaflestad
parent
e80ecf1077
commit
dac3142d16
10
Makefile.am
10
Makefile.am
@ -15,12 +15,14 @@ noinst_LTLIBRARIES = libcpgpreprocess_noinst.la
|
||||
libcpgpreprocess_la_SOURCES = \
|
||||
cgridinterface.c \
|
||||
facetopology.c \
|
||||
uniquepoints.c \
|
||||
facetopology.h \
|
||||
geometry.c \
|
||||
geometry.h \
|
||||
preprocess.c \
|
||||
sparsetable.c \
|
||||
facetopology.h \
|
||||
uniquepoints.h \
|
||||
sparsetable.h
|
||||
sparsetable.h \
|
||||
uniquepoints.c \
|
||||
uniquepoints.h
|
||||
|
||||
libcpgpreprocess_noinst_la_SOURCES = \
|
||||
$(libcpgpreprocess_la_SOURCES)
|
||||
|
||||
@ -1,28 +1,30 @@
|
||||
#include <assert.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
#include <grid.h>
|
||||
|
||||
#include "geometry.h"
|
||||
#include "cgridinterface.h"
|
||||
|
||||
|
||||
|
||||
|
||||
static int *compute_cell_facepos(grid_t *g)
|
||||
{
|
||||
int i,j,k;
|
||||
int *facepos = malloc((g->number_of_cells + 1) * sizeof *facepos);
|
||||
int *fcells = g->face_cells;
|
||||
|
||||
|
||||
for (i=0; i<g->number_of_cells; ++i) {
|
||||
facepos [i] = 0;
|
||||
}
|
||||
|
||||
|
||||
for (i=0; i<2*g->number_of_faces; ++i) {
|
||||
if (*fcells != -1) {
|
||||
(facepos[*fcells])++;
|
||||
}
|
||||
fcells++;
|
||||
}
|
||||
|
||||
|
||||
/* cumsum */
|
||||
j=0;
|
||||
for (i=0; i<g->number_of_cells; ++i) {
|
||||
@ -45,10 +47,10 @@ static int *compute_cell_faces(grid_t *g)
|
||||
for(i=0; i<g->number_of_cells; ++i) {
|
||||
work[i] = 0;
|
||||
}
|
||||
|
||||
|
||||
for (i=0; i<g->number_of_faces; ++i) {
|
||||
for (k=0;k<2; ++k) {
|
||||
|
||||
|
||||
if (*fcells != -1) {
|
||||
cell = *fcells;
|
||||
cfaces[g->cell_facepos[cell] + work[cell]] = i;
|
||||
@ -57,13 +59,13 @@ static int *compute_cell_faces(grid_t *g)
|
||||
fcells++;
|
||||
}
|
||||
}
|
||||
free(work);
|
||||
free(work);
|
||||
|
||||
return cfaces;
|
||||
}
|
||||
|
||||
void preprocess (const struct grdecl *in,
|
||||
double tol,
|
||||
void preprocess (const struct grdecl *in,
|
||||
double tol,
|
||||
struct CornerpointGrid *G)
|
||||
{
|
||||
struct processed_grid pg;
|
||||
@ -73,8 +75,8 @@ void preprocess (const struct grdecl *in,
|
||||
|
||||
process_grdecl(in, tol, &pg);
|
||||
|
||||
/*
|
||||
* General grid interface
|
||||
/*
|
||||
* General grid interface
|
||||
*/
|
||||
base->dimensions = 3;
|
||||
|
||||
@ -87,11 +89,11 @@ void preprocess (const struct grdecl *in,
|
||||
base->face_nodes = pg.face_nodes;
|
||||
base->face_nodepos = pg.face_ptr;
|
||||
base->face_cells = pg.face_neighbors;
|
||||
|
||||
|
||||
base->face_centroids = NULL;
|
||||
base->face_normals = NULL;
|
||||
base->face_areas = NULL;
|
||||
|
||||
|
||||
/* NB: compute_cell_facepos must be called before compute_cell_faces */
|
||||
base->cell_facepos = compute_cell_facepos(base);
|
||||
base->cell_faces = compute_cell_faces (base);
|
||||
@ -99,8 +101,8 @@ void preprocess (const struct grdecl *in,
|
||||
base->cell_volumes = NULL;
|
||||
|
||||
|
||||
/*
|
||||
* Cornerpoint grid interface
|
||||
/*
|
||||
* Cornerpoint grid interface
|
||||
*/
|
||||
G->cartdims[0] = pg.dimensions[0];
|
||||
G->cartdims[1] = pg.dimensions[1];
|
||||
@ -114,7 +116,7 @@ void preprocess (const struct grdecl *in,
|
||||
|
||||
G->index_map = pg.local_cell_index;
|
||||
}
|
||||
|
||||
|
||||
void free_cornerpoint_grid(struct CornerpointGrid *G)
|
||||
{
|
||||
free(G->grid.face_nodes);
|
||||
@ -132,3 +134,53 @@ void free_cornerpoint_grid(struct CornerpointGrid *G)
|
||||
|
||||
free(G->index_map);
|
||||
}
|
||||
|
||||
static int
|
||||
allocate_geometry(struct CornerpointGrid *g)
|
||||
{
|
||||
int ok;
|
||||
size_t nc, nf, nd;
|
||||
|
||||
assert (g->grid.dimensions == 3);
|
||||
|
||||
nc = g->grid.number_of_cells;
|
||||
nf = g->grid.number_of_faces;
|
||||
nd = 3;
|
||||
|
||||
g->grid.face_areas = malloc(nf * 1 * sizeof *g->grid.face_areas);
|
||||
g->grid.face_centroids = malloc(nf * nd * sizeof *g->grid.face_centroids);
|
||||
g->grid.face_normals = malloc(nf * nd * sizeof *g->grid.face_normals);
|
||||
|
||||
g->grid.cell_volumes = malloc(nc * 1 * sizeof *g->grid.cell_volumes);
|
||||
g->grid.cell_centroids = malloc(nc * nd * sizeof *g->grid.cell_centroids);
|
||||
|
||||
ok = g->grid.face_areas != NULL;
|
||||
ok += g->grid.face_centroids != NULL;
|
||||
ok += g->grid.face_normals != NULL;
|
||||
|
||||
ok += g->grid.cell_volumes != NULL;
|
||||
ok += g->grid.cell_centroids != NULL;
|
||||
|
||||
return ok == 5;
|
||||
}
|
||||
|
||||
void compute_geometry(struct CornerpointGrid *g)
|
||||
{
|
||||
int ok;
|
||||
|
||||
ok = allocate_geometry(g);
|
||||
|
||||
if (ok) {
|
||||
compute_face_geometry(nd, g->grid.node_coordinates, nf,
|
||||
g->grid.face_nodepos, g->grid.face_nodes,
|
||||
g->grid.face_normals, g->grid.face_centroids,
|
||||
g->grid.face_areas);
|
||||
|
||||
compute_cell_geometry(nd, g->grid.node_coordinates,
|
||||
g->grid.face_nodepos, g->grid.face_nodes,
|
||||
g->grid.face_cells, g->grid.face_normals,
|
||||
g->grid.face_centroids, nc,
|
||||
g->grid.cell_facepos, g->grid.cell_faces,
|
||||
g->grid.cell_centroids, g->grid.cell_volumes);
|
||||
}
|
||||
}
|
||||
|
||||
@ -42,6 +42,8 @@ extern "C" {
|
||||
double tol,
|
||||
struct CornerpointGrid *out);
|
||||
|
||||
void compute_geometry (struct CornerpointGrid *g);
|
||||
|
||||
void free_cornerpoint_grid(struct CornerpointGrid *g);
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
|
||||
235
geometry.c
Normal file
235
geometry.c
Normal file
@ -0,0 +1,235 @@
|
||||
/*
|
||||
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
|
||||
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
|
||||
*/
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include "geometry.h"
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static void
|
||||
cross(const double u[3], const double v[3], double w[3])
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
w[0] = u[1]*v[2]-u[2]*v[1];
|
||||
w[1] = u[2]*v[0]-u[0]*v[2];
|
||||
w[2] = u[0]*v[1]-u[1]*v[0];
|
||||
}
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static double
|
||||
norm(const double w[3])
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
return sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
void
|
||||
compute_face_geometry(int ndims, double *coords, int nfaces,
|
||||
int *nodepos, int *facenodes, double *fnormals,
|
||||
double *fcentroids, double *fareas)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
|
||||
/* Assume 3D for now */
|
||||
int f;
|
||||
double x[3];
|
||||
double u[3];
|
||||
double v[3];
|
||||
double w[3];
|
||||
|
||||
int i,k;
|
||||
int node;
|
||||
|
||||
double cface[3] = {0};
|
||||
double n[3] = {0};
|
||||
const double twothirds = 0.666666666666666666666666666667;
|
||||
for (f=0; f<nfaces; ++f)
|
||||
{
|
||||
int num_face_nodes;
|
||||
double area = 0.0;
|
||||
for(i=0; i<ndims; ++i) x[i] = 0.0;
|
||||
for(i=0; i<ndims; ++i) n[i] = 0.0;
|
||||
for(i=0; i<ndims; ++i) cface[i] = 0.0;
|
||||
|
||||
/* average node */
|
||||
for(k=nodepos[f]; k<nodepos[f+1]; ++k)
|
||||
{
|
||||
node = facenodes[k];
|
||||
for (i=0; i<ndims; ++i) x[i] += coords[3*node+i];
|
||||
}
|
||||
num_face_nodes = nodepos[f+1] - nodepos[f];
|
||||
for(i=0; i<ndims; ++i) x[i] /= num_face_nodes;
|
||||
|
||||
|
||||
|
||||
/* compute first vector u (to the last node in the face) */
|
||||
node = facenodes[nodepos[f+1]-1];
|
||||
for(i=0; i<ndims; ++i) u[i] = coords[3*node+i] - x[i];
|
||||
|
||||
|
||||
/* Compute triangular contrib. to face normal and face centroid*/
|
||||
for(k=nodepos[f]; k<nodepos[f+1]; ++k)
|
||||
{
|
||||
double a;
|
||||
|
||||
node = facenodes[k];
|
||||
for (i=0; i<ndims; ++i) v[i] = coords[3*node+i] - x[i];
|
||||
|
||||
cross(u,v,w);
|
||||
a = 0.5*norm(w);
|
||||
area += a;
|
||||
if(!(a>0))
|
||||
{
|
||||
fprintf(stderr, "Internal error in compute_face_geometry.");
|
||||
}
|
||||
|
||||
/* face normal */
|
||||
for (i=0; i<ndims; ++i) n[i] += w[i];
|
||||
|
||||
/* face centroid */
|
||||
for (i=0; i<ndims; ++i)
|
||||
cface[i] += a*(x[i]+twothirds*0.5*(u[i]+v[i]));
|
||||
|
||||
/* Store v in u for next iteration */
|
||||
for (i=0; i<ndims; ++i) u[i] = v[i];
|
||||
}
|
||||
|
||||
/* Store face normal and face centroid */
|
||||
for (i=0; i<ndims; ++i)
|
||||
{
|
||||
/* normal is scaled with face area */
|
||||
fnormals [3*f+i] = 0.5*n[i];
|
||||
fcentroids[3*f+i] = cface[i]/area;
|
||||
}
|
||||
fareas[f] = area;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
void
|
||||
compute_cell_geometry(int ndims, double *coords,
|
||||
int *nodepos, int *facenodes, int *neighbors,
|
||||
double *fnormals,
|
||||
double *fcentroids,
|
||||
int ncells, int *facepos, int *cellfaces,
|
||||
double *ccentroids, double *cvolumes)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
|
||||
int i,k, f,c;
|
||||
int face,node;
|
||||
double x[3];
|
||||
double u[3];
|
||||
double v[3];
|
||||
double w[3];
|
||||
double xcell[3];
|
||||
double ccell[3];
|
||||
double cface[3] = {0};
|
||||
const double twothirds = 0.666666666666666666666666666667;
|
||||
for (c=0; c<ncells; ++c)
|
||||
{
|
||||
int num_faces;
|
||||
double volume = 0.0;
|
||||
|
||||
for(i=0; i<ndims; ++i) xcell[i] = 0.0;
|
||||
for(i=0; i<ndims; ++i) ccell[i] = 0.0;
|
||||
|
||||
|
||||
/*
|
||||
* Approximate cell center as average of face centroids
|
||||
*/
|
||||
for(f=facepos[c]; f<facepos[c+1]; ++f)
|
||||
{
|
||||
face = cellfaces[f];
|
||||
for (i=0; i<ndims; ++i) xcell[i] += fcentroids[3*face+i];
|
||||
}
|
||||
num_faces = facepos[c+1] - facepos[c];
|
||||
|
||||
for(i=0; i<ndims; ++i) xcell[i] /= num_faces;
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* For all faces, add tetrahedron's volume and centroid to
|
||||
* 'cvolume' and 'ccentroid'.
|
||||
*/
|
||||
for(f=facepos[c]; f<facepos[c+1]; ++f)
|
||||
{
|
||||
int num_face_nodes;
|
||||
|
||||
for(i=0; i<ndims; ++i) x[i] = 0.0;
|
||||
for(i=0; i<ndims; ++i) cface[i] = 0.0;
|
||||
|
||||
face = cellfaces[f];
|
||||
|
||||
/* average face node x */
|
||||
for(k=nodepos[face]; k<nodepos[face+1]; ++k)
|
||||
{
|
||||
node = facenodes[k];
|
||||
for (i=0; i<ndims; ++i) x[i] += coords[3*node+i];
|
||||
}
|
||||
num_face_nodes = nodepos[face+1] - nodepos[face];
|
||||
for(i=0; i<ndims; ++i) x[i] /= num_face_nodes;
|
||||
|
||||
|
||||
|
||||
/* compute first vector u (to the last node in the face) */
|
||||
node = facenodes[nodepos[face+1]-1];
|
||||
for(i=0; i<ndims; ++i) u[i] = coords[3*node+i] - x[i];
|
||||
|
||||
|
||||
/* Compute triangular contributions to face normal and face centroid */
|
||||
for(k=nodepos[face]; k<nodepos[face+1]; ++k)
|
||||
{
|
||||
double tet_volume, subnormal_sign;
|
||||
node = facenodes[k];
|
||||
for (i=0; i<ndims; ++i) v[i] = coords[3*node+i] - x[i];
|
||||
|
||||
cross(u,v,w);
|
||||
|
||||
|
||||
|
||||
tet_volume = 0;
|
||||
for(i=0; i<ndims; ++i){
|
||||
tet_volume += 0.5/3 * w[i]*(x[i]-xcell[i]);
|
||||
/*tet_volume += 0.5/3 * w[i]*(x[i]-xcell[i]);*/
|
||||
}
|
||||
/*tet_volume = fabs(tet_volume);*/
|
||||
|
||||
subnormal_sign=0.0;
|
||||
for(i=0; i<ndims; ++i){
|
||||
subnormal_sign += w[i]*fnormals[3*face+i];
|
||||
}
|
||||
|
||||
if(subnormal_sign<0){
|
||||
tet_volume*=-1.0;
|
||||
}
|
||||
if(!(neighbors[2*face+0]==c)){
|
||||
tet_volume*=-1.0;
|
||||
}
|
||||
|
||||
volume += tet_volume;
|
||||
if(!(volume>0)){
|
||||
fprintf(stderr, "Internal error in mex_compute_geometry: negative volume\n");
|
||||
}
|
||||
|
||||
/* face centroid of triangle */
|
||||
for (i=0; i<ndims; ++i) cface[i] = (x[i]+twothirds*0.5*(u[i]+v[i]));
|
||||
|
||||
/* Cell centroid */
|
||||
for (i=0; i<ndims; ++i) ccell[i] += tet_volume * 3/4.0*(cface[i] - xcell[i]);
|
||||
|
||||
|
||||
/* Store v in u for next iteration */
|
||||
for (i=0; i<ndims; ++i) u[i] = v[i];
|
||||
}
|
||||
}
|
||||
for (i=0; i<ndims; ++i) ccentroids[3*c+i] = xcell[i] + ccell[i]/volume;
|
||||
cvolumes[c] = volume;
|
||||
}
|
||||
}
|
||||
19
geometry.h
Normal file
19
geometry.h
Normal file
@ -0,0 +1,19 @@
|
||||
/*
|
||||
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
|
||||
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
|
||||
*/
|
||||
#ifndef MRST_GEOMETRY_H_INCLUDED
|
||||
#define MRST_GEOMETRY_H_INCLUDED
|
||||
|
||||
void compute_face_geometry(int ndims, double *coords, int nfaces,
|
||||
int *nodepos, int *facenodes,
|
||||
double *fnormals, double *fcentroids,
|
||||
double *fareas);
|
||||
void compute_cell_geometry(int ndims, double *coords,
|
||||
int *nodepos, int *facenodes, int *neighbours,
|
||||
double *fnormals,
|
||||
double *fcentroids, int ncells,
|
||||
int *facepos, int *cellfaces,
|
||||
double *ccentroids, double *cvolumes);
|
||||
|
||||
#endif /* MRST_GEOMETRY_H_INCLUDED */
|
||||
Loading…
Reference in New Issue
Block a user