when inserting mexPrintf commands I get different outputs

12 Jul. 2011
4 Antworten
4 Ansichten (30 Tage)

0 Stimmen

[EDIT: 20110711 21:16 CDT - merge Answer, reformat - WDR]
Hi,
I'm writing a mex file and I have observed something very weird. When I'm inserting code like the following:
mexPrintf("yi_k %e\n",yi_k[2]);
in various places to print the values of the variables, sometimes I get different results although this has practically no effect on the code.
Has anyone any idea what to look for?
Thank you in advance
[EDIT: bring code up from Answer]
Here is the code which is very long
/* Version 3 forces velocity of top layer to be zero when
*recharge is zero
*
*/
#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include "mex.h"
/* Input Arguments */
#define A_IN prhs[0]
#define B_IN prhs[1]
#define ND_IN prhs[2]
#define VX_IN prhs[3]
#define VY_IN prhs[4]
#define VZ_IN prhs[5]
#define SX_IN prhs[6]
#define SY_IN prhs[7]
#define SZ_IN prhs[8]
#define ID_IN prhs[9]
#define STEP_IN prhs[10]
#define CS_IN prhs[11]
#define CS1_IN prhs[12]
#define VTOP_IN prhs[13]
/* Output Arguments */
#define XYZV_OUT plhs[0]
#define ISIDE_OUT plhs[1]
/*****************************************************
*
* PlaneLineCross
*
******************************************************/
void PlaneLineCross(double T[],double PL[], double L[])
{
/* L[0] L[2] L[4]
* L[1] L[3] L[5]
* PL[0] PL[3] PL[6]
* PL[1] PL[4] PL[7]
* PL[2] PL[5] PL[8]
*/
double a,b,c,d,e,f,g,h,i;
double A,B,C,D,E,F,G,H,I;
double x1,x2,x3;
a=L[0]-L[1]; b=PL[1]-PL[0]; c=PL[2]-PL[0];
d=L[2]-L[3]; e=PL[4]-PL[3]; f=PL[5]-PL[3];
g=L[4]-L[5]; h=PL[7]-PL[6]; i=PL[8]-PL[6];
A=(e*i - f*h)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
B=-(b*i - c*h)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
C=(b*f - c*e)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
D=-(d*i - f*g)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
E=(a*i - c*g)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
F=-(a*f - c*d)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
G=(d*h - e*g)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
H=-(a*h - b*g)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
I=(a*e - b*d)/(a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g);
x1=L[0]-PL[0];x2=L[2]-PL[3];x3=L[4]-PL[6];
/* Matrix multiplication
* A B C x1
* D E F * x2
* G H I x3
*/
T[0]=A*x1+B*x2+C*x3;
T[1]=D*x1+E*x2+F*x3;
T[2]=G*x1+H*x2+I*x3;
}
/*****************************************************
*
* in3dpolygon
*
******************************************************/
bool in3dpolygon(double A[],double ND[],int ii ,double p[])
{
int i,index;
double mnz=10000000000;
double mxz=-10000000000;
double tempnd[18];
double PL[9];
double Tb[3];
double Tt[3];
double L[6];
bool in=false;
//mexPrintf("p %e",p[0]);mexPrintf("p %e",p[1]);mexPrintf("p %e\n",p[2]);
for(i=0;i<6;i++)
{
index=(int)A[i+18*(ii-1)];
//mexPrintf("elindex %d\n",index);
tempnd[i+6*0]=ND[0+3*(index-1)];
tempnd[i+6*1]=ND[1+3*(index-1)];
tempnd[i+6*2]=ND[2+3*(index-1)];
}
for (i=0;i<6;i++)
{
if (tempnd[i+6*2]>mxz)
mxz=tempnd[i+6*2];
if (tempnd[i+6*2]<mnz)
mnz=tempnd[i+6*2];
}
if (p[2]>mxz || p[2]<mnz)
return(in);
L[0]=p[0];L[2]=p[1];L[4]=p[2];
L[1]=p[0];L[3]=p[1];L[5]=mxz+(mxz-mnz);
PL[0]=tempnd[0];PL[3]=tempnd[6];PL[6]=tempnd[12];
PL[1]=tempnd[1];PL[4]=tempnd[7];PL[7]=tempnd[13];
PL[2]=tempnd[2];PL[5]=tempnd[8];PL[8]=tempnd[14];
PlaneLineCross(Tb,PL,L);
//mexPrintf("Tb %e",Tb[0]);mexPrintf("Tb %e",Tb[1]);mexPrintf("Tb %e\n",Tb[2]);
PL[0]=tempnd[3];PL[3]=tempnd[9];PL[6]=tempnd[15];
PL[1]=tempnd[4];PL[4]=tempnd[10];PL[7]=tempnd[16];
PL[2]=tempnd[5];PL[5]=tempnd[11];PL[8]=tempnd[17];
PlaneLineCross(Tt, PL,L);
//mexPrintf("Tt %e",Tt[0]);mexPrintf("Tt %e",Tt[1]);mexPrintf("Tt %e\n",Tt[2]);
if ((Tt[0]>=0 && Tt[0]<=1) && Tb[0]<0 && (Tt[1]>=0 && Tt[1]<=1) && (Tt[2]>=0 && Tt[2]<=1) && (Tt[1]+Tt[2]<=1))
in=true;
//mexPrintf("p %e",p[0]);mexPrintf("p %e",p[1]);mexPrintf("p %e\n",p[2]);
return(in);
}
/*****************************************************
*
* listtemp
*
******************************************************/
void listtemp(double B[], double lst[],int inds[])
{
int cnt=1;
int istr;
int iend;
istr=inds[0];
iend=inds[1];
bool dontuse=false;
int ii, i, k;
for(ii=istr-1;ii<iend;ii++)
{
for (i=0;i<5;i++)
{
if (B[i+5*((int)lst[ii]-1)]!=0)
{
for (k=0;k<iend+cnt-1;k++)
{
if (lst[k]==B[i+5*((int)lst[ii]-1)])
{
dontuse=true;
break;
}
}
if (dontuse==true)
dontuse=false;
else
{
lst[iend+cnt-1]=B[i+5*((int)lst[ii]-1)];
cnt++;
}
}
}
}
istr=iend+1;
iend=cnt+iend-1;
inds[0]=istr;
inds[1]=iend;
//mexPrintf("istr_1 %d ",istr);mexPrintf("iend_1 %d\n",iend);
}
/*****************************************************
*
* ListNeighborElem
*
******************************************************/
void ListNeighborElem(double B[],double lst[], int nGen)
{
int igen=1;
int inds[2];
inds[0]=1;inds[1]=1;
//mexPrintf("lst(1) %f\n",lst[0]);
//mexPrintf("lst(2) %f\n",lst[1]);
while (igen<nGen)
{
//mexPrintf("igen %d\n",igen);
//mexPrintf("istr %d ",inds[0]);mexPrintf("iend %d\n",inds[1]);
listtemp(B,lst,inds);
igen++;
}
}
/******************************************************
*
* CalcVel
*
******************************************************/
void CalcVel(double Vel[],double X[],double V[],double x[], double vtop)
{
//mexPrintf("x %e",x[0]);mexPrintf("x %e",x[1]);mexPrintf("x %e\n",x[2]);
/*
x=X(0:17)
y=Y(18:35)
z=Z(36:53)
vx=V(0:17)
vy=V(18:35)
vz=V(36:53)
nodes 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
--------------------------------------------------------------
x or vx 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
y or vy 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
z or vz 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
*/
double CT[9];
double z[3];
double DD, zt, zb, lamd;
// calculate the triangular coordinates from cartesian
//D=1/((x2*y3-x3*y2)+x1*(y2-y3)+y1*(x3-x2)) D is twice the are of the tringle projected in 2D
DD=(X[1]*X[20]-X[2]*X[19])+X[0]*(X[19]-X[20])+X[18]*(X[2]-X[1]);
DD=1/DD;
//CT=[x2*y3-x3*y2 y2-y3 x3-x2;
// x3*y1-x1*y3 y3-y1 x1-x3;
// x1*y2-x2*y1 y1-y2 x2-x1]
CT[0]=X[1]*X[20]-X[2]*X[19]; CT[3]=X[19]-X[20]; CT[6]=X[2]-X[1];
CT[1]=X[2]*X[18]-X[0]*X[20]; CT[4]=X[20]-X[18]; CT[7]=X[0]-X[2];
CT[2]=X[0]*X[19]-X[1]*X[18]; CT[5]=X[18]-X[19]; CT[8]=X[1]-X[0];
// z contains the triangular coordinates
z[0]=DD*(CT[0]+CT[3]*x[0]+CT[6]*x[1]);
z[1]=DD*(CT[1]+CT[4]*x[0]+CT[7]*x[1]);
z[2]=DD*(CT[2]+CT[5]*x[0]+CT[8]*x[1]);
zt=z[0]*X[39]+z[1]*X[40]+z[2]*X[41];
zb=z[0]*X[36]+z[1]*X[37]+z[2]*X[38];
lamd=2*(x[2]-zb)/(zt-zb)-1;
// interpolate vx
Vel[0]=((z[0]*(2*z[0]-1)*(pow(lamd,2)-lamd))/2)*V[0]+
((z[1]*(2*z[1]-1)*(pow(lamd,2)-lamd))/2)*V[1]+
((z[2]*(2*z[2]-1)*(pow(lamd,2)-lamd))/2)*V[2]+
((z[0]*(2*z[0]-1)*(pow(lamd,2)+lamd))/2)*V[3]+
((z[1]*(2*z[1]-1)*(pow(lamd,2)+lamd))/2)*V[4]+
((z[2]*(2*z[2]-1)*(pow(lamd,2)+lamd))/2)*V[5]+
2*z[0]*z[1]*(pow(lamd,2)-lamd)*V[6]+
2*z[2]*z[0]*(pow(lamd,2)-lamd)*V[7]+
2*z[1]*z[2]*(pow(lamd,2)-lamd)*V[8]+
z[0]*(2*z[0]-1)*(1-pow(lamd,2))*V[9]+
4*z[0]*z[1]*(1-pow(lamd,2))*V[10]+
z[1]*(2*z[1]-1)*(1-pow(lamd,2))*V[11]+
4*z[2]*z[0]*(1-pow(lamd,2))*V[12]+
4*z[1]*z[2]*(1-pow(lamd,2))*V[13]+
z[3]*(2*z[3]-1)*(1-pow(lamd,2))*V[14]+
2*z[0]*z[1]*(pow(lamd,2)+lamd)*V[15]+
2*z[2]*z[0]*(pow(lamd,2)+lamd)*V[16]+
2*z[1]*z[2]*(pow(lamd,2)+lamd)*V[17];
// interpolate vy
Vel[1]=((z[0]*(2*z[0]-1)*(pow(lamd,2)-lamd))/2)*V[18]+
((z[1]*(2*z[1]-1)*(pow(lamd,2)-lamd))/2)*V[19]+
((z[2]*(2*z[2]-1)*(pow(lamd,2)-lamd))/2)*V[20]+
((z[0]*(2*z[0]-1)*(pow(lamd,2)+lamd))/2)*V[21]+
((z[1]*(2*z[1]-1)*(pow(lamd,2)+lamd))/2)*V[22]+
((z[2]*(2*z[2]-1)*(pow(lamd,2)+lamd))/2)*V[23]+
2*z[0]*z[1]*(pow(lamd,2)-lamd)*V[24]+
2*z[2]*z[0]*(pow(lamd,2)-lamd)*V[25]+
2*z[1]*z[2]*(pow(lamd,2)-lamd)*V[26]+
z[0]*(2*z[0]-1)*(1-pow(lamd,2))*V[27]+
4*z[0]*z[1]*(1-pow(lamd,2))*V[28]+
z[1]*(2*z[1]-1)*(1-pow(lamd,2))*V[29]+
4*z[2]*z[0]*(1-pow(lamd,2))*V[30]+
4*z[1]*z[2]*(1-pow(lamd,2))*V[31]+
z[3]*(2*z[3]-1)*(1-pow(lamd,2))*V[32]+
2*z[0]*z[1]*(pow(lamd,2)+lamd)*V[33]+
2*z[2]*z[0]*(pow(lamd,2)+lamd)*V[34]+
2*z[1]*z[2]*(pow(lamd,2)+lamd)*V[35];
// interpolate vz
if (vtop==1)
{
Vel[2]=((z[0]*(2*z[0]-1)*(pow(lamd,2)-lamd))/2)*V[36]+
((z[1]*(2*z[1]-1)*(pow(lamd,2)-lamd))/2)*V[37]+
((z[2]*(2*z[2]-1)*(pow(lamd,2)-lamd))/2)*V[38]+
((z[0]*(2*z[0]-1)*(pow(lamd,2)+lamd))/2)*V[39]+
((z[1]*(2*z[1]-1)*(pow(lamd,2)+lamd))/2)*V[40]+
((z[2]*(2*z[2]-1)*(pow(lamd,2)+lamd))/2)*V[41]+
2*z[0]*z[1]*(pow(lamd,2)-lamd)*V[42]+
2*z[2]*z[0]*(pow(lamd,2)-lamd)*V[43]+
2*z[1]*z[2]*(pow(lamd,2)-lamd)*V[44]+
z[0]*(2*z[0]-1)*(1-pow(lamd,2))*V[45]+
4*z[0]*z[1]*(1-pow(lamd,2))*V[46]+
z[1]*(2*z[1]-1)*(1-pow(lamd,2))*V[47]+
4*z[2]*z[0]*(1-pow(lamd,2))*V[48]+
4*z[1]*z[2]*(1-pow(lamd,2))*V[49]+
z[3]*(2*z[3]-1)*(1-pow(lamd,2))*V[50]+
2*z[0]*z[1]*(pow(lamd,2)+lamd)*V[51]+
2*z[2]*z[0]*(pow(lamd,2)+lamd)*V[52]+
2*z[1]*z[2]*(pow(lamd,2)+lamd)*V[53];
}
else
{
Vel[2]=0;
}
//mexPrintf("Vel %e",Vel[0]);mexPrintf("Vel %e",Vel[1]);mexPrintf("Vel %e\n",Vel[2]);
}
/******************************************************
******************************************************
**
** FindNextPoint
** FindNextPoint(p,vp,A,B,ND,vx,vy,vz,cs,cs1,step,sides);
*****************************************************
*****************************************************/
void FindNextPoint(double p[], double vp[], double A[], double B[],
double ND[], double vx[], double vy[], double vz[],
double cs[], double cs1[], double step[], double sides[], double Vtop[])
{
int i, j, index, cnt;
double thres=0.0001;
double err;
double hhh;
double tempND[54];
double tempV[54];
double yi[3];
double vi[3];
double yi_k[3];
double vi_k[3];
double yi_kk[3];
bool exitdom=false;
bool in,found;
double lst[1000];
int idel;
double idtemp[1];
double vtop;
idtemp[0]=sides[0];
yi[0]=p[0];yi[1]=p[1];yi[2]=p[2];
vi[0]=vp[0];vi[1]=vp[1];vi[2]=vp[2];
//mexPrintf("yi %e",yi[0]);mexPrintf("yi %e",yi[1]);mexPrintf("yi %e\n",yi[2]);
//mexPrintf("vi %e",vi[0]);mexPrintf("vi %e",vi[1]);mexPrintf("vi %e\n",vi[2]);
hhh=step[0]/sqrt(pow(vi[0],2)+pow(vi[1],2)+pow(vi[2],2));
//mexPrintf("h %e\n",hhh);
//take the first step
for(i=0;i<3;i++)
{
yi_k[i]=yi[i]+vi[i]*hhh;
}
//mexPrintf("yi_k %e",yi_k[0]);mexPrintf("yi_k %e",yi_k[1]);mexPrintf("yi_k %e\n",yi_k[2]);
err=10000000;cnt=0;
//Euler predictor corrector scheme
while (err>thres)
{
//mexPrintf("cnt %d\n",cnt);
idel=(int)idtemp[0];
//mexPrintf("idel %d\n",idel);
in=in3dpolygon(A,ND,idel,yi_k);
//mexPrintf("in out %d\n",in);
if (in==false)
{
found=false;
// make list of neibors
//set lst equal to zero
for(i=0;i<1000;i++)
{
lst[i]=0;
}
lst[0]=idtemp[0];
//mexPrintf("lst(1) %f\n",lst[0]);
ListNeighborElem(B,lst,6);
for(i=1;i<1000;i++)
{
//mexPrintf("lst %d %f\n",i, lst[i]);
if (lst[i]==0)
break;
idel=(int)lst[i];
//mexPrintf("idel_11 %d\n",idel);
in=in3dpolygon(A,ND,idel,yi_k);
//mexPrintf("in %d\n",in);
if (in==true)
{
idtemp[0]=lst[i];
//mexPrintf("newElem %d\n",idtemp[0]);
found=true;
break;
}
}
if (found==false)
{
exitdom=true;
break;
}
}
if (exitdom==false)
{
//mexPrintf("idtemp %e\n",idtemp[0]);
for(i=0;i<18;i++)
{
index=(int)A[i+18*(int)(idtemp[0]-1)];
tempND[i+18*0]=ND[0+3*(index-1)];
tempND[i+18*1]=ND[1+3*(index-1)];
tempND[i+18*2]=ND[2+3*(index-1)];
tempV[i+18*0]=vx[0+1*(index-1)];
tempV[i+18*1]=vy[0+1*(index-1)];
tempV[i+18*2]=vz[0+1*(index-1)];
}
vtop=Vtop[(int)(idtemp[0]-1)];
//mexPrintf("yi_k %e",yi_k[0]);mexPrintf("yi_k %e",yi_k[1]);mexPrintf("yi_k %e\n",yi_k[2]);
CalcVel(vi_k,tempND,tempV,yi_k,vtop);
//mexPrintf("vi_k %e",vi_k[0]);mexPrintf("vi %e",vi_k[1]);mexPrintf("vi %e\n",vi_k[2]);
for (j=0;j<3;j++)
{
yi_kk[j]=yi[j]+(0.5*(vi[j]+vi_k[j]))*hhh;//(step[0]/sqrt(pow(0.5*(vi[0]+vi_k[0]),2)+pow(0.5*(vi[1]+vi_k[1]),2)+pow(0.5*(vi[2]+vi_k[2]),2)));
}
//mexPrintf("yi %e",yi[0]);mexPrintf("yi %e",yi[1]);mexPrintf("yi %e\n",yi[2]);
//mexPrintf("vi %e",vi[0]);mexPrintf("vi %e",vi[1]);mexPrintf("vi %e\n",yi[2]);
//mexPrintf("vi_k %e",vi_k[0]);mexPrintf("vi_k %e",vi_k[1]);mexPrintf("vi_k %e\n",vi_k[2]);
//mexPrintf("yi_kk %e",yi_kk[0]);mexPrintf("yi_kk %e",yi_kk[1]);mexPrintf("yi_kk %e\n",yi_kk[2]);
err=sqrt(pow((yi_kk[0]-yi_k[0]),2)+pow((yi_kk[1]-yi_k[1]),2)+pow((yi_kk[2]-yi_k[2]),2));
//mexPrintf("err %e\n",err);
yi_k[0]=yi_kk[0];yi_k[1]=yi_kk[1];yi_k[2]=yi_kk[2];
//if the error is larger than 10^-4 after 20 iteretions
//use the explicit euler with x times smaller than step
if (cnt>20)
{
//mexPrintf("cnt %d\n",cnt);
for (j=0;j<3;j++)
{
yi_k[j]=yi[j]+0.5*(vi[j]+vi_k[j])*(hhh);//(step[0]/2)/sqrt(pow(0.5*(vi[0]+vi_k[0]),2)+pow(0.5*(vi[1]+vi_k[1]),2)+pow(0.5*(vi[2]+vi_k[2]),2));
}
err=0;
}
}
cnt++;
}
//mexPrintf("sides[0] %d\n",sides[0]);
if (exitdom==true)
{
p[0]=yi_k[0];p[1]=yi_k[1];p[2]=yi_k[2];
vp[0]=vi_k[0];vp[1]=vi_k[1];vp[2]=vi_k[2];
sides[2]=sides[0];
sides[0]=0;
}
else
{
//mexPrintf("idel %d\n",idel);
idel=(int)idtemp[0];
in=in3dpolygon(A,ND,idel,yi_k);
// its rather ulikely to be out of polygon now but just in case
//make a check and correct again if needed
if (in==false)
{
found=false;
// make list of neibors
//set lst equal to zero
for(i=0;i<1000;i++)
{
lst[i]=0;
}
lst[0]=idtemp[0];
ListNeighborElem(B,lst,6);
for(i=1;i<1000;i++)
{
if (lst[i]==0)
break;
idel=(int)lst[i];
in=in3dpolygon(A,ND,idel,yi_k);
if (in==true)
{
if (idtemp[0]!=sides[0])
{
sides[2]=sides[0];
sides[0]=lst[i];
}
found=true;
break;
}
}
if (found==false)
{
exitdom=true;
}
}
if (exitdom==false)
{
// mexPrintf("idtemp[0] %f\n",idtemp[0]);
//mexPrintf("sides[0] %f\n",sides[0]);
if (idtemp[0]!=sides[0])
{
sides[2]=sides[0];
sides[0]=idtemp[0];
//mexPrintf("sides[0] %f\n",sides[0]);
}
for(i=0;i<18;i++)
{
index=(int)A[i+18*(int)(sides[0]-1)];
tempND[i+18*0]=ND[0+3*(index-1)];
tempND[i+18*1]=ND[1+3*(index-1)];
tempND[i+18*2]=ND[2+3*(index-1)];
tempV[i+18*0]=vx[0+1*(index-1)];
tempV[i+18*1]=vy[0+1*(index-1)];
tempV[i+18*2]=vz[0+1*(index-1)];
}
vtop=Vtop[(int)(sides[0]-1)];
CalcVel(vi_k,tempND,tempV,yi_k,vtop);
p[0]=yi_k[0];p[1]=yi_k[1];p[2]=yi_k[2];
vp[0]=vi_k[0];vp[1]=vi_k[1];vp[2]=vi_k[2];
//mexPrintf("p %e",p[0]);mexPrintf("p %e",p[1]);mexPrintf("p %e\n",p[2]);
//mexPrintf("vp %e",vp[0]);mexPrintf("vp %e",vp[1]);mexPrintf("vp %e\n",vp[2]);
}
else
{
p[0]=yi_k[0];p[1]=yi_k[1];p[2]=yi_k[2];
vp[0]=vi_k[0];vp[1]=vi_k[1];vp[2]=vi_k[2];
sides[2]=sides[0];
sides[0]=0;
}
}
}
/******************************************************
******************************************************
*
** movebackwards
**
******************************************************
******************************************************/
int
movebackwards(double XYZV[], double sides[], double A[], double B[], double ND[],
double vx[], double vy[], double vz[], double sx[], double sy[], double sz[],
double step[], double cs[], double cs1[], double id[], double Vtop[])
{
int iterations = 0;
bool prtclexit = false;
double vp[3]; // velocity vector
double tempND[54];
double tempV[54];
int i, index, cnt;
double p[3]; //position vevtor
double pprv[3]; // previous position vector
int idside;
double vtop;
p[0]=sx[0];
p[1]=sy[0];
p[2]=sz[0];
cnt=0;
sides[0]=id[0];sides[1]=9;sides[2]=9;sides[3]=-999;
// Calculate initial velocity
for(i=0;i<18;i++)
{
index=(int)A[i+18*(int)(sides[0]-1)];
tempND[i+18*0]=ND[0+3*(index-1)];
tempND[i+18*1]=ND[1+3*(index-1)];
tempND[i+18*2]=ND[2+3*(index-1)];
tempV[i+18*0]=vx[0+1*(index-1)];
tempV[i+18*1]=vy[0+1*(index-1)];
tempV[i+18*2]=vz[0+1*(index-1)];
}
//mexPrintf("X %e",tempND[0]);mexPrintf("X %e",tempND[1]);mexPrintf("X %e\n",tempND[2]);
vtop=Vtop[(int)(sides[0]-1)];
CalcVel(vp,tempND,tempV,p,vtop);
//mexPrintf("vp %e",vp[0]);mexPrintf("vp %e",vp[1]);mexPrintf("vp %e\n",vp[2]);
XYZV[6*iterations + 0] = p[0];
XYZV[6*iterations + 1] = p[1];
XYZV[6*iterations + 2] = p[2];
XYZV[6*iterations + 3] = vp[0];
XYZV[6*iterations + 4] = vp[1];
XYZV[6*iterations + 5] = vp[2];
iterations++;
while(prtclexit==false)
{
FindNextPoint(p,vp,A,B,ND,vx,vy,vz,cs,cs1,step,sides,Vtop);
//mexPrintf("p %e",p[0]);mexPrintf("p %e",p[1]);mexPrintf("p %e\n",p[2]);
//mexPrintf("vp %e",vp[0]);mexPrintf("vp %e",vp[1]);mexPrintf("vp %e\n",vp[2]);
//mexPrintf("sides[0] %f\n",sides[0]);
XYZV[6*iterations + 0] = p[0];
XYZV[6*iterations + 1] = p[1];
XYZV[6*iterations + 2] = p[2];
XYZV[6*iterations + 3] = vp[0];
XYZV[6*iterations + 4] = vp[1];
XYZV[6*iterations + 5] = vp[2];
iterations++;
//mexPrintf("X %e\n",p[0]);
if (sides[0]==0)
{
/*
mexPrintf("iterations %d\n",iterations);
//find from which side the particle exit the element
//take a step first
pprv[0]=XYZV[4*(iterations-2) + 0];
pprv[1]=XYZV[4*(iterations-2) + 1];
pprv[2]=XYZV[4*(iterations-2) + 2];
mexPrintf("pprv %e",pprv[0]);mexPrintf("pprv %e",pprv[1]);mexPrintf("pprv %e\n",pprv[2]);
mexPrintf("p %e",p[0]);mexPrintf("p %e",p[1]);mexPrintf("p %e\n",p[2]);
//for(i=0;i<3;i++)
//{
// pl[i]=p[i]+vp[i]*step[0]/sqrt(pow(vp[0],2)+pow(vp[1],2)+pow(vp[2],2));
//}
idside=PlaneLineCrossF(A,ND,cs,cs1,pprv,p,sides);
sides[1]=idside;
*/
prtclexit=true;
}
if (iterations>19000)
prtclexit=true;
}
return iterations;
}
/******************************************************
* MAIN
*****************************************************/
void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{
double *A, *B, *ND, *vx, *vy, *vz, *sx, *sy, *sz, *id, *step, *cs, *cs1, *TopId, *Vtop;
double *XYZV, *ISIDE;
int iter;
/* get the INPUTS */
A = mxGetPr( A_IN );
B = mxGetPr( B_IN );
ND = mxGetPr( ND_IN );
vx = mxGetPr( VX_IN );
vy = mxGetPr( VY_IN );
vz = mxGetPr( VZ_IN );
sx = mxGetPr( SX_IN );
sy = mxGetPr( SY_IN );
sz = mxGetPr( SZ_IN );
id = mxGetPr( ID_IN );
step = mxGetPr( STEP_IN );
cs = mxGetPr( CS_IN );
cs1 = mxGetPr( CS1_IN );
Vtop = mxGetPr( VTOP_IN );
/* Create matrices for the return arguments */
XYZV_OUT = mxCreateDoubleMatrix( 6, 20000, mxREAL );
ISIDE_OUT = mxCreateDoubleMatrix( 1, 4, mxREAL );
/* Assign pointers to the various parameters */
XYZV = mxGetPr( XYZV_OUT );
ISIDE = mxGetPr( ISIDE_OUT );
iter = movebackwards(XYZV, ISIDE, A, B, ND, vx, vy, vz, sx, sy, sz, step, cs, cs1,id, Vtop);
ISIDE[3]=iter;
}

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (4)

