OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Return pointers from getXX(...). Transpose centroids and normals.

Browse Source 
And coordinates.
This commit is contained in:
Jostein R. Natvig
2010-08-05 12:13:59 +00:00
parent 53de3c03d1
commit 97565f3bb7
2 changed files with 194 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files

226
mrst_api.c
View File

@@ -1,10 +1,7 @@
/* "API" */
#include <stdlib.h>
#include <mex.h>
#include "mrst_api.h"
#define Assert mxAssert
/* ------------------------------------------------------------------ */
static mxArray*
getField(const mxArray *a, char *field, char *subfield)
@@ -25,8 +22,8 @@ static int *
extractIntMatrix(const mxArray *a)
/* ------------------------------------------------------------------ */
{
int n = mxGetNumberOfElements(a);
int *q = malloc(n * sizeof *q);
int n = mxGetNumberOfElements(a);
int *q = mxMalloc(n * sizeof *q);
if (q != NULL)
{
if (mxIsInt32(a))
@@ -35,8 +32,6 @@ extractIntMatrix(const mxArray *a)
int i,j;
for (i=0; i<n; ++i)
{
mxAssert (p[i] <= INT_MAX,
"Matrix entry exceeds INT_MAX");
q[i] = p[i]-1;
}
}
@@ -46,7 +41,7 @@ extractIntMatrix(const mxArray *a)
int i,j;
for (i=0; i<n; ++i)
{
mxAssert (p[i] <= INT_MAX,
mxAssert ((1 <= p[i]) && (p[i] <= INT_MAX),
"Matrix entry exceeds INT_MAX");
q[i] = p[i]-1;
}
@@ -61,48 +56,80 @@ static int *
extractIntMatrixTranspose(const mxArray *a)
/* ------------------------------------------------------------------ */
{
int *q = extractIntMatrix(a);
int *p = extractIntMatrix(a);
int M = mxGetM(a);
int N = mxGetN(a);
int i,j;
for(i=0; i<M; ++i)
int *q = mxMalloc(M * N * sizeof *q);
if (q != NULL)
{
for(j=0; j<i; ++j)
int i,j;
for(i=0; i<M; ++i)
{
int tmp = q[i+M*j];
q[i+M*j] = q[j+N*i];
q[j+N*i] = tmp;
for(j=0; j<N; ++j)
{
q[i*N+j] = p[i+M*j];
}
}
}
mxFree(p);
return q;
}
/* ------------------------------------------------------------------ */
static double *
extractDoubleMatrixTranspose(const mxArray *a)
/* ------------------------------------------------------------------ */
{
int M = mxGetM(a);
int N = mxGetN(a);
double *q = mxMalloc(M * N * sizeof *q);
if (q != NULL)
{
double *p = mxGetPr(a);
int i,j;
for(i=0; i<M; ++i)
{
for(j=0; j<N; ++j)
{
q[i*N+j] = p[i+M*j];
}
}
}
return q;
}
/* ------------------------------------------------------------------ */
int getNumberOfDimensions(const mxArray *G)
int
getNumberOfDimensions(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return mxGetN(getField(G, "nodes", "coords"));
}
/* ------------------------------------------------------------------ */
void getLocal2GlobalCellMap(const mxArray *G)
void
getLocal2GlobalCellMap(const mxArray *G)
/* ------------------------------------------------------------------ */
{
Assert(0, "Not implemented!");
mxAssert(0, "Not implemented!");
}
/* ------------------------------------------------------------------ */
int getNumberOfNodes(const mxArray *G)
int
getNumberOfNodes(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return mxGetM(getField(G, "nodes", "coords"));
}
/* ------------------------------------------------------------------ */
double *getNodeCoordinates(const mxArray *G)
double *
getNodeCoordinates(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return extractDoubleMatrixTranspose(getField(G, "nodes", "coords"));
#if 0
mxArray *p1, *p2;
p1 = mxGetField(G , 0, "nodes" );
@@ -111,7 +138,7 @@ double *getNodeCoordinates(const mxArray *G)
const int n = getNumberOfNodes(G);
const int d = getNumberOfDimensions(G);
double *v = malloc(n * d * sizeof *v);
double *v = mxMalloc(n * d * sizeof *v);
if (v != NULL)
{
@@ -126,31 +153,36 @@ double *getNodeCoordinates(const mxArray *G)
}
}
return v;
#endif
}
/* ------------------------------------------------------------------ */
int getNumberOfFaces(const mxArray *G)
int
getNumberOfFaces(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return mxGetNumberOfElements(getField(G, "faces", "nodePos"))-1;
}
/* ------------------------------------------------------------------ */
int *getFaceNodePos(const mxArray *G)
int *
getFaceNodePos(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return extractIntMatrix(getField(G, "faces", "nodePos"));
}
/* ------------------------------------------------------------------ */
int getNumberOfFaceNodes(const mxArray *G)
int
getNumberOfFaceNodes(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return mxGetNumberOfElements(getField(G, "faces", "nodes"));
}
/* ------------------------------------------------------------------ */
int *getFaceNodes(const mxArray *G)
int *
getFaceNodes(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return extractIntMatrix(getField(G, "faces", "nodes"));
@@ -158,74 +190,184 @@ int *getFaceNodes(const mxArray *G)
/* ------------------------------------------------------------------ */
int *getFaceCellNeighbors(const mxArray *G)
int *
getFaceCellNeighbors(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return extractIntMatrixTranspose(getField(G, "faces", "neighbors"));
}
/* ------------------------------------------------------------------ */
void getFaceAreas(const mxArray *G, double **v)
double *
getFaceAreas(const mxArray *G)
/* ------------------------------------------------------------------ */
{
*v = mxGetPr(getField(G, "faces", "areas"));
return mxGetPr(getField(G, "faces", "areas"));
}
/* ------------------------------------------------------------------ */
void getFaceNormals(const mxArray *G, double **v)
double *
getFaceNormals(const mxArray *G)
/* ------------------------------------------------------------------ */
{
*v = mxGetPr(getField(G, "faces", "normals"));
return extractDoubleMatrixTranspose(getField(G, "faces", "normals"));
}
/* ------------------------------------------------------------------ */
void getFaceCentroids(const mxArray *G, double **v)
double *
getFaceCentroids(const mxArray *G)
/* ------------------------------------------------------------------ */
{
*v = mxGetPr(getField(G, "faces", "centroids"));
return extractDoubleMatrixTranspose(getField(G, "faces", "centroids"));
}
/* ------------------------------------------------------------------ */
int getNumberOfCells(const mxArray *G)
int
getNumberOfCells(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return mxGetNumberOfElements(getField(G, "cells", "facePos"))-1;
return mxGetNumberOfElements(getField(G, "cells", "facePos"))-1;
}
/* ------------------------------------------------------------------ */
int *getCellFacePos(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return extractIntMatrix(getField(G, "cells", "facePos"));
return extractIntMatrix(getField(G, "cells", "facePos"));
}
/* ------------------------------------------------------------------ */
int getNumberOfCellFaces(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return mxGetNumberOfElements(getField(G, "cells", "faces"))-1;
return mxGetNumberOfElements(getField(G, "cells", "faces"))-1;
}
/* ------------------------------------------------------------------ */
int *getCellFaces(const mxArray *G)
/* ------------------------------------------------------------------ */
{
return extractIntMatrix(getField(G, "cells", "faces"));
return extractIntMatrix(getField(G, "cells", "faces"));
}
/* ------------------------------------------------------------------ */
void getCellVolumes(const mxArray *G, double **v)
double *
getCellVolumes(const mxArray *G)
/* ------------------------------------------------------------------ */
{
*v = mxGetPr(getField(G, "cells", "volumes"));
return mxGetPr(getField(G, "cells", "volumes"));
}
/* ------------------------------------------------------------------ */
void getCellCentroids(const mxArray *G, double **v)
double *
getCellCentroids(const mxArray *G)
/* ------------------------------------------------------------------ */
{
*v = mxGetPr(getField(G, "cells", "centroids"));
return extractDoubleMatrixTranspose(getField(G, "cells", "centroids"));
}
/*
*
*
*
*/
/* ------------------------------------------------------------------ */
static int
allocate_perm_data(int ncells, int d, double **K)
/* ------------------------------------------------------------------ */
{
int ret;
double *perm;
perm = mxMalloc(ncells * d * d * sizeof *perm);
if (perm != NULL) {
*K = perm;
ret = 1;
} else {
ret = 0;
}
return ret;
}
/* ------------------------------------------------------------------ */
double *
getPermeability(const mxArray *perm, int d)
/* ------------------------------------------------------------------ */
{
int ncells, ncomp, alloc_ok;
int c, i, off;
double *k, *tensor;
ncells = mxGetM(perm);
ncomp = mxGetN(perm);
alloc_ok = allocate_perm_data(ncells, d, &k);
if (alloc_ok) {
for (i = 0; i < ncells * d * d; i++) {
k[i] = 0.0;
}
tensor = mxGetPr(perm);
if (ncomp == 1) {
/* Isotropic (scalar) tensor */
for (c = 0; c < ncells; c++) {
off = c * d * d;
for (i = 0; i < d; i++) {
k[i*(d + 1) + off] = tensor[c];
}
}
} else if (ncomp == d) {
/* Diagonal tensor */
for (c = 0; c < ncells; c++) {
off = c * d * d;
for (i = 0; i < d; i++) {
k[i*(d + 1) + off] = tensor[c + i*ncells];
}
}
} else if (d == 2) {
/* Full 2D tensor */
mxAssert (ncomp == 3, "");
for (c = 0; c < ncells; c++) {
off = c * d * d;
k[0 + off] = tensor[c + 0*ncells];
k[1 + off] = tensor[c + 1*ncells];
k[2 + off] = tensor[c + 1*ncells];
k[3 + off] = tensor[c + 2*ncells];
}
} else {
/* Full 3D tensor */
mxAssert ((d == 3) && (ncomp == 6), "");
for (c = 0; c < ncells; c++) {
off = c * d * d;
k[0 + off] = tensor[c + 0*ncells];
k[1 + off] = tensor[c + 1*ncells];
k[2 + off] = tensor[c + 2*ncells];
k[3 + off] = tensor[c + 1*ncells];
k[4 + off] = tensor[c + 3*ncells];
k[5 + off] = tensor[c + 4*ncells];
k[6 + off] = tensor[c + 2*ncells];
k[7 + off] = tensor[c + 4*ncells];
k[8 + off] = tensor[c + 5*ncells];
}
}
}
return k;
}

15
mrst_api.h
View File

@@ -20,9 +20,9 @@ int *getFaceNodes (const mxArray *G); /* copy */
int *getFaceCellNeighbors (const mxArray *G); /* copy */
/* Face geometry */
void getFaceAreas (const mxArray *G, double **v);
void getFaceNormals (const mxArray *G, double **v);
void getFaceCentroids (const mxArray *G, double **v);
double *getFaceAreas (const mxArray *G);
double *getFaceNormals (const mxArray *G); /* copy */
double *getFaceCentroids (const mxArray *G); /* copy */
/* Cell topology */
int getNumberOfCells (const mxArray *G);
@@ -31,7 +31,12 @@ int *getCellFacePos (const mxArray *G); /* copy */
int *getCellFaces (const mxArray *G); /* copy */
/* Cell geometry */
void getCellVolumes (const mxArray *G, double **v);
void getCellCentroids (const mxArray *G, double **v);
double *getCellVolumes (const mxArray *G);
double *getCellCentroids (const mxArray *G); /* copy */
/* Rock properties */
double *
getPermeability(const mxArray *perm, int d); /* copy */
#endif /* MRST_API_H_INCLUDED */