By jrn's blessing, enable the alternative implementation of 'processgrid'. Also remove the existing implementation.

While here, assert copyright for 2011.

Signed-off-by: Bård Skaflestad <Bard.Skaflestad@sintef.no>
This commit is contained in:
Bård Skaflestad 2011-01-09 15:15:43 +00:00 committed by Bård Skaflestad
parent 714e052822
commit 691c3a3b7b

View File

@ -13,8 +13,8 @@
//===========================================================================*/ //===========================================================================*/
/* /*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. Copyright 2009, 2010, 2011 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA. Copyright 2009, 2010, 2011 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS). This file is part of The Open Reservoir Simulator Project (OpenRS).
@ -46,7 +46,6 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
#include "mxgrdecl.h" #include "mxgrdecl.h"
#if 0
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
static mxArray * static mxArray *
allocate_nodes(size_t nnodes) allocate_nodes(size_t nnodes)
@ -454,237 +453,3 @@ mexFunction(int nlhs, mxArray *plhs[],
mexErrMsgTxt(errmsg); mexErrMsgTxt(errmsg);
} }
} }
#else
void fill_grid(mxArray **out, struct processed_grid *grid)
{
const char *names[] = {"nodes", "faces", "cells", "cartDims", "type"};
mxArray *G = mxCreateStructMatrix(1,1,sizeof names / sizeof names[0],names);
int i,j;
double *ptr;
/* nodes */
const char *n2[] = {"num", "coords"};
mxArray *nodes = mxCreateStructMatrix(1,1,sizeof n2 / sizeof n2[0],n2);
mxSetField(nodes, 0, "num", mxCreateDoubleScalar(grid->number_of_nodes));
mxArray *nodecoords = mxCreateDoubleMatrix(grid->number_of_nodes, 3, mxREAL);
ptr = mxGetPr(nodecoords);
for (j=0;j<3;++j)
{
for(i=0; i<grid->number_of_nodes; ++i)
{
ptr[i+grid->number_of_nodes*j] = grid->node_coordinates[3*i+j];
}
}
mxSetField(nodes, 0, "coords", nodecoords);
mxSetField(G, 0, "nodes", nodes);
/* faces */
const char *n3[] = {"num", "neighbors", "nodes", "numNodes", "nodePos", "tag"};
mxArray *faces = mxCreateStructMatrix(1,1,sizeof n3 / sizeof n3[0], n3);
mxSetField(faces, 0, "num", mxCreateDoubleScalar(grid->number_of_faces));
mxArray *faceneighbors = mxCreateNumericMatrix(grid->number_of_faces, 2,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(faceneighbors);
for(j=0; j<2; ++j)
{
for (i=0; i<grid->number_of_faces; ++i)
{
int ix = grid->face_neighbors[2*i+j];
if (ix == -1)
{
iptr[i+grid->number_of_faces*j] = 0;
}
else
{
iptr[i+grid->number_of_faces*j] = ix+1;
}
}
}
}
mxSetField(faces, 0, "neighbors", faceneighbors);
mxArray *numnodes = mxCreateNumericMatrix(grid->number_of_faces, 1,
mxINT32_CLASS, mxREAL);
mxArray *nodepos = mxCreateNumericMatrix(grid->number_of_faces+1, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr1 = mxGetData(numnodes);
int *iptr2 = mxGetData(nodepos);
iptr2[0] = 1;
for (i=0; i<grid->number_of_faces; ++i)
{
iptr1[i] = grid->face_ptr[i+1] - grid->face_ptr[i];
iptr2[i+1] = iptr2[i] + iptr1[i];
}
}
mxSetField(faces, 0, "numNodes", numnodes);
mxSetField(faces, 0, "nodePos", nodepos);
mxArray *tags = mxCreateNumericMatrix(grid->number_of_faces, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(tags);
for (i = 0; i < grid->number_of_faces; ++i)
{
iptr[i] = grid->face_tag[i] + 1;
}
}
mxSetField(faces, 0, "tag", tags);
const char *n4[] = {"num", "faces", "facePos", "indexMap"};
mxArray *cells = mxCreateStructMatrix(1,1,sizeof n4 / sizeof n4[0], n4);
mxSetField(cells, 0, "num", mxCreateDoubleScalar(grid->number_of_cells));
mxArray *map = mxCreateNumericMatrix(grid->number_of_cells, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(map);
for(i=0; i<grid->number_of_cells; ++i)
{
iptr[i] = grid->local_cell_index[i]+1;
}
}
mxSetField(cells, 0, "indexMap", map);
mxArray *facepos = mxCreateNumericMatrix(grid->number_of_cells+1, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(facepos);
for(i=0; i<grid->number_of_cells+1; ++i)
{
iptr[i] = 0;
}
for (i=0; i<2*grid->number_of_faces; ++i)
{
int c=grid->face_neighbors[i];
if(c != -1)
{
iptr[c+1]++;
}
}
iptr[0] = 1;
for(i=0; i<grid->number_of_cells; ++i)
{
iptr[i+1] += iptr[i];
}
}
mxSetField(cells, 0, "facePos", facepos);
int *counter = calloc(grid->number_of_cells, sizeof(*counter));
int num_half_faces = 0;
{
int *iptr = mxGetData(facepos);
for(i=0; i<grid->number_of_cells; ++i)
{
counter[i] = num_half_faces;
num_half_faces += iptr[i+1]-iptr[i];
}
}
mxArray *cellfaces = mxCreateNumericMatrix(num_half_faces, 2,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(cellfaces);
int c1, c2, cf_tag;
for (i=0; i<grid->number_of_faces; ++i)
{
cf_tag = 2 * grid->face_tag[i] + 1; /* [1, 3, 5] */
c1 = grid->face_neighbors[2*i+0];
c2 = grid->face_neighbors[2*i+1];
if(c1 != -1) {
iptr[counter[c1] + 0*num_half_faces] = i+1;
iptr[counter[c1] + 1*num_half_faces] = cf_tag + 1; /* out */
counter[c1] += 1;
}
if(c2 != -1) {
iptr[counter[c2] + 0*num_half_faces] = i+1;
iptr[counter[c2] + 1*num_half_faces] = cf_tag + 0; /* in */
counter[c2] += 1;
}
}
}
mxSetField(cells, 0, "faces", cellfaces);
mxSetField(G, 0, "cells", cells);
int n = grid->face_ptr[grid->number_of_faces];
mxArray *facenodes = mxCreateNumericMatrix(n, 1,
mxINT32_CLASS, mxREAL);
{
int *iptr = mxGetData(facenodes);
for (i=0; i<n; ++i)
{
iptr[i] = grid->face_nodes[i]+1;
}
}
mxSetField(faces, 0, "nodes", facenodes);
mxSetField(G, 0, "faces", faces);
free(counter);
/* cartDims */
mxArray *cartDims = mxCreateDoubleMatrix(1, 3, mxREAL);
ptr = mxGetPr(cartDims);
ptr[0] = grid->dimensions[0];
ptr[1] = grid->dimensions[1];
ptr[2] = grid->dimensions[2];
mxSetField(G, 0, "cartDims", cartDims);
/* type */
mxArray *type = mxCreateCellMatrix(1, 1);
mxSetCell(type, 0, mxCreateString("processgrid"));
mxSetField(G, 0, "type", type);
out[0] = G;
}
/* Gateway routine for Matlab mex function. */
/*-------------------------------------------------------*/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/* Set up data passed from Matlab */
struct grdecl g;
struct processed_grid out;
double tolerance = 0.0;
mx_init_grdecl(&g, prhs[0]);
if (nrhs == 2)
{
tolerance = mxGetScalar (prhs[1]);
}
process_grdecl(&g, tolerance, &out);
if (plhs >0)
{
/* write to matlab struct */
fill_grid(plhs, &out);
}
/* Free whatever was allocated in initGrdecl. */
free_processed_grid(&out);
}
#endif