Jan
Jan am 12 Jul. 2011

1 Stimme

It would be helpful, if you find out which mexPrintf causes the difference. At first glance it seems, like e.g.:
mexPrintf("sides[0] %d\n", sides[0]);
could cause troubles, because "sides" is a DOUBLE array, but %d expects an integer. This will confuse the stack. And it does confuse me: Your code is hard to read. I'd be too old to maintain and debug this code.
The massively repeated calculations, e.g. of "(pow(lamd,2)-lamd)", wastes a lot of energy.
Walter Roberson
Walter Roberson am 12 Jul. 2011

0 Stimmen

You indicated that sometimes you get different answers, but you did not describe how the answers differ.
One thing you should keep in mind is that you have programmed in mex file in C, which has different rules than MATLAB does. One of the rules of C is that the C optimizer is allowed to produce code that is only correct if you have not violated any of the many subtle rules of C. In theory, if you have a semantic error anywhere in your code that is certain to be reached, then the C compiler is allowed to generate code that does anything, including creating nasal demons.
There are also sometimes compiler bugs having to do with "register spill". I seem to recall reading a number of years ago that lcc32 had some register spill problems. I also seem to recall that some of the Microsoft C compilers violated IEEE 754 by using chained 80 bit arithmetic operations (a violation of the C rules that the operations must be computed "as if" the values were stored to memory after each operation.) I did not, however, keep track of version numbers for either problem, as I do not build code on MS Windows.
James Tursa
James Tursa am 12 Jul. 2011

