Suggest an alternative implementation of the 'processgrid' function. This, disabled, implementation follows a two-stage allocate+fill algorithm that is slightly easier to control in memory-tight environments. We still have to maintain two copies of the grid, though.
Signed-off-by: Bård Skaflestad <Bard.Skaflestad@sintef.no>
This commit is contained in:
Bård Skaflestad
Bård Skaflestad
parent
fbe8157332
commit
b23ad048a2
408
processgrid.c
408
processgrid.c
@ -46,6 +46,412 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
|
||||
#include "mxgrdecl.h"
|
||||
|
||||
|
||||
#if 0
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static mxArray *
|
||||
allocate_nodes(size_t nnodes)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t nflds;
|
||||
const char *fields[] = { "num", "coords" };
|
||||
|
||||
mxArray *nodes, *num, *coords;
|
||||
|
||||
nflds = sizeof(fields) / sizeof(fields[0]);
|
||||
|
||||
nodes = mxCreateStructMatrix(1, 1, nflds, fields);
|
||||
|
||||
num = mxCreateDoubleScalar(nnodes);
|
||||
coords = mxCreateDoubleMatrix(nnodes, 3, mxREAL);
|
||||
|
||||
if ((nodes != NULL) && (num != NULL) && (coords != NULL)) {
|
||||
mxSetField(nodes, 0, "num" , num);
|
||||
mxSetField(nodes, 0, "coords", coords);
|
||||
} else {
|
||||
if (coords != NULL) { mxDestroyArray(coords); }
|
||||
if (num != NULL) { mxDestroyArray(num); }
|
||||
if (nodes != NULL) { mxDestroyArray(nodes); }
|
||||
|
||||
nodes = NULL;
|
||||
}
|
||||
|
||||
return nodes;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
fill_nodes(mxArray *nodes, struct processed_grid *grid)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t i, j, nnodes;
|
||||
double *coords;
|
||||
|
||||
nnodes = grid->number_of_nodes;
|
||||
coords = mxGetPr(mxGetField(nodes, 0, "coords"));
|
||||
|
||||
for (i = 0; i < nnodes; ++i) {
|
||||
for (j = 0; j < 3; ++j) {
|
||||
coords[i + j*nnodes] = grid->node_coordinates[3*i + j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static mxArray *
|
||||
allocate_faces(size_t nf, size_t nfacenodes)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t nflds;
|
||||
const char *fields[] = { "num", "neighbors", "nodePos", "nodes", "tag" };
|
||||
|
||||
mxArray *faces, *num, *neighbors, *nodePos, *nodes, *tag;
|
||||
|
||||
nflds = sizeof(fields) / sizeof(fields[0]);
|
||||
|
||||
faces = mxCreateStructMatrix(1, 1, nflds, fields);
|
||||
|
||||
num = mxCreateDoubleScalar (nf);
|
||||
neighbors = mxCreateNumericMatrix(nf , 2, mxINT32_CLASS, mxREAL);
|
||||
nodePos = mxCreateNumericMatrix(nf + 1 , 1, mxINT32_CLASS, mxREAL);
|
||||
nodes = mxCreateNumericMatrix(nfacenodes, 1, mxINT32_CLASS, mxREAL);
|
||||
tag = mxCreateNumericMatrix(nf , 1, mxINT32_CLASS, mxREAL);
|
||||
|
||||
if ((faces != NULL) && (num != NULL) && (neighbors != NULL) &&
|
||||
(nodePos != NULL) && (nodes != NULL) && (tag != NULL)) {
|
||||
mxSetField(faces, 0, "num" , num);
|
||||
mxSetField(faces, 0, "neighbors", neighbors);
|
||||
mxSetField(faces, 0, "nodePos" , nodePos);
|
||||
mxSetField(faces, 0, "nodes" , nodes);
|
||||
mxSetField(faces, 0, "tag" , tag);
|
||||
} else {
|
||||
if (tag != NULL) { mxDestroyArray(tag); }
|
||||
if (nodes != NULL) { mxDestroyArray(nodes); }
|
||||
if (nodePos != NULL) { mxDestroyArray(nodePos); }
|
||||
if (neighbors != NULL) { mxDestroyArray(neighbors); }
|
||||
if (num != NULL) { mxDestroyArray(num); }
|
||||
if (faces != NULL) { mxDestroyArray(faces); }
|
||||
|
||||
faces = NULL;
|
||||
}
|
||||
|
||||
return faces;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
fill_faces(mxArray *faces, struct processed_grid *grid)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t i, f, nf, nfn;
|
||||
|
||||
int *pi;
|
||||
|
||||
nf = grid->number_of_faces;
|
||||
|
||||
/* Fill faces.neighbors */
|
||||
pi = mxGetData(mxGetField(faces, 0, "neighbors"));
|
||||
for (f = 0; f < nf; f++) {
|
||||
/* Add one for one-based indexing in M */
|
||||
pi[f + 0*nf] = grid->face_neighbors[2*f + 0] + 1;
|
||||
pi[f + 1*nf] = grid->face_neighbors[2*f + 1] + 1;
|
||||
}
|
||||
|
||||
/* Fill faces.nodePos */
|
||||
pi = mxGetData(mxGetField(faces, 0, "nodePos"));
|
||||
for (i = 0; i <= nf; i++) { pi[i] = grid->face_ptr[i] + 1; }
|
||||
|
||||
/* Fill faces.nodes */
|
||||
pi = mxGetData(mxGetField(faces, 0, "nodes"));
|
||||
nfn = grid->face_ptr[nf]; /* Total number of face nodes */
|
||||
for (i = 0; i < nfn; i++) { pi[i] = grid->face_nodes[i] + 1; }
|
||||
|
||||
/* Fill faces.tag */
|
||||
pi = mxGetData(mxGetField(faces, 0, "tag"));
|
||||
for (f = 0; f < nf; f++) { pi[f] = grid->face_tag[i] + 1; }
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static size_t
|
||||
count_halffaces(size_t nf, const int *neighbors)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c1, c2;
|
||||
size_t nhf, f;
|
||||
|
||||
for (f = nhf = 0; f < nf; f++) {
|
||||
c1 = neighbors[2*f + 0];
|
||||
c2 = neighbors[2*f + 1];
|
||||
|
||||
nhf += c1 >= 0;
|
||||
nhf += c2 >= 0;
|
||||
}
|
||||
|
||||
return nhf;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static mxArray *
|
||||
allocate_cells(size_t nc, size_t ncf)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t nflds;
|
||||
const char *fields[] = { "num", "facePos", "faces", "indexMap" };
|
||||
|
||||
mxArray *cells, *num, *facePos, *faces, *indexMap;
|
||||
|
||||
nflds = sizeof(fields) / sizeof(fields[0]);
|
||||
|
||||
cells = mxCreateStructMatrix(1, 1, nflds, fields);
|
||||
|
||||
num = mxCreateDoubleScalar (nc);
|
||||
facePos = mxCreateNumericMatrix(nc + 1, 1, mxINT32_CLASS, mxREAL);
|
||||
faces = mxCreateNumericMatrix(ncf , 2, mxINT32_CLASS, mxREAL);
|
||||
indexMap = mxCreateNumericMatrix(nc , 1, mxINT32_CLASS, mxREAL);
|
||||
|
||||
if ((cells != NULL) && (num != NULL) && (facePos != NULL) &&
|
||||
(faces != NULL) && (indexMap != NULL)) {
|
||||
mxSetField(cells, 0, "num" , num );
|
||||
mxSetField(cells, 0, "facePos" , facePos );
|
||||
mxSetField(cells, 0, "faces" , faces );
|
||||
mxSetField(cells, 0, "indexMap", indexMap);
|
||||
} else {
|
||||
if (indexMap != NULL) { mxDestroyArray(indexMap); }
|
||||
if (faces != NULL) { mxDestroyArray(faces); }
|
||||
if (facePos != NULL) { mxDestroyArray(facePos); }
|
||||
if (num != NULL) { mxDestroyArray(num); }
|
||||
if (cells != NULL) { mxDestroyArray(cells); }
|
||||
|
||||
cells = NULL;
|
||||
}
|
||||
|
||||
return cells;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
fill_cells(mxArray *cells, struct processed_grid *grid)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t c, nc, f, nf, i;
|
||||
|
||||
int c1, c2, cf_tag, nhf;
|
||||
int *pi1, *pi2;
|
||||
|
||||
nc = grid->number_of_cells;
|
||||
nf = grid->number_of_faces;
|
||||
|
||||
/* Simultaneously fill cells.facePos and cells.faces by transposing the
|
||||
* neighbours mapping. */
|
||||
pi1 = mxGetData(mxGetField(cells, 0, "facePos"));
|
||||
pi2 = mxGetData(mxGetField(cells, 0, "faces" ));
|
||||
for (i = 0; i < nc + 1; i++) { pi1[i] = 0; }
|
||||
|
||||
/* 1) Count connections (i.e., faces per cell). */
|
||||
for (f = 0; f < nf; f++) {
|
||||
c1 = grid->face_neighbors[2*f + 0];
|
||||
c2 = grid->face_neighbors[2*f + 1];
|
||||
|
||||
if (c1 >= 0) { pi1[c1 + 1] += 1; }
|
||||
if (c2 >= 0) { pi1[c2 + 1] += 1; }
|
||||
}
|
||||
|
||||
/* 2) Define start pointers (really, position *end* pointers at start). */
|
||||
for (c = 1; c <= nc; c++) {
|
||||
pi1[0] += pi1[c];
|
||||
pi1[c] = pi1[0] - pi1[c];
|
||||
}
|
||||
|
||||
/* 3) Fill connection structure whilst advancing end pointers. */
|
||||
nhf = pi1[0];
|
||||
pi1[0] = 0;
|
||||
|
||||
mxAssert (((size_t) nhf) == mxGetM(mxGetField(cells, 0, "faces")),
|
||||
"Number of half faces (SIZE(cells.faces,1)) incorrectly "
|
||||
"determined earlier.");
|
||||
|
||||
for (f = 0; f < nf; f++) {
|
||||
cf_tag = 2*grid->face_tag[f] + 1; /* [1, 3, 5] */
|
||||
c1 = grid->face_neighbors[2*f + 0];
|
||||
c2 = grid->face_neighbors[2*f + 1];
|
||||
|
||||
if (c1 >= 0) {
|
||||
pi2[ pi1[ c1 + 1 ] + 0*nhf ] = f + 1;
|
||||
pi2[ pi1[ c1 + 1 ] + 1*nhf ] = cf_tag + 1; /* out */
|
||||
|
||||
pi1[ c1 + 1 ] += 1;
|
||||
}
|
||||
if (c2 >= 0) {
|
||||
pi2[ pi1[ c2 + 1 ] + 0*nhf ] = f + 1;
|
||||
pi2[ pi1[ c2 + 1 ] + 1*nhf ] = cf_tag + 0; /* in */
|
||||
|
||||
pi1[ c2 + 1 ] += 1;
|
||||
}
|
||||
}
|
||||
|
||||
/* Finally, adjust pointer array for one-based indexing in M. */
|
||||
for (i = 0; i < nc + 1; i++) { pi1[i] += 1; }
|
||||
|
||||
/* Fill cells.indexMap. Note that despite the name, 'local_cell_index'
|
||||
* really *is* the (zero-based) indexMap of the 'processed_grid'. */
|
||||
pi1 = mxGetData(mxGetField(cells, 0, "indexMap"));
|
||||
for (c = 0; c < nc; c++) { pi1[c] = grid->local_cell_index[c] + 1; }
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static mxArray *
|
||||
allocate_grid(struct processed_grid *grid)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t nflds, nhf;
|
||||
const char *fields[] = { "nodes", "faces", "cells", "type", "cartDims" };
|
||||
|
||||
mxArray *G, *nodes, *faces, *cells, *type, *cartDims;
|
||||
|
||||
nflds = sizeof(fields) / sizeof(fields[0]);
|
||||
nhf = count_halffaces(grid->number_of_faces, grid->face_neighbors);
|
||||
|
||||
G = mxCreateStructMatrix(1, 1, nflds, fields);
|
||||
|
||||
nodes = allocate_nodes(grid->number_of_nodes);
|
||||
faces = allocate_faces(grid->number_of_faces,
|
||||
grid->face_ptr[ grid->number_of_faces ]);
|
||||
cells = allocate_cells(grid->number_of_cells, nhf);
|
||||
type = mxCreateString("processgrid");
|
||||
cartDims = mxCreateDoubleMatrix(1, 3, mxREAL);
|
||||
|
||||
if ((G != NULL) && (nodes != NULL) && (faces != NULL) &&
|
||||
(cells != NULL) && (type != NULL) && (cartDims != NULL)) {
|
||||
mxSetField(G, 0, "nodes" , nodes );
|
||||
mxSetField(G, 0, "faces" , faces );
|
||||
mxSetField(G, 0, "cells" , cells );
|
||||
mxSetField(G, 0, "type" , type );
|
||||
mxSetField(G, 0, "cartDims", cartDims);
|
||||
} else {
|
||||
if (cartDims != NULL) { mxDestroyArray(cartDims); }
|
||||
if (type != NULL) { mxDestroyArray(type); }
|
||||
if (cells != NULL) { mxDestroyArray(cells); }
|
||||
if (faces != NULL) { mxDestroyArray(faces); }
|
||||
if (nodes != NULL) { mxDestroyArray(nodes); }
|
||||
if (G != NULL) { mxDestroyArray(G); }
|
||||
|
||||
G = NULL;
|
||||
}
|
||||
|
||||
return G;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
fill_grid(mxArray *G, struct processed_grid *grid)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
double *pr;
|
||||
|
||||
pr = mxGetPr(mxGetField(G, 0, "cartDims"));
|
||||
pr[0] = grid->dimensions[0];
|
||||
pr[1] = grid->dimensions[1];
|
||||
pr[2] = grid->dimensions[2];
|
||||
|
||||
fill_nodes(mxGetField(G, 0, "nodes"), grid);
|
||||
fill_faces(mxGetField(G, 0, "faces"), grid);
|
||||
fill_cells(mxGetField(G, 0, "cells"), grid);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
args_ok(int nlhs, int nrhs, const mxArray *prhs[])
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ok;
|
||||
|
||||
ok = (nlhs == 1) && ((nrhs == 1) || (nrhs == 2));
|
||||
|
||||
ok = ok && !mxIsEmpty(prhs[0]);
|
||||
ok = ok && mxIsStruct(prhs[0]);
|
||||
|
||||
ok = ok && (mxGetFieldNumber(prhs[0], "cartDims") >= 0);
|
||||
ok = ok && (mxGetFieldNumber(prhs[0], "COORD" ) >= 0);
|
||||
ok = ok && (mxGetFieldNumber(prhs[0], "ZCORN" ) >= 0);
|
||||
ok = ok && (mxGetFieldNumber(prhs[0], "ACTNUM" ) >= 0);
|
||||
|
||||
if (nrhs == 2) {
|
||||
ok = mxIsDouble(prhs[1]) && (mxGetNumberOfElements(prhs[1]) == 1);
|
||||
}
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static double
|
||||
define_tolerance(int nrhs, const mxArray *prhs[])
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
double tol;
|
||||
|
||||
tol = 0.0;
|
||||
|
||||
if (nrhs == 2) {
|
||||
tol = mxGetScalar(prhs[1]);
|
||||
}
|
||||
|
||||
return tol;
|
||||
}
|
||||
|
||||
|
||||
/* G = processgrid(grdecl)
|
||||
G = processgrid(grdecl, tolerance)
|
||||
*/
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
mexFunction(int nlhs, mxArray *plhs[],
|
||||
int nrhs, const mxArray *prhs[])
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
double tolerance;
|
||||
char errmsg[1023 + 1];
|
||||
struct grdecl grdecl;
|
||||
struct processed_grid g;
|
||||
|
||||
if (args_ok(nlhs, nrhs, prhs)) {
|
||||
mx_init_grdecl(&grdecl, prhs[0]);
|
||||
tolerance = define_tolerance(nrhs, prhs);
|
||||
|
||||
process_grdecl(&grdecl, tolerance, &g);
|
||||
|
||||
plhs[0] = allocate_grid(&g);
|
||||
|
||||
if (plhs[0] != NULL) {
|
||||
fill_grid(plhs[0], &g);
|
||||
} else {
|
||||
/* Failed to create grid structure. Return empty. */
|
||||
plhs[0] = mxCreateDoubleMatrix(0, 0, mxREAL);
|
||||
}
|
||||
|
||||
free_processed_grid(&g);
|
||||
} else {
|
||||
sprintf(errmsg,
|
||||
"Calling sequence is\n\t"
|
||||
"G = %s(grdecl)\t%%or\n\t"
|
||||
"G = %s(grdecl, tolerance)\n"
|
||||
"The 'grdecl' must be a valid structure with fields\n"
|
||||
"\t'cartDims', 'COORD', 'ZCORN', and 'ACTNUM'",
|
||||
mexFunctionName(), mexFunctionName());
|
||||
mexErrMsgTxt(errmsg);
|
||||
}
|
||||
}
|
||||
|
||||
#else
|
||||
|
||||
void fill_grid(mxArray **out, struct processed_grid *grid)
|
||||
{
|
||||
const char *names[] = {"nodes", "faces", "cells", "cartDims", "type"};
|
||||
@ -247,6 +653,7 @@ void fill_grid(mxArray **out, struct processed_grid *grid)
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* Gateway routine for Matlab mex function. */
|
||||
/*-------------------------------------------------------*/
|
||||
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
||||
@ -275,3 +682,4 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
||||
/* Free whatever was allocated in initGrdecl. */
|
||||
free_processed_grid(&out);
|
||||
}
|
||||
#endif
|
||||
|
||||
Loading…
Reference in New Issue
Block a user