Add rudimentary MEX support for computing non-symmetric (i.e.,
compressible) cell contributions.
This commit is contained in:
Bård Skaflestad
parent
47abc777e7
commit
8df4a5f14b
207
mex_schur_comp.c
Normal file
207
mex_schur_comp.c
Normal file
@ -0,0 +1,207 @@
|
||||
#include <string.h>
|
||||
|
||||
#include <mex.h>
|
||||
#define MAT_SIZE_T mwSignedIndex
|
||||
|
||||
#include "hybsys.h"
|
||||
|
||||
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
verify_args(int nlhs, int nrhs, const mxArray *prhs[])
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ok;
|
||||
|
||||
ok = (nlhs == 5) && (nrhs == 3);
|
||||
ok = ok && mxIsDouble(prhs[0]);
|
||||
ok = ok && (mxIsDouble(prhs[1]) || mxIsInt32(prhs[1]));
|
||||
ok = ok && (mxGetNumberOfElements(prhs[0]) >
|
||||
mxGetNumberOfElements(prhs[1]));
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
count_cellconn(int nc, const int *pconn)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, nconn, max_nconn;
|
||||
|
||||
max_nconn = 0;
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
nconn = pconn[c + 1] - pconn[c];
|
||||
|
||||
max_nconn = MAX(max_nconn, nconn);
|
||||
}
|
||||
|
||||
return max_nconn;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
deallocate_aux_arrays(int *pconn, double *src, double *gflux,
|
||||
double *BIV, double *P)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
/* Apparently mxFree() makes no guarantee regarding NULL arguments.
|
||||
* Institute a belt-and-suspenders approach to releasing resources.
|
||||
*/
|
||||
if (P != NULL) { mxFree(P); }
|
||||
if (BIV != NULL) { mxFree(BIV); }
|
||||
if (gflux != NULL) { mxFree(gflux); }
|
||||
if (src != NULL) { mxFree(src); }
|
||||
if (pconn != NULL) { mxFree(pconn); }
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
allocate_aux_arrays(int nc, int nconn_tot,
|
||||
double **src, double **gflux,
|
||||
double **BIV, double **P)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ret;
|
||||
double *s, *g, *B, *p;
|
||||
|
||||
s = mxMalloc(nc * sizeof *s);
|
||||
g = mxMalloc(nconn_tot * sizeof *g);
|
||||
B = mxMalloc(nconn_tot * sizeof *B);
|
||||
p = mxMalloc(nc * sizeof *p);
|
||||
|
||||
if ((s == NULL) || (g == NULL) || (B == NULL) || (p == NULL)) {
|
||||
deallocate_aux_arrays(NULL, s, g, B, p);
|
||||
|
||||
*src = NULL;
|
||||
*gflux = NULL;
|
||||
*BIV = NULL;
|
||||
*P = NULL;
|
||||
|
||||
ret = 0;
|
||||
} else {
|
||||
*src = s;
|
||||
*gflux = g;
|
||||
*BIV = B;
|
||||
*P = p;
|
||||
|
||||
ret = 1;
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
get_pconn(const mxArray *M_pconn, int **pconn)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ret;
|
||||
|
||||
size_t e, ne;
|
||||
|
||||
int *pi;
|
||||
double *pd;
|
||||
|
||||
ne = mxGetNumberOfElements(M_pconn);
|
||||
|
||||
*pconn = mxMalloc(ne * sizeof **pconn);
|
||||
|
||||
if (*pconn != NULL) {
|
||||
if (mxIsDouble(M_pconn)) {
|
||||
pd = mxGetPr(M_pconn);
|
||||
|
||||
for (e = 0; e < ne; e++) { (*pconn)[e] = pd[e] - 1; }
|
||||
} else {
|
||||
pi = mxGetData(M_pconn);
|
||||
|
||||
for (e = 0; e < ne; e++) { (*pconn)[e] = pi[e] - 1; };
|
||||
}
|
||||
|
||||
ret = ne - 1;
|
||||
} else {
|
||||
ret = -1;
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* [S, r, F1, F2, L] = mex_schur_comp(BI, connPos, conns)
|
||||
*/
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
mexFunction(int nlhs, mxArray *plhs[],
|
||||
int nrhs, const mxArray *prhs[])
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ok, nc, nconn_tot, max_nconn, sum_nconn2;
|
||||
int p2, c, i, nconn, *pconn;
|
||||
double *Binv, *ptr1, *ptr2, *src, *gpress, *BIV, *P;
|
||||
struct hybsys *sys;
|
||||
|
||||
ok = verify_args(nlhs, nrhs, prhs);
|
||||
|
||||
if (ok) {
|
||||
nc = get_pconn(prhs[1], &pconn);
|
||||
|
||||
nconn_tot = pconn[nc];
|
||||
max_nconn = count_cellconn(nc, pconn);
|
||||
|
||||
allocate_aux_arrays(nc, nconn_tot, &src, &gpress, &BIV, &P);
|
||||
|
||||
sum_nconn2 = mxGetNumberOfElements(prhs[0]);
|
||||
plhs[0] = mxCreateDoubleMatrix(sum_nconn2, 1, mxREAL);
|
||||
plhs[1] = mxCreateDoubleMatrix(nconn_tot, 1, mxREAL);
|
||||
plhs[2] = mxCreateDoubleMatrix(nconn_tot, 1, mxREAL);
|
||||
plhs[3] = mxCreateDoubleMatrix(nconn_tot, 1, mxREAL);
|
||||
plhs[4] = mxCreateDoubleMatrix(nc, 1, mxREAL);
|
||||
|
||||
sys = hybsys_allocate_unsymm(max_nconn, nc, nconn_tot);
|
||||
hybsys_init(max_nconn, sys);
|
||||
|
||||
for (i = 0; i < nc; i++) { src[i] = 0.0; /* No sources */
|
||||
P[i] = 0.0; } /* No accumulation */
|
||||
for (i = 0; i < nconn_tot; i++) { gpress[i] = 0.0; /* No gravity */
|
||||
BIV[i] = 0.0; } /* No V^2 */
|
||||
|
||||
Binv = mxGetPr(prhs[0]);
|
||||
|
||||
hybsys_schur_comp_unsymm(nc, pconn, Binv, BIV, P, sys);
|
||||
|
||||
ptr1 = mxGetPr(plhs[0]);
|
||||
ptr2 = mxGetPr(plhs[1]);
|
||||
p2 = 0;
|
||||
for (c = 0; c < nc; c++) {
|
||||
nconn = pconn[c + 1] - pconn[c];
|
||||
|
||||
hybsys_cellcontrib_unsymm(c, nconn, pconn[c], p2,
|
||||
gpress, src, Binv, sys);
|
||||
|
||||
memcpy(ptr1 + p2 , sys->S, nconn * nconn * sizeof *ptr1);
|
||||
memcpy(ptr2 + pconn[c], sys->r, nconn * sizeof *ptr2);
|
||||
|
||||
p2 += nconn * nconn;
|
||||
}
|
||||
|
||||
ptr1 = mxGetPr(plhs[2]);
|
||||
memcpy(ptr1, sys->F1, nconn_tot * sizeof *ptr1);
|
||||
|
||||
ptr1 = mxGetPr(plhs[3]);
|
||||
memcpy(ptr1, sys->F2, nconn_tot * sizeof *ptr1);
|
||||
|
||||
ptr1 = mxGetPr(plhs[4]);
|
||||
memcpy(ptr1, sys->L , nc * sizeof *ptr1);
|
||||
|
||||
hybsys_free(sys);
|
||||
deallocate_aux_arrays(pconn, src, gpress, BIV, P);
|
||||
}
|
||||
}
|
||||
64
mex_schur_comp.m
Normal file
64
mex_schur_comp.m
Normal file
@ -0,0 +1,64 @@
|
||||
function varargout = mex_schur_comp_symm(varargin)
|
||||
%Compute hybrid system component matrices using compiled C code.
|
||||
%
|
||||
% SYNOPSIS:
|
||||
% [S, r, F1, F2, L] = mex_schur_comp_symm(BI, connPos, conns)
|
||||
%
|
||||
% PARAMETERS:
|
||||
% BI - Inner product values.
|
||||
%
|
||||
% connPos - Indirection map of size [G.cells.num,1] into 'conns' table
|
||||
% (i.e., the connections or DOFs). Specifically, the DOFs
|
||||
% connected to cell 'i' are found in the submatrix
|
||||
%
|
||||
% conns(connPos(i) : connPos(i + 1) - 1)
|
||||
%
|
||||
% conns - A (connPos(end)-1)-by-1 array of cell connections
|
||||
% (local-to-global DOF mapping in FEM parlance).
|
||||
%
|
||||
% NOTE:
|
||||
% The (connPos,conns) array pair is expected to be the output of function
|
||||
% 'mex_ip_simple'.
|
||||
%
|
||||
% RETURNS:
|
||||
% S - A SUM(nconn .^ 2)-by-1 array of unassembled system matrix values,
|
||||
% ordered by cells.
|
||||
%
|
||||
% r - A SUM(nconn)-by-1 array of unassemble system rhs values, ordered by
|
||||
% cells.
|
||||
%
|
||||
% F1 - A SUM(nconn)-by-1 array of C'*inv(B) values, ordered by cells.
|
||||
%
|
||||
% F2 - A SUM(nconn)-by-1 array of C'*inv(B) values, ordered by cells.
|
||||
%
|
||||
% L - A G.cells.num-by-1 array of C'*inv(B)*C values, ordered by cells.
|
||||
%
|
||||
% SEE ALSO:
|
||||
% mex_ip_simple, mex_compute_press_flux.
|
||||
|
||||
%{
|
||||
#COPYRIGHT#
|
||||
%}
|
||||
|
||||
% $Date$
|
||||
% $Revision$
|
||||
|
||||
buildmex CFLAGS="\$CFLAGS -Wall -Wextra -ansi -pedantic ...
|
||||
-Wformat-nonliteral -Wcast-align -Wpointer-arith ...
|
||||
-Wbad-function-cast -Wmissing-prototypes -Wstrict-prototypes ...
|
||||
-Wmissing-declarations -Winline -Wundef -Wnested-externs ...
|
||||
-Wcast-qual -Wshadow -Wconversion -Wwrite-strings ...
|
||||
-Wno-conversion -Wchar-subscripts -Wredundant-decls" ...
|
||||
...
|
||||
-O -largeArrayDims -DCOMPILING_FOR_MATLAB=1 ...
|
||||
-DASSEMBLE_AND_SOLVE_UMFPACK=0 ...
|
||||
-I/work/include/SuiteSparse/ ...
|
||||
-I/usr/include/suitesparse/ ...
|
||||
...
|
||||
mex_schur_comp.c hybsys.c call_umfpack.c ...
|
||||
...
|
||||
-lmwumfpack -lmwamd -lmwlapack -lmwblas
|
||||
|
||||
% Call MEX'ed edition.
|
||||
[varargout{1:nargout}] = mex_schur_comp_symm(varargin{:});
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user