| 1 |
#include <math.h> |
|---|
| 2 |
#include <stdio.h> |
|---|
| 3 |
#include "mex.h" |
|---|
| 4 |
#include "opnml_mex5_allocs.c" |
|---|
| 5 |
|
|---|
| 6 |
/************************************************************ |
|---|
| 7 |
|
|---|
| 8 |
#### ## ##### ###### # # ## # # |
|---|
| 9 |
# # # # # # # # # # # # |
|---|
| 10 |
# # # # ##### # # # # # |
|---|
| 11 |
# ### ###### # # # ## # ###### # |
|---|
| 12 |
# # # # # # ## ## # # # |
|---|
| 13 |
#### # # # ###### # # # # # |
|---|
| 14 |
|
|---|
| 15 |
************************************************************/ |
|---|
| 16 |
|
|---|
| 17 |
/* ---- AA,BB,TT are defined to perform the following array |
|---|
| 18 |
element extractions AA(i,j,m) AA[i+m*j] ---------------- */ |
|---|
| 19 |
#define AA(i,j,m) A[i+m*j] |
|---|
| 20 |
#define BB(i,j,m) B[i+m*j] |
|---|
| 21 |
#define TT(i,j,m) T[i+m*j] |
|---|
| 22 |
|
|---|
| 23 |
void mexFunction(int nlhs, |
|---|
| 24 |
mxArray *plhs[], |
|---|
| 25 |
int nrhs, |
|---|
| 26 |
const mxArray *prhs[]) |
|---|
| 27 |
{ |
|---|
| 28 |
|
|---|
| 29 |
/* ---- findelemex will be called as : |
|---|
| 30 |
j_el=findelemex(xp,yp,AR,A,B,T); ---------------------------- */ |
|---|
| 31 |
/* ---- xp,yp are NOT nodal coordinates; they are the points we are |
|---|
| 32 |
finding elements for. Nodal coordinates have already been |
|---|
| 33 |
accounted for in A,B,T ----- */ |
|---|
| 34 |
|
|---|
| 35 |
int ip,j,np,nl,nh,ne; |
|---|
| 36 |
double *xp, *yp; |
|---|
| 37 |
double *AR,*A,*B,*T; |
|---|
| 38 |
double *fnd; |
|---|
| 39 |
double NaN=mxGetNaN(); |
|---|
| 40 |
double fac,S1,S2,S3,ONE,ZERO; |
|---|
| 41 |
|
|---|
| 42 |
/* ---- check I/O arguments ----------------------------------------- */ |
|---|
| 43 |
if (nrhs != 6) |
|---|
| 44 |
mexErrMsgTxt("findelemex requires 6 input arguments."); |
|---|
| 45 |
else if (nlhs != 1) |
|---|
| 46 |
mexErrMsgTxt("findelemex requires 1 output arguments."); |
|---|
| 47 |
|
|---|
| 48 |
/* ---- dereference input arrays ------------------------------------ */ |
|---|
| 49 |
xp =mxGetPr(prhs[0]); |
|---|
| 50 |
yp =mxGetPr(prhs[1]); |
|---|
| 51 |
AR =mxGetPr(prhs[2]); |
|---|
| 52 |
A =mxGetPr(prhs[3]); |
|---|
| 53 |
B =mxGetPr(prhs[4]); |
|---|
| 54 |
T =mxGetPr(prhs[5]); |
|---|
| 55 |
|
|---|
| 56 |
np=mxGetM(prhs[0]); |
|---|
| 57 |
ne=mxGetM(prhs[2]); |
|---|
| 58 |
|
|---|
| 59 |
/* ---- allocate space for list containing element numbers following NRC |
|---|
| 60 |
allocation style |
|---|
| 61 |
double *mxDvector(int nl,int nh) |
|---|
| 62 |
fnd= (double *) mxDvector(0,np); ---------------------------- */ |
|---|
| 63 |
fnd= (double *) mxDvector(0,np); |
|---|
| 64 |
for (ip=0;ip<np;ip++)fnd[ip]=-1.; |
|---|
| 65 |
ONE=1.+1e-8; |
|---|
| 66 |
ZERO=0.-1e-8; |
|---|
| 67 |
for (j=0;j<ne;j++){ |
|---|
| 68 |
for (ip=0;ip<np;ip++){ |
|---|
| 69 |
if(fnd[ip]<(double)0){ |
|---|
| 70 |
fac=.5/AR[j]; |
|---|
| 71 |
S1=(TT(j,0,ne)+BB(j,0,ne)*xp[ip]+AA(j,0,ne)*yp[ip])*fac; |
|---|
| 72 |
if (S1>ONE|S1<ZERO)goto l20; |
|---|
| 73 |
S2=(TT(j,1,ne)+BB(j,1,ne)*xp[ip]+AA(j,1,ne)*yp[ip])*fac; |
|---|
| 74 |
if (S2>ONE|S2<ZERO)goto l20; |
|---|
| 75 |
S3=(TT(j,2,ne)+BB(j,2,ne)*xp[ip]+AA(j,2,ne)*yp[ip])*fac; |
|---|
| 76 |
if (S3>ONE|S3<ZERO)goto l20; |
|---|
| 77 |
fnd[ip]=(double)(j+1); |
|---|
| 78 |
} |
|---|
| 79 |
l20: continue; |
|---|
| 80 |
} |
|---|
| 81 |
} |
|---|
| 82 |
for (ip=0;ip<np;ip++) if(fnd[ip]<(double)0)fnd[ip]=NaN; |
|---|
| 83 |
|
|---|
| 84 |
/* ---- Set elements of return matrix, pointed to by plhs[0] -------- */ |
|---|
| 85 |
plhs[0]=mxCreateDoubleMatrix(np,1,mxREAL); |
|---|
| 86 |
mxSetPr(plhs[0],fnd); |
|---|
| 87 |
|
|---|
| 88 |
/* ---- No need to free memory allocated with "mxCalloc"; MATLAB |
|---|
| 89 |
does this automatically. The CMEX allocation functions in |
|---|
| 90 |
"opnml_allocs.c" use mxCalloc. ----------------------------------- */ |
|---|
| 91 |
return; |
|---|
| 92 |
} |
|---|
