Add preliminary support for wells in 'mex_ip_simple'. This is
untested, so disable direct MEX building for the time being.
This commit is contained in:
Bård Skaflestad
parent
caea28b2fb
commit
fc082bcc9a
327
mex_ip_simple.c
327
mex_ip_simple.c
@ -1,3 +1,4 @@
|
||||
#include <assert.h>
|
||||
#include <string.h>
|
||||
#include <mex.h>
|
||||
|
||||
@ -46,24 +47,26 @@ static int
|
||||
verify_grid_structure(const mxArray *G)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int nodes_ok, faces_ok, cells_ok, field_no;
|
||||
int nodes_ok = 0, faces_ok = 0, cells_ok = 0, field_no;
|
||||
|
||||
mxArray *pm;
|
||||
|
||||
nodes_ok = mxGetFieldNumber(G, "nodes") >= 0;
|
||||
if (mxIsStruct(G)) {
|
||||
nodes_ok = mxGetFieldNumber(G, "nodes") >= 0;
|
||||
|
||||
field_no = mxGetFieldNumber(G, "faces");
|
||||
faces_ok = field_no >= 0;
|
||||
if (faces_ok) {
|
||||
pm = mxGetFieldByNumber(G, 0, field_no);
|
||||
faces_ok = verify_faces_structure(pm);
|
||||
}
|
||||
field_no = mxGetFieldNumber(G, "faces");
|
||||
faces_ok = field_no >= 0;
|
||||
if (faces_ok) {
|
||||
pm = mxGetFieldByNumber(G, 0, field_no);
|
||||
faces_ok = verify_faces_structure(pm);
|
||||
}
|
||||
|
||||
field_no = mxGetFieldNumber(G, "cells");
|
||||
cells_ok = field_no >= 0;
|
||||
if (cells_ok) {
|
||||
pm = mxGetFieldByNumber(G, 0, field_no);
|
||||
cells_ok = verify_cells_structure(pm);
|
||||
field_no = mxGetFieldNumber(G, "cells");
|
||||
cells_ok = field_no >= 0;
|
||||
if (cells_ok) {
|
||||
pm = mxGetFieldByNumber(G, 0, field_no);
|
||||
cells_ok = verify_cells_structure(pm);
|
||||
}
|
||||
}
|
||||
|
||||
return nodes_ok && faces_ok && cells_ok;
|
||||
@ -72,35 +75,46 @@ verify_grid_structure(const mxArray *G)
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static int
|
||||
verify_structural_consistency(int nlhs, int nrhs, const mxArray *prhs[])
|
||||
verify_well_structure(const mxArray *W)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int rock_ok;
|
||||
int grid_ok;
|
||||
int conn_ok;
|
||||
int well_ok = mxIsStruct(W);
|
||||
|
||||
mxAssert((nlhs == 1) && (nrhs == 4),
|
||||
"Must be called with exactly four inputs and one output.");
|
||||
if (well_ok && (mxGetNumberOfElements(W) > 0)) {
|
||||
/* All STRUCT array elements have the same fields... */
|
||||
well_ok = (mxGetFieldNumber(W, "cells") >= 0);
|
||||
well_ok = well_ok && (mxGetFieldNumber(W, "WI" ) >= 0);
|
||||
}
|
||||
|
||||
mxAssert(mxIsStruct(prhs[0]) && mxIsStruct(prhs[1]),
|
||||
"First two inputs must be STRUCTs.");
|
||||
return well_ok;
|
||||
}
|
||||
|
||||
mxAssert((mxGetNumberOfElements(prhs[0]) == 1) &&
|
||||
(mxGetNumberOfElements(prhs[1]) == 1),
|
||||
"First two inputs must be 1-by-1 STRUCT arrays.");
|
||||
|
||||
grid_ok = verify_grid_structure(prhs[0]);
|
||||
/* ------------------------------------------------------------------ */
|
||||
static int
|
||||
args_ok(int nlhs, int nrhs, const mxArray *prhs[])
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int grid_ok = 0, rock_ok = 0, well_ok;
|
||||
|
||||
/* rock must contain a field 'perm' */
|
||||
rock_ok = mxGetFieldNumber(prhs[1], "perm") >= 0;
|
||||
if (nrhs > 0) {
|
||||
grid_ok = verify_grid_structure(prhs[0]);
|
||||
}
|
||||
|
||||
/* connection structure */
|
||||
conn_ok = (mxIsDouble(prhs[2]) || mxIsInt32(prhs[2])) &&
|
||||
(mxIsDouble(prhs[3]) || mxIsInt32(prhs[3])) &&
|
||||
(mxGetNumberOfElements(prhs[2]) <
|
||||
mxGetNumberOfElements(prhs[3]));
|
||||
if (nrhs > 1) {
|
||||
/* rock must contain a field 'perm' */
|
||||
rock_ok = mxIsStruct(prhs[1]) &&
|
||||
(mxGetFieldNumber(prhs[1], "perm") >= 0);
|
||||
}
|
||||
|
||||
return grid_ok && rock_ok && conn_ok;
|
||||
/* well structure */
|
||||
if (nrhs == 2) {
|
||||
well_ok = 1;
|
||||
} else {
|
||||
well_ok = verify_well_structure(prhs[2]);
|
||||
}
|
||||
|
||||
return (nlhs == 3) && grid_ok && rock_ok && well_ok;
|
||||
}
|
||||
|
||||
|
||||
@ -191,48 +205,163 @@ extract_connection_data(const mxArray *M_nconn, const mxArray *M_conn,
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static int
|
||||
extract_cell_data(const mxArray *G,
|
||||
const mxArray *M_nconn, const mxArray *M_conn,
|
||||
int *max_ncf, int *sum_nconn, int *sum_nconn2,
|
||||
int **ncfaces, int **nconn, int **conn,
|
||||
double **ccentroids, double **cvolumes)
|
||||
static void
|
||||
extract_grid_connections(const mxArray *G, int *nc, int *nf,
|
||||
int **ncf, int **pconn, int **conn)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int ncells, i, n;
|
||||
int i;
|
||||
|
||||
ncells = getNumberOfCells(G);
|
||||
*ncfaces = getCellFacePos(G);
|
||||
*nc = getNumberOfCells(G);
|
||||
*nf = getNumberOfFaces(G);
|
||||
|
||||
extract_connection_data(M_nconn, M_conn, ncells,
|
||||
nconn, conn, sum_nconn, sum_nconn2);
|
||||
/* Yes, two calls to getCellFacePos() */
|
||||
*ncf = getCellFacePos(G);
|
||||
*pconn = getCellFacePos(G);
|
||||
|
||||
*max_ncf = 0;
|
||||
*conn = getCellFaces(G);
|
||||
|
||||
for(i=0; i<ncells; ++i)
|
||||
{
|
||||
n = (*ncfaces)[i+1] - (*ncfaces)[i];
|
||||
(*ncfaces)[i] = n;
|
||||
|
||||
*max_ncf = MAX(*max_ncf, n);
|
||||
}
|
||||
|
||||
*ccentroids = getCellCentroids(G);
|
||||
*cvolumes = getCellVolumes(G);
|
||||
|
||||
return ncells;
|
||||
for (i = 0; i < *nc; i++) {
|
||||
(*ncf)[i] = (*ncf)[i + 1] - (*ncf)[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static void
|
||||
deallocate_cell_data(int *ncfaces, int *nconn, int *conn, double *ccentroids)
|
||||
add_well_connections(const mxArray *W, int nc, int nf,
|
||||
int *nw, int *pconn, int **conn)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int c, i, np, fld, neconn, w, p, n, dst, src;
|
||||
int *pi, *cwork, *c2w, *w2c, *tmp;
|
||||
double *pd;
|
||||
|
||||
mxArray *M_perfs;
|
||||
|
||||
if ((*nw = mxGetNumberOfElements(W)) > 0) {
|
||||
fld = mxGetFieldNumber(W, "cells");
|
||||
}
|
||||
|
||||
neconn = 0; /* Number of additional/extra conns */
|
||||
for (w = 0; w < *nw; w++) {
|
||||
neconn += mxGetNumberOfElements(mxGetFieldByNumber(W, w, fld));
|
||||
}
|
||||
|
||||
if (neconn > 0) {
|
||||
cwork = mxCalloc(nc , sizeof *cwork); /* all(cwork == 0) */
|
||||
c2w = mxMalloc(neconn * sizeof *c2w );
|
||||
w2c = mxMalloc(neconn * sizeof *w2c );
|
||||
|
||||
i = 0;
|
||||
for (w = 0; w < *nw; w++) {
|
||||
M_perfs = mxGetFieldByNumber(W, w, fld);
|
||||
np = mxGetNumberOfElements(M_perfs);
|
||||
|
||||
if (mxIsInt32(M_perfs)) {
|
||||
pi = mxGetData(M_perfs);
|
||||
|
||||
for (p = 0; p < np; p++, i++) {
|
||||
c = pi[p] - 1;
|
||||
|
||||
w2c[i] = c;
|
||||
c2w[i] = w;
|
||||
|
||||
cwork[c]++;
|
||||
}
|
||||
} else {
|
||||
pd = mxGetPr(M_perfs);
|
||||
|
||||
for (p = 0; p < np; p++, i++) {
|
||||
c = pd[p] - 1;
|
||||
|
||||
w2c[i] = c;
|
||||
c2w[i] = w;
|
||||
|
||||
cwork[c]++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
tmp = mxRealloc(*conn, (pconn[nc] + neconn) * sizeof *tmp);
|
||||
if (tmp != NULL) {
|
||||
n = pconn[nc] - pconn[nc - 1]; /* # exisiting conns */
|
||||
pconn[nc] += neconn; /* ubnd, *LAST* bin */
|
||||
|
||||
for (c = nc - 1; c > 0; c--) {
|
||||
dst = pconn[c + 1] - (n + cwork[c]);
|
||||
src = pconn[c + 0];
|
||||
|
||||
/* Move existing connections to start of new bin.
|
||||
* Note: Bins may overlap so use memmove(). */
|
||||
memmove(tmp + dst, tmp + src, n * sizeof *tmp);
|
||||
|
||||
/* Set cell's pointer for new connections. */
|
||||
cwork[c] = dst + n;
|
||||
|
||||
n = pconn[c] - pconn[c - 1];
|
||||
pconn[c] = dst;
|
||||
}
|
||||
|
||||
|
||||
/* Assign well connections.
|
||||
* Note: Generally, new DOF is nf + # ext. wells + c2w[i].
|
||||
* In this case, though, # ext. wells == 0 ... */
|
||||
for (i = 0; i < neconn; i++) {
|
||||
assert (cwork[w2c[i]] < pconn[w2c[i] + 1]);
|
||||
tmp[cwork[w2c[i]]++] = nf + c2w[i];
|
||||
}
|
||||
|
||||
/* Finally, assign updated connection structure. */
|
||||
*conn = tmp;
|
||||
}
|
||||
|
||||
mxFree(w2c); mxFree(c2w); mxFree(cwork);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static void
|
||||
count_connections(int nc, const int *pconn,
|
||||
int *max_nconn, int *sum_nconn2)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int c, nconn;
|
||||
|
||||
*max_nconn = *sum_nconn2 = 0;
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
nconn = pconn[c + 1] - pconn[c];
|
||||
|
||||
*max_nconn = MAX(*max_nconn, nconn);
|
||||
*sum_nconn2 += nconn * nconn;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static int
|
||||
extract_cell_data(const mxArray *G,
|
||||
double **ccentroids, double **cvolumes)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
*ccentroids = getCellCentroids(G);
|
||||
*cvolumes = getCellVolumes(G);
|
||||
|
||||
return getNumberOfCells(G);
|
||||
}
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static void
|
||||
deallocate_cell_data(int *ncf, int *pconn, int *conn, double *ccentroids)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
mxFree(ccentroids);
|
||||
mxFree(conn);
|
||||
mxFree(nconn);
|
||||
mxFree(ncfaces);
|
||||
mxFree(pconn);
|
||||
mxFree(ncf);
|
||||
}
|
||||
|
||||
|
||||
@ -245,8 +374,22 @@ deallocate_perm_data(double *K)
|
||||
}
|
||||
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
static void
|
||||
set_connections_M(int nc, const int *pconn, const int *conn,
|
||||
int *pconn_M, int *conn_M)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0; i <= nc ; i++) { pconn_M[i] = pconn[i] + 1; }
|
||||
for (i = 0; i < pconn[nc]; i++) { conn_M [i] = conn [i] + 1; }
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* BI = mex_ip_simple(G, rock, nconn, conn)
|
||||
* [BI, connPos, conns] = mex_ip_simple(G, rock)
|
||||
* [BI, connPos, conns] = mex_ip_simple(G, rock, W)
|
||||
*/
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
@ -255,6 +398,64 @@ mexFunction(int nlhs, mxArray *plhs[],
|
||||
int nrhs, const mxArray *prhs[])
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int d, nc, nf, nw = 0, max_nconn, sum_nconn2;
|
||||
int *pconn, *conn, *ncf, *fneighbour;
|
||||
double *fcentroid, *fnormal, *farea, *ccentroids, *cvolumes, *perm;
|
||||
char errmsg[1023 + 1];
|
||||
|
||||
double *Binv;
|
||||
|
||||
if (args_ok(nlhs, nrhs, prhs)) {
|
||||
extract_grid_connections(prhs[0], &nc, &nf, &ncf, &pconn, &conn);
|
||||
|
||||
if (nrhs == 3) {
|
||||
/* [...] = mex_ip_simple(G, rock, W) */
|
||||
add_well_connections(prhs[2], nc, nf, &nw, pconn, &conn);
|
||||
} else {
|
||||
nw = 0;
|
||||
}
|
||||
|
||||
count_connections(nc, pconn, &max_nconn, &sum_nconn2);
|
||||
|
||||
extract_face_data(prhs[0], &fneighbour, &farea,
|
||||
&fnormal, &fcentroid);
|
||||
|
||||
extract_cell_data(prhs[0], &ccentroids, &cvolumes);
|
||||
|
||||
d = getNumberOfDimensions(prhs[0]);
|
||||
perm = getPermeability(mxGetField(prhs[1], 0, "perm"), d);
|
||||
|
||||
/* Create return values */
|
||||
plhs[0] = mxCreateDoubleMatrix (sum_nconn2, 1, mxREAL);
|
||||
plhs[1] = mxCreateNumericMatrix(nc + 1 , 1, mxINT32_CLASS, mxREAL);
|
||||
plhs[2] = mxCreateNumericMatrix(pconn[nc] , 1, mxINT32_CLASS, mxREAL);
|
||||
|
||||
/* Compute IP for all cells (reservoir) */
|
||||
Binv = mxGetPr(plhs[0]); /* mxCreateDoubleMatrix == ZEROS */
|
||||
|
||||
/* mim_ip_simple_all() requires zeroed target array... */
|
||||
mim_ip_simple_all(nc, d, max_nconn, ncf, pconn, conn,
|
||||
fneighbour, fcentroid, fnormal, farea,
|
||||
ccentroids, cvolumes, perm, Binv);
|
||||
|
||||
set_connections_M(nc, pconn, conn,
|
||||
mxGetData(plhs[1]),
|
||||
mxGetData(plhs[2]));
|
||||
|
||||
deallocate_perm_data(perm);
|
||||
deallocate_cell_data(ncf, pconn, conn, ccentroids);
|
||||
deallocate_face_data(fneighbour, fnormal, fcentroid);
|
||||
} else {
|
||||
sprintf(errmsg,
|
||||
"Calling sequence is\n"
|
||||
"\t[BI, connPos, conn] = %s(G, rock)\nor\n"
|
||||
"\t[BI, connPos, conn] = %s(G, rock, W)\n",
|
||||
mexFunctionName(), mexFunctionName());
|
||||
|
||||
mexErrMsgTxt(errmsg);
|
||||
}
|
||||
}
|
||||
#if 0
|
||||
int ncells, nfaces, max_ncf, sum_ncf, sum_ncf2, d, structure_ok;
|
||||
|
||||
int *fneighbour;
|
||||
@ -269,6 +470,8 @@ mexFunction(int nlhs, mxArray *plhs[],
|
||||
structure_ok = verify_structural_consistency(nlhs, nrhs, prhs);
|
||||
|
||||
if (structure_ok) {
|
||||
extract_grid_connections(
|
||||
|
||||
d = getNumberOfDimensions(prhs[0]);
|
||||
|
||||
nfaces = extract_face_data(prhs[0], &fneighbour, &farea,
|
||||
@ -298,4 +501,4 @@ mexFunction(int nlhs, mxArray *plhs[],
|
||||
} else {
|
||||
plhs[0] = mxCreateDoubleScalar(mxGetNaN());
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -65,6 +65,8 @@ function varargout = mex_ip_simple(varargin)
|
||||
% $Date$
|
||||
% $Revision$
|
||||
|
||||
error('Currently being re-worked');
|
||||
|
||||
buildmex CFLAGS="\$CFLAGS -Wall -Wextra -ansi -pedantic ...
|
||||
-Wformat-nonliteral -Wcast-align -Wpointer-arith ...
|
||||
-Wbad-function-cast -Wmissing-prototypes -Wstrict-prototypes ...
|
||||
|
||||
27
mimetic.c
27
mimetic.c
@ -16,27 +16,27 @@
|
||||
|
||||
/* ------------------------------------------------------------------ */
|
||||
void
|
||||
mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
|
||||
int *nconn, int *cf,
|
||||
mim_ip_simple_all(int ncells, int d, int max_nconn, int *ncf,
|
||||
int *pconn, int *conn,
|
||||
int *fneighbour, double *fcentroid, double *fnormal,
|
||||
double *farea, double *ccentroid, double *cvol,
|
||||
double *perm, double *Binv)
|
||||
/* ------------------------------------------------------------------ */
|
||||
{
|
||||
int i, j, c, f, nf, fpos, fpos2, lwork;
|
||||
int i, j, c, f, nf, nconn, fpos2, lwork;
|
||||
|
||||
double *C, *N, *A, *work, s;
|
||||
|
||||
double cc[3] = { 0.0 }; /* No more than 3 space dimensions */
|
||||
|
||||
lwork = 64 * (max_ncf * d); /* 64 from ILAENV() */
|
||||
C = malloc((max_ncf * d) * sizeof *C);
|
||||
N = malloc((max_ncf * d) * sizeof *N);
|
||||
A = malloc(max_ncf * sizeof *A);
|
||||
work = malloc(lwork * sizeof *work);
|
||||
lwork = 64 * (max_nconn * d); /* 64 from ILAENV() */
|
||||
C = malloc((max_nconn * d) * sizeof *C);
|
||||
N = malloc((max_nconn * d) * sizeof *N);
|
||||
A = malloc(max_nconn * sizeof *A);
|
||||
work = malloc(lwork * sizeof *work);
|
||||
|
||||
if ((C != NULL) && (N != NULL) && (A != NULL) && (work != NULL)) {
|
||||
fpos = fpos2 = 0;
|
||||
fpos2 = 0;
|
||||
|
||||
for (c = 0; c < ncells; c++) {
|
||||
for (j = 0; j < d; j++) {
|
||||
@ -46,7 +46,7 @@ mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
|
||||
nf = ncf[c];
|
||||
|
||||
for (i = 0; i < nf; i++) {
|
||||
f = cf[fpos + i];
|
||||
f = conn[pconn[c] + i];
|
||||
s = 2.0*(fneighbour[2 * f] == c) - 1.0;
|
||||
|
||||
A[i] = farea[f];
|
||||
@ -57,11 +57,12 @@ mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
|
||||
}
|
||||
}
|
||||
|
||||
mim_ip_simple(nf, nconn[c], d, cvol[c], &perm[c * d * d],
|
||||
nconn = pconn[c + 1] - pconn[c];
|
||||
|
||||
mim_ip_simple(nf, nconn, d, cvol[c], &perm[c * d * d],
|
||||
C, A, N, &Binv[fpos2], work, lwork);
|
||||
|
||||
fpos += nconn[c];
|
||||
fpos2 += nconn[c] * nconn[c];
|
||||
fpos2 += nconn * nconn;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -21,7 +21,7 @@ void mim_ip_simple(int nf, int nconn, int d,
|
||||
double *work, int lwork);
|
||||
|
||||
void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
|
||||
int *nconn, int *conn,
|
||||
int *pconn, int *conn,
|
||||
int *fneighbour, double *fcentroid, double *fnormal,
|
||||
double *farea, double *ccentroid, double *cvol,
|
||||
double *perm, double *Binv);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user