| 1 |
#include <math.h> |
|---|
| 2 |
#include <stdio.h> |
|---|
| 3 |
#include "mex.h" |
|---|
| 4 |
#include "matrix.h" |
|---|
| 5 |
#include "opnml_mex5_allocs.c" |
|---|
| 6 |
/* Array Access Definitions */ |
|---|
| 7 |
#define MN m*n |
|---|
| 8 |
#define MNP m*n*p |
|---|
| 9 |
#define ELE(i,j,m) ele[i+m*j] |
|---|
| 10 |
#define PHI(i,j,m) phi[i+m*j] |
|---|
| 11 |
#define JXY(i,j,m) jxy[i+m*j] |
|---|
| 12 |
#define SFEM(i,j,m) sfem[i+m*j] /* SFEM(i,j,nn) */ |
|---|
| 13 |
#define FDZgrid(i,j,k,m,n) fdzgrid[i+m*j+m*n*k] /* FDZgrid(i,j,k,m,n) */ |
|---|
| 14 |
#define N1(i,j,k,m,n) n1[i+m*j+m*n*k] /* N1(i,j,k,m,n) */ |
|---|
| 15 |
#define N2(i,j,k,m,n) n2[i+m*j+m*n*k] /* N2(i,j,k,m,n) */ |
|---|
| 16 |
#define B1(i,j,k,m,n) b1[i+m*j+m*n*k] /* B1(i,j,k,m,n) */ |
|---|
| 17 |
#define B2(i,j,k,m,n) b2[i+m*j+m*n*k] /* B2(i,j,k,m,n) */ |
|---|
| 18 |
#define SFD(i,j,k,m,n) sfd[i+m*j+m*n*k] /* SFD(i,j,k,m,n) */ |
|---|
| 19 |
|
|---|
| 20 |
/************************************************************ |
|---|
| 21 |
|
|---|
| 22 |
#### ## ##### ###### # # ## # # |
|---|
| 23 |
# # # # # # # # # # # # |
|---|
| 24 |
# # # # ##### # # # # # |
|---|
| 25 |
# ### ###### # # # ## # ###### # |
|---|
| 26 |
# # # # # # ## ## # # # |
|---|
| 27 |
#### # # # ###### # # # # # |
|---|
| 28 |
|
|---|
| 29 |
************************************************************/ |
|---|
| 30 |
|
|---|
| 31 |
|
|---|
| 32 |
void mexFunction(int nlhs, |
|---|
| 33 |
mxArray *plhs[], |
|---|
| 34 |
int nrhs, |
|---|
| 35 |
const mxArray *prhs[]) |
|---|
| 36 |
{ |
|---|
| 37 |
/* 0 1 2 3 4 5 6 7 8 9 |
|---|
| 38 |
/* S=map_scalar_mex5(phi,N1,N2,B1,B2,FDZgrid,zprof,scalar,e,jXY); */ |
|---|
| 39 |
|
|---|
| 40 |
int m,n,p,i,j,k,ne,nn,nnv, *ele, *n1,*n2,*jxy; |
|---|
| 41 |
int node1,node2,node3,elem,lup,llow; |
|---|
| 42 |
double s_upper,s_lower; |
|---|
| 43 |
const int *dims; |
|---|
| 44 |
double *phi,*b1,*b2,*temp,*temp1,*temp2,*temp3,*dele,*djxy,*dN1,*dN2; |
|---|
| 45 |
double *sfem, *fdzgrid,*zprof,*sfd; |
|---|
| 46 |
double phi1,phi2,phi3; |
|---|
| 47 |
mxArray *S; |
|---|
| 48 |
double NaN=mxGetNaN(); |
|---|
| 49 |
unsigned char *start_of_array; |
|---|
| 50 |
size_t bytes_to_copy; |
|---|
| 51 |
int number_of_dims; |
|---|
| 52 |
|
|---|
| 53 |
/* get m,n,p from RHS arg #1 */ |
|---|
| 54 |
dims = mxGetDimensions(prhs[1]); |
|---|
| 55 |
m=dims[0];n=dims[1];p=dims[2]; |
|---|
| 56 |
|
|---|
| 57 |
/* phi is (m*n)X3 */ |
|---|
| 58 |
temp=mxGetPr(prhs[0]); |
|---|
| 59 |
phi=(double *)mxDvector(0,3*m*n); |
|---|
| 60 |
for (i=0;i<3*MN;i++) |
|---|
| 61 |
phi[i]=temp[i]; |
|---|
| 62 |
|
|---|
| 63 |
/* The scalar to be mapped is in prhs[10]. get nn,nnv,ne */ |
|---|
| 64 |
nn=mxGetM(prhs[7]); |
|---|
| 65 |
nnv=mxGetN(prhs[7]); |
|---|
| 66 |
ne=mxGetM(prhs[8]); |
|---|
| 67 |
|
|---|
| 68 |
/* Get the FEM-based scalar field */ |
|---|
| 69 |
sfem=mxGetPr(prhs[7]); |
|---|
| 70 |
|
|---|
| 71 |
/* ---- allocate space for int representation of dele & |
|---|
| 72 |
convert double element representation to int & |
|---|
| 73 |
shift node numbers toward 0 by 1 for proper indexing -------- */ |
|---|
| 74 |
dele=mxGetPr(prhs[8]); |
|---|
| 75 |
ele=(int *)mxIvector(0,3*ne); |
|---|
| 76 |
for (i=0;i<3*ne;i++){ |
|---|
| 77 |
ele[i]=(int)dele[i]; |
|---|
| 78 |
ele[i]=ele[i]-1; |
|---|
| 79 |
} |
|---|
| 80 |
|
|---|
| 81 |
djxy=mxGetPr(prhs[9]); |
|---|
| 82 |
jxy=(int *)mxIvector(0,MN*3); |
|---|
| 83 |
for (i=0;i<3*MN;i++){ |
|---|
| 84 |
jxy[i]=(int)djxy[i]; |
|---|
| 85 |
jxy[i]=jxy[i]-1; |
|---|
| 86 |
} |
|---|
| 87 |
|
|---|
| 88 |
/* N1,etc is mXnXp */ |
|---|
| 89 |
dN1=mxGetPr(prhs[1]); |
|---|
| 90 |
dN2=mxGetPr(prhs[2]); |
|---|
| 91 |
b1=mxGetPr(prhs[3]); /* B1 */ |
|---|
| 92 |
b2=mxGetPr(prhs[4]); /* B2 */ |
|---|
| 93 |
fdzgrid=mxGetPr(prhs[5]); /* FDZgrid */ |
|---|
| 94 |
n1=(int *)mxIvector(0,MNP); |
|---|
| 95 |
n2=(int *)mxIvector(0,MNP); |
|---|
| 96 |
for (i=0;i<MNP;i++){ |
|---|
| 97 |
n1[i]=(int)dN1[i]-1; |
|---|
| 98 |
n2[i]=(int)dN2[i]-1; |
|---|
| 99 |
} |
|---|
| 100 |
|
|---|
| 101 |
|
|---|
| 102 |
/* Get the FD depth profile */ |
|---|
| 103 |
zprof=mxGetPr(prhs[6]); |
|---|
| 104 |
|
|---|
| 105 |
S=mxCreateNumericArray(3,dims,mxDOUBLE_CLASS,mxREAL); |
|---|
| 106 |
sfd=(double *)mxDvector(0,MNP); |
|---|
| 107 |
|
|---|
| 108 |
/* time to map the scalar fem-based field to the FD array */ |
|---|
| 109 |
for (i=0;i<m;i++){ |
|---|
| 110 |
for (j=0;j<n;j++){ |
|---|
| 111 |
phi1=PHI(i+j*n,0,MN); |
|---|
| 112 |
phi2=PHI(i+j*n,1,MN); |
|---|
| 113 |
phi3=PHI(i+j*n,2,MN); |
|---|
| 114 |
elem=JXY(i,j,m); |
|---|
| 115 |
node1=ELE(elem,0,ne); |
|---|
| 116 |
node2=ELE(elem,1,ne); |
|---|
| 117 |
node3=ELE(elem,2,ne); |
|---|
| 118 |
if (mxIsNaN(FDZgrid(i,j,0,m,n))){ |
|---|
| 119 |
for (k=0;k<p;k++) |
|---|
| 120 |
SFD(i,j,k,m,n)=NaN; /* This is an "outside the domain" check!!*/ |
|---|
| 121 |
} |
|---|
| 122 |
else { |
|---|
| 123 |
for (k=0;k<p;k++){ |
|---|
| 124 |
if(zprof[k]<FDZgrid(i,j,0,m,n)) |
|---|
| 125 |
SFD(i,j,k,m,n)=NaN; |
|---|
| 126 |
else{ |
|---|
| 127 |
llow=N1(i,j,k,m,n); |
|---|
| 128 |
if (llow==nnv-1) |
|---|
| 129 |
SFD(i,j,k,m,n)=SFEM(node1,llow,nn)*phi1+ |
|---|
| 130 |
SFEM(node2,llow,nn)*phi2+ |
|---|
| 131 |
SFEM(node3,llow,nn)*phi3; |
|---|
| 132 |
|
|---|
| 133 |
else{ |
|---|
| 134 |
lup=N2(i,j,k,m,n); |
|---|
| 135 |
s_upper=SFEM(node1,lup,nn)*phi1+ |
|---|
| 136 |
SFEM(node2,lup,nn)*phi2+ |
|---|
| 137 |
SFEM(node3,lup,nn)*phi3; |
|---|
| 138 |
s_lower=SFEM(node1,llow,nn)*phi1+ |
|---|
| 139 |
SFEM(node2,llow,nn)*phi2+ |
|---|
| 140 |
SFEM(node3,llow,nn)*phi3; |
|---|
| 141 |
SFD(i,j,k,m,n)=s_upper*B1(i,j,k,m,n) + s_lower*B2(i,j,k,m,n); |
|---|
| 142 |
} |
|---|
| 143 |
} |
|---|
| 144 |
} |
|---|
| 145 |
} |
|---|
| 146 |
} |
|---|
| 147 |
} |
|---|
| 148 |
|
|---|
| 149 |
|
|---|
| 150 |
/* Set pointer to lhs */ |
|---|
| 151 |
start_of_array=(unsigned char *)mxGetPr(S); |
|---|
| 152 |
bytes_to_copy=MNP*mxGetElementSize(S); |
|---|
| 153 |
memcpy(start_of_array,sfd,bytes_to_copy); |
|---|
| 154 |
plhs[0]=S; |
|---|
| 155 |
return; |
|---|
| 156 |
} |
|---|