0 Stimmen

I haven't examined your code, but I would look for missing prototypes for the functions you use. E.g., it could be that a function returns a char but you didn't prototype it. So the compiler thinks that it returns an int and grabs 32-bits from a register for the return value instead of 8-bits, some of which are garbage. But the mexPrintf function returns an int so when you use it, it effectively clears out all of those other garbage bits, thus changing an 8-bit return value downstream from what it was without the mexPrintf. A nice compiler will warn you about missing prototypes, but compilers are not required to do this so don't necessarily count on getting these warning messages.

1 Kommentar
-1 ältere Kommentare anzeigen -1 ältere Kommentare ausblenden

James Tursa
James Tursa am 12 Jul. 2011
E.g., try putting in a different function rather than mexPrintf (any arbitrary function will do), and experiment with making its return value a char, a short, and an int, and see what happens.

Melden Sie sich an, um zu kommentieren.

Giorgos Kourakos
Giorgos Kourakos am 12 Jul. 2011

0 Stimmen

Thank you all for the responses
What I get from all the answers is that It would be extremely difficult to identify the error in this code. Therefore I wrote individual routines in cpp which perform the most time consuming and repetitive tasks and I kept the main program with the while loops in Matlab. This way it seems to be a bit slower, but I feel I have the control of the program again. In addition it does the calculations I want.
Jan Simon: thanks for your comment about the pow() I've already change that
Giorgos

1 Kommentar
-1 ältere Kommentare anzeigen -1 ältere Kommentare ausblenden

Jan
Jan am 13 Jul. 2011
@Giorgros: Does fixing the mentioned line to "mexPrintf("sides[0] %f\n", sides[0])" solve the problem?

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Web Services finden Sie in Hilfe-Center und File Exchange

Tags

Gefragt:

Giorgos Kourakos
Giorgos Kourakos
am 12 Jul. 2011

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by