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

Add system assembly.

Browse Source
This commit is contained in:
Jostein R. Natvig
2010-08-06 10:53:07 +00:00
parent f7e93692f9
commit 7e3fe07012
4 changed files with 249 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files
hybsys.c
+184
View File
@@ -215,3 +215,187 @@ hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
p2 += nconn[c] * nconn[c];
}
}
/*
* Routines to assemble global matrix
*
*/
/* ---------------------------------------------------------------------- */
static int *
hybsys_build_ia(int nc, int nf, int *nconn, int *conn)
/* ---------------------------------------------------------------------- */
{
int *ia = malloc((nf+1) * sizeof *ia);
int i;
for(i=0; i<nf+1; ++i)
{
ia[i] = 0;
}
/*
* Compute rowsizes
*/
int c, pos = 0;
for(c=0; c<nc; ++c)
{
int n = nconn[c];
for (i=pos; i<pos+n; ++i)
{
mxAssert(conn[i]<nf, "conn out of bounds");
ia[1+conn[i]] += n - 1;
}
pos += n;
}
/*
* cumulative sum...
*/
for(i=1; i<nf+1; ++i)
{
ia[i] = ia[i-1] + ia[i]+1;
}
return ia;
}
/* ---------------------------------------------------------------------- */
static int*
hybsys_build_ja(int nc, int nf, int *nconn, int *conn,
int *ia, int *work)
/* ---------------------------------------------------------------------- */
{
int *ja = malloc(ia[nf] * sizeof *ja);
int i,j;
/*
* For each row, diagonal entries are positioned first.
*/
for(i=0; i<nf; ++i)
{
work[i] = 1;
ja[ia[i]] = i;
}
int c, pos = 0;
for(c=0; c<nc; ++c)
{
int n = nconn[c];
for (i=pos; i<pos+n; ++i)
{
int fi = conn[i];
/* mxAssert(fi<nf, "fi out of bounds"); */
for (j=pos; j<i; ++j)
{
int fj = conn[j];
/* mxAssert(fj<nf, "fj out of bounds"); */
/*
* No conditionals since off-diagonals entries are
* visited only once.
*/
ja[ia[fi] + work[fi]++] = fj;
ja[ia[fj] + work[fj]++] = fi;
}
}
pos += n;
}
return ja;
}
/* ---------------------------------------------------------------------- */
static double*
hybsys_build_sa(int nc, int nf, int *nconn, int *conn, int *ia,
double *S, double *R, int *work)
/* ---------------------------------------------------------------------- */
{
double *sa = malloc(ia[nf] * sizeof *sa);
int i,j;
/*
* Clear diagonal and work array
*/
for(i=0; i<nf; ++i)
{
work[i] = 1;
sa[ia[i]] = 0;
R[i] = 0;
}
double *s = S;
double *r = R;
int c, pos = 0;
for(c=0; c<nc; ++c)
{
int n = nconn[c];
for (i=pos; i<pos+n; ++i)
{
int fi = conn[i];
int ii = i-pos;
for (j=pos; j<i; ++j)
{
int fj = conn[j];
int jj = j-pos;
/*
* We can use assignment since off-diagonal entries are
* visited only once.
*/
sa[ia[fi] + work[fi]++] = s[ii + jj*n];
sa[ia[fj] + work[fj]++] = s[jj + ii*n];
}
/*
* Diagonal entries are more than once.
*/
sa[ia[fi]] += s[ii + ii*n];
R[fi] += r[ii];
}
s += n*n;
pos += n;
}
return sa;
}
/* ---------------------------------------------------------------------- */
void hybsys_build_matrix_structure(int nc, int nf, int *nconn, int *conn,
int **ia, int **ja)
/* ---------------------------------------------------------------------- */
{
int *work = malloc(nf * sizeof *work);
*ia = hybsys_build_ia(nc, nf, nconn, conn);
*ja = hybsys_build_ja(nc, nf, nconn, conn, *ia, work);
free(work);
}
/* ---------------------------------------------------------------------- */
void hybsys_assemble_global_system(int nc, int nf, int *nconn, int *conn,
double *S, double *R,
double **sa, int *ia)
/* ---------------------------------------------------------------------- */
{
int *work = malloc(nf * sizeof *work);
*sa = hybsys_build_sa(nc, nf, nconn, conn, ia, S, R, work);
free(work);
}
/* ---------------------------------------------------------------------- */
struct Sparse*
hybsys_assemble(int nc, int nf, int *nconn, int *conn, double *S, double *R)
/* ---------------------------------------------------------------------- */
{
struct Sparse *A = malloc(sizeof *A);
A->m = A->n = nf;
hybsys_build_matrix_structure(nc, nf, nconn, conn, &A->ia, &A->ja);
hybsys_assemble_global_system(nc, nf, nconn, conn, S, R, &A->sa, A->ia);
return A;
}
hybsys.h
+10 -1
View File
@@ -8,7 +8,14 @@ struct hybsys {
double *S; /* system matrix in single cell */
double *one; /* ones(max_ncf, 1) */
};
struct Sparse
{
int m;
int n;
int *ia;
int *ja;
double *sa;
};
struct hybsys *
hybsys_allocate(int max_ncf, int nc, int ncf_tot);
@@ -38,4 +45,6 @@ hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
const double *pi, double *press, double *flux,
double *work, const int lwork);
struct Sparse*
hybsys_assemble(int nc, int nf, int *nconn, int *conn, double *S, double *R);
#endif /* HYBSYS_H_INCLUDED */
mex_schur_comp_symm.c
+54 -1
View File
@@ -14,7 +14,7 @@ verify_args(int nlhs, int nrhs, const mxArray *prhs[])
{
int ok;
ok = (nlhs == 4) && (nrhs == 2);
ok = (nlhs == 4) && (nrhs == 3);
ok = ok && mxIsDouble(prhs[0]);
ok = ok && (mxIsDouble(prhs[1]) || mxIsInt32(prhs[1]));
ok = ok && (mxGetNumberOfElements(prhs[0]) >
@@ -124,6 +124,51 @@ get_nconn(const mxArray *M_nconn, int *nconn)
}
}
/* ---------------------------------------------------------------------- */
static void
get_conn(const mxArray *M_conn, int *conn)
/* ---------------------------------------------------------------------- */
{
size_t nel, i;
int *pi;
double *pd;
nel = mxGetNumberOfElements(M_conn);
if (mxIsDouble(M_conn)) {
pd = mxGetPr(M_conn);
for (i = 0; i < nel; i++) { conn[i] = pd[i] - 1; }
} else {
pi = mxGetData(M_conn);
for (i = 0; i < nel; i++) { conn[i] = pi[i] - 1; }
}
}
/* ---------------------------------------------------------------------- */
static int
get_number_of_faces(int nc, int *nconn, int *conn)
/* ---------------------------------------------------------------------- */
{
int N = 0;
int i;
for (i=0; i<nc; ++i)
{
N += nconn[i];
}
int nf=0;
for(i=0; i<N; ++i)
{
nf = MAX(nf, conn[i]+1);
}
return nf;
}
/*
* [S, r, F, L] = mex_schur_comp_symm(BI, nconn)
@@ -173,7 +218,14 @@ mexFunction(int nlhs, mxArray *plhs[],
p1 += nconn[c];
p2 += nconn[c] * nconn[c];
}
int *conn= mxMalloc(mxGetNumberOfElements(prhs[2])*sizeof *conn);
get_conn(prhs[2], conn);
int nf = get_number_of_faces(nc, nconn, conn);
struct Sparse *A = hybsys_assemble(nc, nf, nconn, conn, ptr, sys->r);
free(A->ia); free(A->ja); free(A->sa); free(A);
ptr = mxGetPr(plhs[1]);
memcpy(ptr, sys->r, ncf_tot * sizeof *ptr);
@@ -185,5 +237,6 @@ mexFunction(int nlhs, mxArray *plhs[],
hybsys_free(sys);
deallocate_aux_arrays(nconn, src, gflux);
mxFree(conn);
}
}
test_mex_schur_comp_symm.m
+1 -1
View File
@@ -5,7 +5,7 @@ BI = mex_ip_simple(G, rock);
nconn = diff(G.cells.facePos);
conn = G.cells.faces(:, 1);
[S, r, F, L] = mex_schur_comp_symm(BI, nconn);
[S, r, F, L] = mex_schur_comp_symm(BI, nconn, conn);
[i, j] = blockDiagIndex(nconn, nconn);