c语言调试过程中遇到的问题,程序调试中遇到的问题

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

#include

#include

#define PI 3.1415926

void main()

{

int i,j,k,l,m,n,p,q,u;//p,q为i点和j点的读取指针

int t,step=1;//定义时间步迭代变量，及总的迭代时间步数，

double L,B,H;//定义长宽高的迭代数

double Length=0.1,Broad=0.05,High=0.005,density=2400,dt=0.0000001;//定义材料的长宽高及密度,dt为时间间隔

double inx=0.05,iny=0.025,inz=0.005;//定义入射点坐标

double Grid=0.0005,Hor=0.0015,R=0.005,r;//定义晶格长度和近场范围半径,以及冲击物的半径,r为距离冲击物中心的距离

double Fr,sumf[3];//定义规定的面力密度,及各点的初始冲击力

double x,y,z,xj,yj,zj;//定义i点及j点的实际坐标形式,。

static double a[300000][3],v[300000][3],s[300000][3];//a为加速度，v为速度,s为位移

double K=20*pow(10,9),c,s0=0.005,vr,ks=pow(10,17);//K为体积模量，c为弹性常量,ks为力常数

double rp[3],rd[3],ap,s1,ad;//定义相对位置及相对位移,以及变形后的距离,s1为伸长率,ad为相对位置

static double sc[300000][3];

c=18*K/(PI*pow(Hor,4));

L=Length/Grid;B=Broad/Grid;H=High/Grid;

for(t=0;t<=step;t++)

{

p=0;

for(i=0;i<=(int)(L-1);i++)

for(j=0;j<=(int)(B-1);j++)

for(k=0;k<=(int)(H-1);k++)

{

x=i*Grid+Grid/2.0;y=j*Grid+Grid/2.0;z=k*Grid+Grid/2.0;//计算i点的实际坐标

sumf[0]=0;sumf[1]=0;sumf[2]=0;

if(t==0)

{

r=sqrt(pow(x-inx,2)+pow(y-iny,2)+pow(z-inz,2));//计算计算点离入射点的距离

if(r

Fr=-ks*pow(r-R,2);

else

Fr=0;

a[p][0]=0;a[p][1]=0;a[p][2]=Fr/(density*pow(Grid,3));

v[p][0]=0;v[p][1]=0;v[p][2]=a[p][2]*dt;

s[p][0]=0;s[p][1]=0;s[p][2]=v[p][2]*dt/2.0;

p++;

}

else

{

q=0;

for(l=0;l<=L-1;l++)

for(m=0;m<=B-1;m++)

for(n=0;n<=H-1;n++)

{

xj=l*Grid+Grid/2.0;yj=m*Grid+Grid/2.0;zj=n*Grid+Grid/2.0;//计算j点的实际左边

rp[0]=xj-x;rp[1]=yj-y;rp[2]=zj-z;//计算相对位置

rd[0]=s[q][0]-s[p][0];rd[1]=s[q][1]-s[p][1];rd[2]=s[q][2]-s[p][2];//计算相对位移

ap=sqrt(pow(rp[0]+rd[0],2)+pow(rp[1]+rd[1],2)+pow(rp[2]+rd[2],2));

ad=sqrt(rp[0]*rp[0]+rp[1]*rp[1]+rp[2]*rp[2]);

s1=(ap-ad)/ad;

if(s1

u=1;

else

u=0;

if(ad<=Hor-0.5*Grid)

vr=1;

else

if(ad>(Hor-0.5*Grid)&&ad<=(Hor+0.5*Grid))

vr=(Hor+0.5*Grid-ad)/Grid;

else

vr=0;

sumf[0]=sumf[0]+c*s1*u*(rd[0]+rp[0])*vr*pow(Grid,3)/ap;

sumf[1]=sumf[1]+c*s1*u*(rd[1]+rp[1])*vr*pow(Grid,3)/ap;

sumf[2]=sumf[2]+c*s1*u*(rd[2]+rp[2])*vr*pow(Grid,3)/ap;

q++;

}

a[p][0]=sumf[0]/density;a[p][1]=sumf[1]/density;a[p][2]=sumf[2]/density;

s[p][0]=s[p][0]+v[p][0]*dt+a[p][0]*dt*dt/2.0;s[p][1]=s[p][1]+v[p][1]*dt+a[p][1]*dt*dt/2.0;s[p][2]=s[p][2]+v[p][2]*dt+a[p][2]*dt*dt/2.0;

v[p][0]=v[p][0]+a[p][0]*dt;v[p][1]=v[p][1]+a[p][1]*dt;v[p][2]=v[p][2]+a[p][2]*dt;

p++;

}

}

}

p=0;

for(i=0;i<=L-1;i++)

for(j=0;j<=B-1;j++)

for(k=0;k<=H-1;k++)

{

x=i*Grid+Grid/2.0;y=j*Grid+Grid/2.0;z=k*Grid+Grid/2.0;

sc[p][0]=x+s[p][0];sc[p][1]=y+s[p][1];sc[p][2]=z+s[p][2];

printf("%.16lf\t%.16lf\t%.16lf\n",sc[p][0],sc[p][1],sc[p][2]);

p++;

}

}