| 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 |
/* ---- findelemex52 will be called as : |
|---|
| 30 |
j_el=findelemex52(xp,yp,AR,A,B,T,jsearch); ---------------------------- */ |
|---|
| 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. jsearch os a list of elements to |
|---|
| 34 |
search in. |
|---|
| 35 |
----- */ |
|---|
| 36 |
|
|---|
| 37 |
int ip,j,jj,np,nl,nh,ne,njsearch; |
|---|
| 38 |
double *xp, *yp,*djsearch; |
|---|
| 39 |
double *AR,*A,*B,*T; |
|---|
| 40 |
double *fnd; |
|---|
| 41 |
double NaN=mxGetNaN(); |
|---|
| 42 |
double fac,S1,S2,S3,ONE,ZERO; |
|---|
| 43 |
double tol=1.e-6; |
|---|
| 44 |
|
|---|
| 45 |
|
|---|
| 46 |
|
|---|
| 47 |
/* ---- check I/O arguments ----------------------------------------- */ |
|---|
| 48 |
if (nrhs != 7) |
|---|
| 49 |
mexErrMsgTxt("findelemex52 requires 7 input arguments."); |
|---|
| 50 |
else if (nlhs != 1) |
|---|
| 51 |
mexErrMsgTxt("findelemex52 requires 1 output arguments."); |
|---|
| 52 |
|
|---|
| 53 |
/* ---- dereference input arrays ------------------------------------ */ |
|---|
| 54 |
xp =mxGetPr(prhs[0]); |
|---|
| 55 |
yp =mxGetPr(prhs[1]); |
|---|
| 56 |
AR =mxGetPr(prhs[2]); |
|---|
| 57 |
A =mxGetPr(prhs[3]); |
|---|
| 58 |
B =mxGetPr(prhs[4]); |
|---|
| 59 |
T =mxGetPr(prhs[5]); |
|---|
| 60 |
djsearch=mxGetPr(prhs[6]); |
|---|
| 61 |
|
|---|
| 62 |
np=mxGetM(prhs[0]); |
|---|
| 63 |
ne=mxGetM(prhs[2]); |
|---|
| 64 |
njsearch=mxGetM(prhs[6]); |
|---|
| 65 |
|
|---|
| 66 |
/* ---- allocate space for list containing element numbers following NRC |
|---|
| 67 |
allocation style |
|---|
| 68 |
double *mxDvector(int nl,int nh) |
|---|
| 69 |
fnd= (double *) mxDvector(0,np); ---------------------------- */ |
|---|
| 70 |
|
|---|
| 71 |
fnd= (double *) mxDvector(0,np); |
|---|
| 72 |
for (ip=0;ip<np;ip++)fnd[ip]=-1.; |
|---|
| 73 |
ONE=1.+tol; |
|---|
| 74 |
ZERO=0.-tol; |
|---|
| 75 |
for (jj=0;jj<njsearch;jj++){ |
|---|
| 76 |
j=(int)djsearch[jj]; |
|---|
| 77 |
for (ip=0;ip<np;ip++){ |
|---|
| 78 |
if(fnd[ip]<(double)0){ |
|---|
| 79 |
fac=.5/AR[j]; |
|---|
| 80 |
S1=(TT(j,0,ne)+BB(j,0,ne)*xp[ip]+AA(j,0,ne)*yp[ip])*fac; |
|---|
| 81 |
if (S1>ONE|S1<ZERO)goto l20; |
|---|
| 82 |
S2=(TT(j,1,ne)+BB(j,1,ne)*xp[ip]+AA(j,1,ne)*yp[ip])*fac; |
|---|
| 83 |
if (S2>ONE|S2<ZERO)goto l20; |
|---|
| 84 |
S3=(TT(j,2,ne)+BB(j,2,ne)*xp[ip]+AA(j,2,ne)*yp[ip])*fac; |
|---|
| 85 |
if (S3>ONE|S3<ZERO)goto l20; |
|---|
| 86 |
fnd[ip]=(double)(j+1); |
|---|
| 87 |
} |
|---|
| 88 |
l20: continue; |
|---|
| 89 |
} |
|---|
| 90 |
} |
|---|
| 91 |
for (ip=0;ip<np;ip++) if(fnd[ip]<(double)0)fnd[ip]=NaN; |
|---|
| 92 |
|
|---|
| 93 |
/* ---- Set elements of return matrix, pointed to by plhs[0] -------- */ |
|---|
| 94 |
plhs[0]=mxCreateDoubleMatrix(np,1,mxREAL); |
|---|
| 95 |
mxSetPr(plhs[0],fnd); |
|---|
| 96 |
|
|---|
| 97 |
/* ---- No need to free memory allocated with "mxCalloc"; MATLAB |
|---|
| 98 |
does this automatically. The CMEX allocation functions in |
|---|
| 99 |
"opnml_allocs.c" use mxCalloc. ----------------------------------- */ |
|---|
| 100 |
return; |
|---|
| 101 |
} |
|---|
