C++编译时不报错 但是运行时显示停止工作 怎么检查错误在哪

[qq_32787919]

#pragma hdrstop

//---------------------------------------------------------------------------

#pragma argsused
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <conio.h>
#define PI 3.1415926535
float **space2d(int nr,int nc);
float ***space3d(int nr, int ny, int nc);
void free_space2d(float **a, int nr);
void free_space3d(float ***b,int nr,int ny);
void wfile(char filename[], float **data, int nr, int nc);
void wfile3d(char filename[],float ***data,int nr,int ny,int nc);
void create_model(float ***vp, float ***vs, float ***rho,float ***vf,float ***rhof, float ***lamda, float ***lamda2u, float ***mu,float***por,int nr,int ny, int nc);

float ***extmodel(float ***init_model,int nr,int ny,int nc,int np);


int main()
{
//给定参数

int NX=400; //x方向网格点数%%%%%
int NY=600; //y方向网格点数%%%%%
int NZ=800; //z方向网格点数%%%%
int NP=20; //pml层网格点数

int sx=NX/2+NP; //震源坐标点号
int sy=NY/2+NP;
int sz=NZ/2+NP;

int NX_ext=NX+2*NP; //
int NY_ext=NY+2*NP;
int NZ_ext=NZ+2*NP;

int NT=7000; //时间层数%%%%
// double H=1.0; //空间步长5%%%%
double H=0.005;
double RC=0.0001;
double DP=NP*H;
// double DT=0.0002; //时间步长
double DT=0.5*pow(10,-6); //时间步长%%%%
// double F0=30.0; //震源主频
double F0=3.0*pow(10,3); //震源主频,主频太高会造成频散
double T0=1.2/F0;
double Vpmax=4270.0; //模型最大纵波速度,用于稳定性计算
double Vpmin=1500.0; //模型最小纵波速度,用于控制数值频散
double Vsmax=2650.0;
double Vsmin=1500.0;

float ***vs; //初始模型横波速度
float ***vp; //初始模型纵波速度
float ***rho; //初始模型密度
float ***vf; //孔隙流体速度
float ***rhof; //孔隙流体密度
float ***mu; //剪切模量
float ***lamda; //拉梅系数
float ***lamda2u;
//float **Q;
//float **R;
float ***por;
// float **D11,**D12,**D22;
// float **b;
float **sis_x; //地震记录x分量 波形记录-二维数据
float **sis_y; //地震记录y分量
float **sis_z; //地震记录z分量


float ***vx_x;//速度分量
float ***vx_y;
float ***vx_z;
float ***vx;
float ***vy_x;
float ***vy_y;
float ***vy_z;
float ***vy;
float ***vz_x;
float ***vz_y;
float ***vz_z;
float ***vz;
/*float **Vx_x;//流相速度
float **Vx_z;
float **Vx;
float **Vz_x;
float **Vz_z;
float **Vz;
*/
float ***txx_x;//应力分量
float ***txx_y;
float ***txx_z;
float ***tyy_x;//应力分量
float ***tyy_y;
float ***tyy_z;

float ***tzz_x;
float ***tzz_y;
float ***tzz_z;

float ***txy_x;
float ***txy_y;
float ***txz_x;
float ***txz_z;
float ***tyz_y;
float ***tyz_z;
/*float **ssx;//流体应力
float **ssz;
float **dxi,**dxi2;
float **dzj,**dzj2;
*/
// double Kb,Ks,Kf,a,D;
double tt,x,y,z,xoleft,xoright,yoleft,yoright,zoleft,zoright,d0;
double v0,y5,y6,y7,y8,y9,y10,y11,y12,y13,rho_tempx,rho_tempy,rho_tempz,muxz,muxy,muyz,lamda2u_temp,lamda_temp;
// double D11_tempx,D12_tempx,D22_tempx,R_temp,D11_tempz,D12_tempz,D22_tempz;
double xx1,xy1,xz1,yx2,yy2,yz2,zx3,zy3,zz3,best_dt;
int ix,iy,iz,it;

vs=space3d(NZ,NY,NX);
vp=space3d(NZ,NY,NX);
rho=space3d(NZ,NY,NX);
vf=space3d(NZ,NY,NX);
rhof=space3d(NZ,NY,NX);
//Q=space2d(NZ,NX);
//R=space2d(NZ,NX);
por=space3d(NZ,NY,NX);
//D11=space2d(NZ,NX);D12=space2d(NZ,NX);D22=space2d(NZ,NX);
// b=space2d(NZ,NX);
mu=space3d(NZ,NY,NX);
lamda=space3d(NZ,NY,NX);
lamda2u=space3d(NZ,NY,NX);
create_model(vp, vs, rho, vf, rhof, lamda, lamda2u, mu,por, NZ,NY, NX);

char vs_name[]="vs_ext.dat";
float ***vs_ext;
float ***vp_ext;
float ***rho_ext;
float ***mu_ext;
float ***lamda_ext;
float ***lamda2u_ext;
// float ***Q_ext;
// float ***R_ext;
float ***vf_ext;
float ***rhof_ext;
float ***por_ext;
float ***dxi;
float ***dxi2;
float ***dyj;
float ***dyj2;
float ***dzk;
float ***dzk2;
//float ***b_ext;
// float **D11_ext,**D12_ext,**D22_ext;

//数据扩充，加上吸收边界
vs_ext=extmodel(vs,NZ,NY,NX,NP);
vp_ext=extmodel(vp,NZ,NY,NX,NP);
rho_ext=extmodel(rho,NZ,NY,NX,NP);
mu_ext=extmodel(mu,NZ,NY,NX,NP);
lamda_ext=extmodel(lamda,NZ,NY,NX,NP);
lamda2u_ext=extmodel(lamda2u,NZ,NY,NX,NP);
vf_ext=extmodel(vf,NZ,NY,NX,NP);
rhof_ext=extmodel(rhof,NZ,NY,NX,NP);
//Q_ext=extmodel(Q,NZ,NX,NP);
//R_ext=extmodel(R,NZ,NX,NP);
por_ext=extmodel(por,NZ,NY,NX,NP);
// b_ext=extmodel(b,NZ,NX,NP);
// D11_ext=extmodel(D11,NZ,NX,NP);D12_ext=extmodel(D12,NZ,NX,NP);
// D22_ext=extmodel(D22,NZ,NX,NP);

wfile3d(vs_name, vs_ext, NZ_ext, NY_ext, NX_ext); //建立文件

//申请空间
vx_x=space3d(NZ_ext,NY_ext,NX_ext);
vx_y=space3d(NZ_ext,NY_ext,NX_ext);
vx_z=space3d(NZ_ext,NY_ext,NX_ext);
vx=space3d(NZ_ext,NY_ext,NX_ext);

vy_x=space3d(NZ_ext,NY_ext,NX_ext);
vy_y=space3d(NZ_ext,NY_ext,NX_ext);
vy_z=space3d(NZ_ext,NY_ext,NX_ext);
vy=space3d(NZ_ext,NY_ext,NX_ext);

vz_x=space3d(NZ_ext,NY_ext,NX_ext);
vz_y=space3d(NZ_ext,NY_ext,NX_ext);
vz_z=space3d(NZ_ext,NY_ext,NX_ext);
vz=space3d(NZ_ext,NY_ext,NX_ext);
/*Vx_x=space2d(NZ_ext,NX_ext);
Vx_z=space2d(NZ_ext,NX_ext);
Vx=space2d(NZ_ext,NX_ext);
Vz_x=space2d(NZ_ext,NX_ext);
Vz_z=space2d(NZ_ext,NX_ext);
Vz=space2d(NZ_ext,NX_ext);*/
txx_x=space3d(NZ_ext,NY_ext,NX_ext);
txx_y=space3d(NZ_ext,NY_ext,NX_ext);
txx_z=space3d(NZ_ext,NY_ext,NX_ext);

tyy_x=space3d(NZ_ext,NY_ext,NX_ext);
tyy_y=space3d(NZ_ext,NY_ext,NX_ext);
tyy_z=space3d(NZ_ext,NY_ext,NX_ext);
tzz_x=space3d(NZ_ext,NY_ext,NX_ext);
tzz_y=space3d(NZ_ext,NY_ext,NX_ext);
tzz_z=space3d(NZ_ext,NY_ext,NX_ext);

txy_x=space3d(NZ_ext,NY_ext,NX_ext);
txy_y=space3d(NZ_ext,NY_ext,NX_ext);
txz_x=space3d(NZ_ext,NY_ext,NX_ext);
txz_z=space3d(NZ_ext,NY_ext,NX_ext);
tyz_y=space3d(NZ_ext,NY_ext,NX_ext);
tyz_z=space3d(NZ_ext,NY_ext,NX_ext);
/* tzz_x=space2d(NZ_ext,NX_ext);
tzz_z=space2d(NZ_ext,NX_ext);
txz_x=space2d(NZ_ext,NX_ext);
txz_z=space2d(NZ_ext,NX_ext);
ssx=space2d(NZ_ext,NX_ext);
ssz=space2d(NZ_ext,NX_ext);*/

dxi=space3d(NZ_ext,NY_ext,NX_ext);
dxi2=space3d(NZ_ext,NY_ext,NX_ext);
dyj=space3d(NZ_ext,NY_ext,NX_ext);
dyj2=space3d(NZ_ext,NY_ext,NX_ext);
dzk=space3d(NZ_ext,NY_ext,NX_ext);
dzk2=space3d(NZ_ext,NY_ext,NX_ext);

//(xoleft,xoright)模型的物理边界—不算PML
xoleft=DP;
xoright=(NX_ext-1)*H-DP;
yoleft=DP;
yoright=(NY_ext-1)*H-DP;
zoleft=DP;
zoright=(NZ_ext-1)*H-DP;

//用于对vx_x[iz][ix]等的求解，加入吸收边界，使数值衰减
for(iz=0; iz<NZ_ext; iz++)
{
for(iy=0;iy<NY_ext;iy++)
{
for (ix=0; ix<NX_ext; ix++)
{
x=ix*H;
if (x < xoleft)
{
v0=vp_ext[iz][iy][0];
d0=3.0*v0*log(1.0/RC)/(2.0*DP);
dxi[iz][iy][ix]=d0*pow(((xoleft-x)/DP),2);
dxi2[iz][iy][ix]=d0*pow(((xoleft-x-0.5*H)/DP),2);
}
else if (x >=0.9999*xoright )
{
v0=vp_ext[iz][iy][NX_ext-1];
d0=3.0*v0*log(1.0/RC)/(2.0*DP);
dxi[iz][iy][ix]=d0*pow(((x-xoright)/DP),2);
dxi2[iz][iy][ix]=d0*pow(((x+0.5*H-xoright)/DP),2);
}
else
{
dxi[iz][iy][ix]=0.; dxi2[iz][iy][ix]=0.;
}
}
}
}

for (iz=0; iz<NZ_ext; iz++)
{
for(iy=0;iy<NY_ext;iy++)
{
y=iy*H;
for(ix=0; ix<NX_ext; ix++)
{
if (y < yoleft)
{
v0=vp_ext[iz][0][ix];
d0=3.0*v0*log(1.0/RC)/(2.0*DP);
//d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H;
dyj[iz][iy][ix]=d0*pow(((yoleft-y)/DP),2);
dyj2[iz][iy][ix]=d0*pow(((yoleft-y-0.5*H)/DP),2);
}
else if (y >=0.9999*yoright )
{
v0=vp_ext[iz][NY_ext-1][ix];
d0=3.0*v0*log(1.0/RC)/(2.0*DP);
dyj[iz][iy][ix]=d0*pow(((y-yoright)/DP),2);
dyj2[iz][iy][ix]=d0*pow(((y+0.5*H-yoright)/DP),2);
}
else
{
dyj[iz][iy][ix]=0.; dyj2[iz][iy][ix]=0.;
}
}
}
}
for (iz=0; iz<NZ_ext; iz++)
{
z=iz*H;
for(iy=0;iy<NY_ext;iy++)
{

for(ix=0; ix<NX_ext; ix++)
{
if (z < zoleft)
{
v0=vp_ext[0][iy][ix];
d0=3.0*v0*log(1.0/RC)/(2.0*DP);
//d0=3.0*v0*(8.0/15.0-3.0/100.0*NP+NP*NP/1500.0)/H;
dzk[iz][iy][ix]=d0*pow(((zoleft-z)/DP),2);
dzk2[iz][iy][ix]=d0*pow(((zoleft-z-0.5*H)/DP),2);
}
else if (z >=0.9999*zoright )
{
v0=vp_ext[NZ_ext-1][iy][ix];
d0=3.0*v0*log(1.0/RC)/(2.0*DP);
dzk[iz][iy][ix]=d0*pow(((z-zoright)/DP),2);
dzk2[iz][iy][ix]=d0*pow(((z+0.5*H-zoright)/DP),2);
}
else
{
dzk[iz][iy][ix]=0.; dzk2[iz][iy][ix]=0.;
}
}
}
}

//开始时间递推
double a1=9.0/8.0;
double a2= -1.0/24.0;

//这里的代码应该修正一下，应该从数据里面提取最大最小速度

//控制网格频散
if (Vsmin / (F0*H) < 15)
printf("空间步长太大,可能引起明显的网格频散\n");

//检查稳定性条件
best_dt=6.0*H/(7.0*sqrt(2.0)*Vpmax);
if (DT >= best_dt)
printf("时间步长过大，应该小于 %f\n",best_dt);

clock_t start,finish;
start=clock();

sis_x=space2d(NZ,NT);
sis_y=space2d(NZ,NT);
sis_z=space2d(NZ,NT);

for (it=0; it<NT; it++)
{
tt=it*DT;
if(it % 100 == 0) printf("it=%d\n",it);
//计算质点速度v
for (iz=2; iz<NZ_ext-2; iz++)
for(iy=2; iy<NY_ext-2; iy++)
for (ix=2; ix<NX_ext-2; ix++)
{
rho_tempx=2/(rho_ext[iz][iy][ix]+rho_ext[iz][iy][ix+1]);
rho_tempy=2/(rho_ext[iz][iy][ix]+rho_ext[iz][iy+1][ix]);
rho_tempz=2/(rho_ext[iz][iy][ix]+rho_ext[iz+1][iy][ix]);


xx1=(a1*(txx_x[iz][iy][ix+1]-txx_x[iz][iy][ix]+txx_y[iz][iy][ix+1]-txx_y[iz][iy][ix]+txx_z[iz][iy][ix+1]-txx_z[iz][iy][ix])+
a2*(txx_x[iz][iy][ix+2]-txx_x[iz][iy][ix-1]+txx_y[iz][iy][ix+2]-txx_y[iz][iy][ix-1]+txx_z[iz][iy][ix+2]-txx_z[iz][iy][ix-1]));
xy1=(a1*(txy_x[iz][iy][ix]-txy_x[iz][iy-1][ix]+txy_y[iz][iy][ix]-txy_y[iz][iy-1][ix])+
a2*(txy_x[iz][iy+1][ix]-txy_x[iz][iy-2][ix]+txy_y[iz][iy+1][ix]-txy_y[iz][iy-2][ix]));
xz1=(a1*(txz_x[iz][iy][ix]-txz_x[iz-1][iy][ix]+txz_z[iz][iy][ix]-txz_z[iz-1][iy][ix])+
a2*(txz_x[iz+1][iy][ix]-txz_x[iz-2][iy][ix]+txz_z[iz+1][iy][ix]-txz_z[iz-2][iy][ix]));

yx2=(a1*(txy_x[iz][iy][ix]-txy_x[iz][iy][ix-1]+txy_y[iz][iy][ix]-txy_y[iz][iy][ix-1])+
a2*(txy_x[iz][iy][ix+1]-txy_x[iz][iy][ix-2]+txy_y[iz][iy][ix+1]-txy_y[iz][iy][ix-2]));
yy2=(a1*(tyy_x[iz][iy+1][ix]-tyy_x[iz][iy][ix]+tyy_y[iz][iy+1][ix]-tyy_y[iz][iy][ix]+tyy_z[iz][iy+1][ix]-tyy_z[iz][iy][ix])+
a2*(tyy_x[iz][iy+2][ix]-tyy_x[iz][iy-1][ix]+tyy_y[iz][iy+2][ix]-tyy_y[iz][iy-1][ix]+tyy_z[iz][iy+2][ix]-tyy_z[iz][iy-1][ix]));
yz2=(a1*(tyz_y[iz][iy][ix]-tyz_y[iz-1][iy][ix]+tyz_z[iz][iy][ix]-tyz_z[iz-1][iy][ix])+
a2*(tyz_y[iz+1][iy][ix]-tyz_y[iz-2][iy][ix]+tyz_z[iz+1][iy][ix]-tyz_z[iz-2][iy][ix]));

zx3=(a1*(txz_x[iz][iy][ix]-txz_x[iz][iy][ix-1]+txz_z[iz][iy][ix]-txz_z[iz][iy][ix-1])+
a2*(txz_x[iz][iy][ix+1]-txz_x[iz][iy][ix-2]+txz_z[iz][iy][ix+1]-txz_z[iz][iy][ix-2]));
zy3=(a1*(tyz_y[iz][iy][ix]-tyz_y[iz][iy-1][ix]+tyz_z[iz][iy][ix]-tyz_z[iz][iy-1][ix])+
a2*(tyz_y[iz][iy+1][ix]-tyz_y[iz][iy-2][ix]+tyz_z[iz][iy+1][ix]-tyz_z[iz][iy-2][ix]));
zz3=(a1*(tzz_x[iz+1][iy][ix]-tzz_x[iz][iy][ix]+tzz_y[iz+1][iy][ix]-tzz_y[iz][iy][ix]+tzz_z[iz+1][iy][ix]-tzz_z[iz][iy][ix])+
a2*(tzz_x[iz+2][iy][ix]-tzz_x[iz-1

[qq_32787919]

vx_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi2[iz][iy][ix]))*((1-0.5*DT*dxi2[iz][iy][ix])*vx_x[iz][iy][ix]+((DT*rho_tempx/H)*xx1));
vx_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj[iz][iy][ix]))*((1-0.5*DT*dyj[iz][iy][ix])*vx_y[iz][iy][ix]+((DT*rho_tempy/H)*xy1));
vx_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk[iz][iy][ix]))*((1-0.5*DT*dzk[iz][iy][ix])*vx_z[iz][iy][ix]+((DT*rho_tempz/H)*xz1));

vy_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi[iz][iy][ix]))*((1-0.5*DT*dxi[iz][iy][ix])*vy_x[iz][iy][ix]+((DT*rho_tempx/H)*yx2));
vy_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj2[iz][iy][ix]))*((1-0.5*DT*dyj2[iz][iy][ix])*vy_y[iz][iy][ix]+((DT*rho_tempy/H)*yy2));
vy_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk[iz][iy][ix]))*((1-0.5*DT*dzk[iz][iy][ix])*vy_z[iz][iy][ix]+((DT*rho_tempz/H)*yz2));

vz_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi[iz][iy][ix]))*((1-0.5*DT*dxi[iz][iy][ix])*vz_x[iz][iy][ix]+((DT*rho_tempx/H)*zx3));
vz_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj[iz][iy][ix]))*((1-0.5*DT*dyj[iz][iy][ix])*vz_y[iz][iy][ix]+((DT*rho_tempy/H)*zy3));
vz_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk2[iz][iy][ix]))*((1-0.5*DT*dzk2[iz][iy][ix])*vz_z[iz][iy][ix]+((DT*rho_tempz/H)*zz3));

}
//计算应力
for (iz=2; iz<NZ_ext-2; iz++)
for(iy=2; iy<NY_ext-2; iy++)
for (ix=2; ix<NX_ext-2; ix++)
{

if(mu_ext[iz][iy][ix]==0||mu_ext[iz][iy][ix+1]==0||mu_ext[iz+1][iy][ix]==0||mu_ext[iz+1][iy][ix+1]==0)
muxz=0.0;
else
muxz=1/(0.25*(1/mu_ext[iz][iy][ix]+1/mu_ext[iz][iy][ix+1]+1/mu_ext[iz+1][iy][ix]+1/mu_ext[iz+1][iy][ix+1]));
if(mu_ext[iz][iy][ix]==0||mu_ext[iz][iy+1][ix]==0||mu_ext[iz+1][iy][ix]==0||mu_ext[iz+1][iy+1][ix]==0)
muyz=0.0;
else
muyz=1/(0.25*(1/mu_ext[iz][iy][ix]+1/mu_ext[iz][iy+1][ix]+1/mu_ext[iz+1][iy][ix]+1/mu_ext[iz+1][iy+1][ix]));
if(mu_ext[iz][iy][ix]==0||mu_ext[iz][iy][ix+1]==0||mu_ext[iz][iy+1][ix]==0||mu_ext[iz][iy+1][ix+1]==0)
muxy=0.0;
else
muxy=1/(0.25*(1/mu_ext[iz][iy][ix]+1/mu_ext[iz][iy][ix+1]+1/mu_ext[iz][iy+1][ix]+1/mu_ext[iz][iy+1][ix+1]));

lamda2u_temp=lamda2u_ext[iz][iy][ix];
lamda_temp=lamda_ext[iz][iy][ix];
// muxz=0.0;
// Qz=Q_ext[iz][ix];R_temp=R_ext[iz][ix];lamda2u_temp=lamda2u_ext[iz][ix];lamda_temp=lamda_ext[iz][ix];

y5=a1*(vx_x[iz][iy][ix]-vx_x[iz][iy][ix-1]+vx_y[iz][iy][ix]-vx_y[iz][iy][ix-1]+vx_z[iz][iy][ix]-vx_z[iz][iy][ix-1])+
a2*(vx_x[iz][iy][ix+1]-vx_x[iz][iy][ix-2]+vx_y[iz][iy][ix+1]-vx_y[iz][iy][ix-2]+vx_z[iz][iy][ix+1]-vx_z[iz][iy][ix-2]);

y6=a1*(vy_x[iz][iy][ix]-vy_x[iz][iy-1][ix]+vy_y[iz][iy][ix]-vy_y[iz][iy-1][ix]+vy_z[iz][iy][ix]-vy_z[iz][iy-1][ix])+
a2*(vy_x[iz][iy+1][ix]-vy_x[iz][iy-2][ix]+vy_y[iz][iy+1][ix]-vy_y[iz][iy-2][ix]+vy_z[iz][iy+1][ix]-vy_z[iz][iy-2][ix]);

y7=a1*(vz_x[iz][iy][ix]-vz_x[iz-1][iy][ix]+vz_y[iz][iy][ix]-vz_y[iz-1][iy][ix]+vz_z[iz][iy][ix]-vz_z[iz-1][iy][ix])+
a2*(vz_x[iz+1][iy][ix]-vz_x[iz-2][iy][ix]+vz_y[iz+1][iy][ix]-vz_y[iz-2][iy][ix]+vz_z[iz+1][iy][ix]-vz_z[iz-2][iy][ix]);

y8=(muxy*DT/H)*(a1*(vx_x[iz][iy+1][ix]-vx_x[iz][iy][ix]+vx_y[iz][iy+1][ix]-vx_y[iz][iy][ix]+vx_z[iz][iy+1][ix]-vx_z[iz][iy][ix])+
a2*(vx_x[iz][iy+2][ix]-vx_x[iz][iy-1][ix]+vx_y[iz][iy+2][ix]-vx_y[iz][iy-1][ix]+vx_z[iz][iy+2][ix]-vx_z[iz][iy-1][ix]));

y9=(muxy*DT/H)*(a1*(vy_x[iz][iy][ix+1]-vy_x[iz][iy][ix]+vy_y[iz][iy][ix+1]-vy_y[iz][iy][ix]+vy_z[iz][iy][ix+1]-vy_z[iz][iy][ix])+
a2*(vy_x[iz][iy][ix+2]-vy_x[iz][iy][ix-1]+vy_y[iz][iy][ix+2]-vy_y[iz][iy][ix-1]+vy_z[iz][iy][ix+2]-vy_z[iz][iy][ix-1]));

y10=(muxz*DT/H)*(a1*(vx_x[iz+1][iy][ix]-vx_x[iz][iy][ix]+vx_y[iz+1][iy][ix]-vx_y[iz][iy][ix]+vx_z[iz+1][iy][ix]-vx_z[iz][iy][ix])+
a2*(vx_x[iz+2][iy][ix]-vx_x[iz-1][iy][ix]+vx_y[iz+2][iy][ix]-vx_y[iz-1][iy][ix]+vx_z[iz+2][iy][ix]-vx_z[iz-1][iy][ix]));

y11=(muxz*DT/H)*(a1*(vz_x[iz][iy][ix+1]-vz_x[iz][iy][ix]+vz_y[iz][iy][ix+1]-vz_y[iz][iy][ix]+vz_z[iz][iy][ix+1]-vz_z[iz][iy][ix])+
a2*(vz_x[iz][iy][ix+2]-vz_x[iz][iy][ix-1]+vz_y[iz][iy][ix+2]-vz_y[iz][iy][ix-1]+vz_z[iz][iy][ix+2]-vz_z[iz][iy][ix-1]));

y12=(muyz*DT/H)*(a1*(vy_x[iz+1][iy][ix]-vy_x[iz][iy][ix]+vy_y[iz+1][iy][ix]-vy_y[iz][iy][ix]+vy_z[iz+1][iy][ix]-vy_z[iz][iy][ix])+
a2*(vy_x[iz+2][iy][ix]-vy_x[iz-1][iy][ix]+vy_y[iz+2][iy][ix]-vy_y[iz-1][iy][ix]+vy_z[iz+2][iy][ix]-vy_z[iz-1][iy][ix]));

y13=(muyz*DT/H)*(a1*(vz_x[iz][iy+1][ix]-vz_x[iz][iy][ix]+vz_y[iz][iy+1][ix]-vz_y[iz][iy][ix]+vz_z[iz][iy+1][ix]-vz_z[iz][iy][ix])+
a2*(vz_x[iz][iy+2][ix]-vz_x[iz][iy-1][ix]+vz_y[iz][iy+2][ix]-vz_y[iz][iy-1][ix]+vz_z[iz][iy+2][ix]-vz_z[iz][iy-1][ix]));


txx_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi[iz][iy][ix]))*((1-0.5*DT*dxi[iz][iy][ix])*txx_x[iz][iy][ix]+(lamda2u_temp*DT/H)*y5);
txx_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj[iz][iy][ix]))*((1-0.5*DT*dyj[iz][iy][ix])*txx_y[iz][iy][ix]+(lamda_temp*DT/H)*y6);

txx_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk[iz][iy][ix]))*((1-0.5*DT*dzk[iz][iy][ix])*txx_z[iz][iy][ix]+(lamda_temp*DT/H)*y7);

tyy_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi[iz][iy][ix]))*((1-0.5*DT*dxi[iz][iy][ix])*tyy_x[iz][iy][ix]+(lamda_temp*DT/H)*y5);
tyy_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj[iz][iy][ix]))*((1-0.5*DT*dyj[iz][iy][ix])*tyy_y[iz][iy][ix]+(lamda2u_temp*DT/H)*y6);

tzz_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk[iz][iy][ix]))*((1-0.5*DT*dzk[iz][iy][ix])*tyy_z[iz][iy][ix]+(lamda_temp*DT/H)*y7);

tzz_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi[iz][iy][ix]))*((1-0.5*DT*dxi[iz][iy][ix])*tzz_x[iz][iy][ix]+(lamda_temp*DT/H)*y5);
tzz_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj[iz][iy][ix]))*((1-0.5*DT*dyj[iz][iy][ix])*tzz_y[iz][iy][ix]+(lamda_temp*DT/H)*y6);

tzz_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk[iz][iy][ix]))*((1-0.5*DT*dzk[iz][iy][ix])*tzz_z[iz][iy][ix]+(lamda2u_temp*DT/H)*y7);

txy_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi2[iz][iy][ix]))*((1-0.5*DT*dxi2[iz][iy][ix])*txy_x[iz][iy][ix]+y9);
txy_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj2[iz][iy][ix]))*((1-0.5*DT*dyj2[iz][iy][ix])*txy_y[iz][iy][ix]+y8);
txz_x[iz][iy][ix]=(1.0/(1.0+0.5*DT*dxi2[iz][iy][ix]))*((1-0.5*DT*dxi2[iz][iy][ix])*txz_x[iz][iy][ix]+y11);
txz_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk2[iz][iy][ix]))*((1-0.5*DT*dzk2[iz][iy][ix])*txz_z[iz][iy][ix]+y10);
tyz_y[iz][iy][ix]=(1.0/(1.0+0.5*DT*dyj2[iz][iy][ix]))*((1-0.5*DT*dyj2[iz][iy][ix])*tyz_y[iz][iy][ix]+y13);
tyz_z[iz][iy][ix]=(1.0/(1.0+0.5*DT*dzk2[iz][iy][ix]))*((1-0.5*DT*dzk2[iz][iy][ix])*tyz_z[iz][iy][ix]+y12);
}
//加震源
tzz_z[sz][sy][sx]=tzz_z[sz][sy][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0)); //对正应力加爆炸源
tyy_y[sz][sy][sx]=tyy_y[sz][sy][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0));
txx_x[sz][sy][sx]=txx_x[sz][sy][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0));
// tzz_z[sz][sx]=tzz_z[sz][sx]-2*PI*PI*F0*F0*(tt-T0)*exp(-PI*PI*F0*F0*(tt-T0)*(tt-T0));
for (iz=0; iz<NZ_ext; iz++)
for(iy=0;iy<NY_ext;iy++)
{
vx_x[iz][iy][0]=0.; vx_y[iz][iy][0]=0.;vx_z[iz][iy][0]=0.;vx_x[iz][iy][NX_ext-1]=0.; vx_y[iz][iy][NX_ext-1]=0.;vx_z[iz][iy][NX_ext-1]=0.;
vy_x[iz][iy][0]=0.; vy_y[iz][iy][0]=0.;vy_z[iz][iy][0]=0.;vy_x[iz][iy][NX_ext-1]=0.; vy_y[iz][iy][NX_ext-1]=0.;vy_z[iz][iy][NX_ext-1]=0.;
vz_x[iz][iy][0]=0.; vz_y[iz][iy][0]=0.;vz_z[iz][iy][0]=0.;vz_x[iz][iy][NX_ext-1]=0.; vz_y[iz][iy][NX_ext-1]=0.;vz_z[iz][iy][NX_ext-1]=0.;
vx_x[iz][iy][1]=0.; vx_y[iz][iy][1]=0.;vx_z[iz][iy][1]=0.;vx_x[iz][iy][NX_ext-2]=0.; vx_y[iz][iy][NX_ext-2]=0.;vx_z[iz][iy][NX_ext-2]=0.;
vy_x[iz][iy][1]=0.; vy_y[iz][iy][1]=0.;vy_z[iz][iy][1]=0.;vy_x[iz][iy][NX_ext-2]=0.; vy_y[iz][iy][NX_ext-2]=0.;vy_z[iz][iy][NX_ext-2]=0.;
vz_x[iz][iy][1]=0.; vz_y[iz][iy][1]=0.;vz_z[iz][iy][1]=0.;vz_x[iz][iy][NX_ext-2]=0.; vz_y[iz][iy][NX_ext-2]=0.;vz_z[iz][iy][NX_ext-2]=0.;
}
for (ix=0; ix<NX_ext; ix++)
for(iy=0;iy<NY_ext;iy++)
{
vx_x[NZ_ext-1][iy][ix]=0.;vx_y[NZ_ext-1][iy][ix]=0.;vx_z[NZ_ext-1][iy][ix]=0.;
vy_x[NZ_ext-1][iy][ix]=0.;vy_y[NZ_ext-1][iy][ix]=0.;vy_z[NZ_ext-1][iy][ix]=0.;
vz_x[NZ_ext-1][iy][ix]=0.;vz_y[NZ_ext-1][iy][ix]=0.;vz_z[NZ_ext-1][iy][ix]=0.;

vx_x[0][iy][ix]=0.;vx_y[0][iy][ix]=0.;vx_z[0][iy][ix]=0.;
vy_x[0][iy][ix]=0.;vy_y[0][iy][ix]=0.;vy_z[0][iy][ix]=0.;
vz_x[0][iy][ix]=0.;vz_y[0][iy][ix]=0.;vz_z[0][iy][ix]=0.;
vx_x[NZ_ext-2][iy][ix]=0.;vx_y[NZ_ext-2][iy][ix]=0.;vx_z[NZ_ext-2][iy][ix]=0.;
vy_x[NZ_ext-2][iy][ix]=0.;vy_y[NZ_ext-2][iy][ix]=0.;vy_z[NZ_ext-2][iy][ix]=0.;
vz_x[NZ_ext-2][iy][ix]=0.;vz_y[NZ_ext-2][iy][ix]=0.;vz_z[NZ_ext-2][iy][ix]=0.;

vx_x[1][iy][ix]=0.;vx_y[1][iy][ix]=0.;vx_z[1][iy][ix]=0.;
vy_x[1][iy][ix]=0.;vy_y[1][iy][ix]=0.;vy_z[1][iy][ix]=0.;
vz_x[1][iy][ix]=0.;vz_y[1][iy][ix]=0.;vz_z[1][iy][ix]=0.;
}
for (iz=0; iz<NZ_ext; iz++)
for(ix=0;ix<NX_ext;ix++)
{
vx_x[iz][NY_ext-1][ix]=0.;vx_y[iz][NY_ext-1][ix]=0.;vx_z[iz][NY_ext-1][ix]=0.;
vy_x[iz][NY_ext-1][ix]=0.;vy_y[iz][NY_ext-1][ix]=0.;vy_z[iz][NY_ext-1][ix]=0.;
vz_x[iz][NY_ext-1][ix]=0.;vz_y[iz][NY_ext-1][ix]=0.;vz_z[iz][NY_ext-1][ix]=0.;

vx_x[iz][0][ix]=0.;vx_y[iz][0][ix]=0.;vx_z[iz][0][ix]=0.;
vy_x[iz][0][ix]=0.;vy_y[iz][0][ix]=0.;vy_z[iz][0][ix]=0.;
vz_x[iz][0][ix]=0.;vz_y[iz][0][ix]=0.;vz_z[iz][0][ix]=0.;
vx_x[iz][NY_ext-2][ix]=0.;vx_y[iz][NY_ext-2][ix]=0.;vx_z[iz][NY_ext-2][ix]=0.;
vy_x[iz][NY_ext-2][ix]=0.;vy_y[iz][NY_ext-2][ix]=0.;vy_z[iz][NY_ext-2][ix]=0.;
vz_x[iz][NY_ext-2][ix]=0.;vz_y[iz][NY_ext-2][ix]=0.;vz_z[iz][NY_ext-2][ix]=0.;

vx_x[iz][1][ix]=0.;vx_y[iz][1][ix]=0.;vx_z[iz][1][ix]=0.;
vy_x[iz][1][ix]=0.;vy_y[iz][1][ix]=0.;vy_z[iz][1][ix]=0.;
vz_x[iz][1][ix]=0.;vz_y[iz][1][ix]=0.;vz_z[iz][1][ix]=0.;
}
//记录地震记录
for (iz=NP; iz<NP+NZ; iz++)
{
sis_z[iz-NP][it]=vz_x[iz][sy][sx]+vz_y[iz][sy][sx]+vz_z[iz][sy][sx];
sis_y[iz-NP][it]=vy_x[iz][sy][sx]+vy_y[iz][sy][sx]+vy_z[iz][sy][sx];
sis_x[iz-NP][it]=vx_x[iz][sy][sx]+vx_y[iz][sy][sx]+vx_z[iz][sy][sx];
//sis_z[iz-NP][it]=vz_x[iz][NX_ext-iz]+vz_z[iz][NX_ext-iz];
// sis_x[iz-NP][it]=vx_x[iz][NX_ext-iz]+vx_z[iz][NX_ext-iz];
}
}
finish=clock();
printf("%f seconds\n",(double)(finish-start)/CLOCKS_PER_SEC);
//输出波场快照
char vzname[]="vz.dat";
char vxname[]="vx.dat";
char vyname[]="vy.dat";
char txx_xname[]="txx_x.dat";
for (iz=0; iz<NZ_ext; iz++)
{
for (iy=0; iy<NY_ext; iy++)
for (ix=0; ix<NX_ext; ix++)
{
vx[iz][iy][ix]=vx_x[iz][iy][ix]+vx_y[iz][iy][ix]+vx_z[iz][iy][ix];
vy[iz][iy][ix]=vy_x[iz][iy][ix]+vy_y[iz][iy][ix]+vy_z[iz][iy][ix];
vz[iz][iy][ix]=vz_x[iz][iy][ix]+vz_y[iz][iy][ix]+vz_z[iz][iy][ix];

}
}

[qq_32787919]

printf("NZ_ext = %d NY_ext = %d NX_ext = %d\n",NZ_ext,NY_ext,NX_ext);
wfile3d(vxname, vx, NZ_ext, NY_ext, NX_ext);
wfile3d(vzname, vz, NZ_ext, NY_ext, NX_ext);
wfile3d(vyname, vy, NZ_ext, NY_ext, NX_ext);
//wfile(ssxname, ssx, NZ_ext, NX_ext);
//输出地震记录
char sisxname[]="sisx.dat";
char sisyname[]="sisy.dat";
char siszname[]="sisz.dat";
printf("NZ = %d NT = %d\n",NZ,NT);
wfile(sisxname, sis_x, NZ, NT);
wfile(sisyname, sis_y, NZ, NT);
wfile(siszname, sis_z, NZ, NT);
free_space3d(tzz_x,NZ_ext,NY_ext);
free_space3d(tzz_y,NZ_ext,NY_ext);
free_space3d(tzz_z,NZ_ext,NY_ext);
free_space3d(tyy_x,NZ_ext,NY_ext);
free_space3d(tyy_y,NZ_ext,NY_ext);
free_space3d(tyy_z,NZ_ext,NY_ext);
free_space3d(txx_x,NZ_ext,NY_ext);
free_space3d(txx_y,NZ_ext,NY_ext);
free_space3d(txx_z,NZ_ext,NY_ext);
free_space3d(txy_x,NZ_ext,NY_ext);
free_space3d(txy_y,NZ_ext,NY_ext);
free_space3d(txz_x,NZ_ext,NY_ext);
free_space3d(txz_z,NZ_ext,NY_ext);
free_space3d(tyz_y,NZ_ext,NY_ext);
free_space3d(tyz_z,NZ_ext,NY_ext);
free_space3d(vx_x,NZ_ext,NY_ext);
free_space3d(vx_y,NZ_ext,NY_ext);
free_space3d(vx_z,NZ_ext,NY_ext);
free_space3d(vy_x,NZ_ext,NY_ext);
free_space3d(vy_y,NZ_ext,NY_ext);
free_space3d(vy_z,NZ_ext,NY_ext);
free_space3d(vz_x,NZ_ext,NY_ext);
free_space3d(vz_y,NZ_ext,NY_ext);
free_space3d(vz_z,NZ_ext,NY_ext);

free_space3d(vs_ext,NZ_ext,NY_ext);
free_space3d(vp_ext,NZ_ext,NY_ext);
free_space3d(rho_ext,NZ_ext,NY_ext);

free_space3d(vs,NZ,NY);
free_space3d(vp,NZ,NY);
free_space3d(mu,NZ,NY);
free_space3d(mu_ext,NZ_ext,NY_ext);
free_space3d(lamda,NZ,NY);
free_space3d(lamda2u,NZ,NY);
free_space3d(lamda_ext,NZ_ext,NY_ext);
free_space3d(lamda2u_ext,NZ_ext,NY_ext);
free_space3d(rho,NZ,NY);

free_space3d(dxi,NZ_ext,NY_ext);
free_space3d(dxi2,NZ_ext,NY_ext);
free_space3d(dyj,NZ_ext,NY_ext);
free_space3d(dyj2,NZ_ext,NY_ext);
free_space3d(dzk,NZ_ext,NY_ext);
free_space3d(dzk2,NZ_ext,NY_ext);


free_space2d(sis_x,NZ);
free_space2d(sis_z,NZ);
printf("\\Press any key to exit program...");
getch();
return 0;
}

//申请二维动态数组
float **space2d(int nr,int nc)
{

float **a;
int i;
a=(float **)calloc(nr,sizeof(float *));
for (i=0; i<nr; i++)
a[i]=(float *)calloc(nc,sizeof(float));

return a;
}
//申请三维动态数组
float ***space3d(int nr,int ny,int nc)
{
float ***b;
int i,j;
b=(float ***)malloc(sizeof(float **) *nr);
for(i=0;i<nr;i++)
{
b[i]=(float **)malloc(sizeof(float *) *ny);
}
for(i=0;i<nr;i++)
{
for(j=0;j<ny;j++)
{
b[i][j]=(float *)malloc(sizeof(float) nc);
}}
return b;
}
//释放二维动态数组
void free_space2d(float **a,int nr)
{
int i;
for (i=0; i<nr; i++)
free(a[i]);
free(a);
}
//释放三维动态数组
void free_space3d(float ***b,int nr,int ny)
{
int i,j;
for( i=0;i<nr;i++)
{
for( j=0;j<ny;j++)
{
free(b[i][j]);
}
}
for(i=0;i<nr;i++)
{
free(b[i]);
}
free(b);
}

//将二进制数据写入文件—二维

void wfile(char filename[],float **data,int nr,int nc)
{
int i,j;
FILE *fp=fopen(filename,"wt");
/*
for(int i=0;i<nr;i++)
{
fwrite(data[i],sizeof(float),nc,fp);
}
*/

for (i=0; i<nr; i++)
{
for (j=0; j<nc; j++)
{
fprintf(fp,"%e ",data[i][j]);
if((j+1) % nc == 0)
fprintf(fp,"\n");
}
fprintf(fp,"\n");
}
// fwrite(&data[i][j],1,sizeof(float),fp);
fclose(fp);
}

//将二进制数据写入文件—三维

void wfile3d(char filename[],float ***data,int nr,int ny,int nc)
{
int i,j,k;
FILE *fp=fopen(filename,"wt");
/*
for(int i=0;i<nr;i++)
{
fwrite(data[i],sizeof(float),nc,fp);
}
*/

for (i=0; i<nr; i++)
{
for (j=0; j<ny; j++)
{
for(k=0;k<nc;k++)
{
fprintf(fp,"%e ",data[i][j][k]);
if((j+1) % ny == 0||(k+1) % nc == 0)
fprintf(fp,"\n");
}
fprintf(fp,"\n");
}
}
// fwrite(&data[i][j],1,sizeof(float),fp);
fclose(fp);
}
//建立地质模型
void create_model(float ***vp, float ***vs, float ***rho,float ***vf,float ***rhof, float ***lamda, float ***lamda2u, float ***mu,float ***por,int nr,int ny, int nc)
{
//这里进行修改，最好是可视化设计或从文件读入
int ix,iy,iz;
// double por1,Ks1,Kb1,Kf1,a1,D1,tao,rou11,rou12,rou22,eta;
// eta=58.0;
for (iz=0; iz<nr; iz++)
{
for(iy=0;iy<ny;iy++)
{
for (ix=0; ix<nc; ix++)
{
if(ix>=(nc/2-10)&&ix<(nc/2+10)&&iy>=(ny/2-10)&&iy<(ny/2+10))//竖直井孔%%%%%%
// if(iz<(nr-14-ix)||iz>=(nr+14-ix))//45度倾斜井孔
{
vp[iz][iy][ix]=1600.0;
vs[iz][iy][ix]=0.0;
vf[iz][iy][ix]=1600.0;
rho[iz][iy][ix]=1000.0;
rhof[iz][iy][ix]=1000.0;
por[iz][iy][ix]=1.0;
// b[iz][ix]=0.0;

mu[iz][iy][ix]=0.0;
lamda2u[iz][iy][ix]=rho[iz][iy][ix]*vp[iz][iy][ix]*vp[iz][iy][ix];
lamda[iz][iy][ix]=lamda2u[iz][iy][ix];
/* Kb1=rho[iz][ix]*(1-por[iz][ix])*(vp[iz][ix]*vp[iz][ix]-vs[iz][ix]*vs[iz][ix]*4.0/3.0);//骨架压缩模量，声波测井原理与应用，P39
Ks1=rho[iz][ix]*(vp[iz][ix]*vp[iz][ix]-vs[iz][ix]*vs[iz][ix]*4.0/3.0);//岩石固态颗粒的体积模量
Kf1=rhof[iz][ix]*vf[iz][ix]*vf[iz][ix];//孔隙流体的体积压缩模量
a1=1-Kb1/Ks1;
D1=a1-por[iz][ix]+por[iz][ix]*Ks1/Kf1;
Q[iz][ix]=(a1-por[iz][ix])*por[iz][ix]*Ks1/D1;
R[iz][ix]=por[iz][ix]*por[iz][ix]*Ks1/D1;
tao=0.5*(1+1/por[iz][ix]);//孔隙弯曲度
rou11=(1-por[iz][ix])*rho[iz][ix]-(1-tao)*por[iz][ix]*rhof[iz][ix];
rou12=(1-tao)*por[iz][ix]*rhof[iz][ix];
rou22=tao*por[iz][ix]*rhof[iz][ix];
D11[iz][ix]=rou11/(rou11*rou22-rou12*rou12);
D12[iz][ix]=rou12/(rou11*rou22-rou12*rou12);
D22[iz][ix]=rou22/(rou11*rou22-rou12*rou12);*/
}
else
{

vp[iz][iy][ix]=4000.0; //纵波
vs[iz][iy][ix]=2340.0; //横波
rho[iz][iy][ix]=2500.0; //密度
vf[iz][iy][ix]=1500.0;
rhof[iz][iy][ix]=1000.0;
por[iz][iy][ix]=0.2;
//b[iz][ix]=0.2*pow(10,-3);// b=eta*por*por/perm*pow(10,-3);

mu[iz][iy][ix]=rho[iz][iy][ix]*vs[iz][iy][ix]*vs[iz][iy][ix];
lamda2u[iz][iy][ix]=rho[iz][iy][ix]*vp[iz][iy][ix]*vp[iz][iy][ix];
lamda[iz][iy][ix]=lamda2u[iz][iy][ix]-2.0*mu[iz][iy][ix];
/* Ks1=rho[iz][ix]*vp[iz][ix]*vp[iz][ix];//岩石固态颗粒的体积模量
Q[iz][ix]=0.0;
// Q[iz][ix]=0.953*pow(10,9);
R[iz][ix]=Ks1;
// R[iz][ix]=0.331*pow(10,9);
tao=1.0;//孔隙弯曲度
rou11=0.0;
rou12=0.0;
rou22=rhof[iz][ix];
D11[iz][ix]=0.0;
D12[iz][ix]=0.0;
D22[iz][ix]=1/rou22;*/
}
}
}
}
}
​

//将模型扩边,用于PML
//具体的操作过程是将实际模型参数放置在扩边后的数据中央，四周的数据用
//最外缘的数据填充_3D
float ***extmodel(float ***init_model,int nz,int ny,int nx,int np)
{
float ***p;
int i,j,k;
int nx2=nx+2*np;
int ny2=ny+2*np;
int nz2=nz+2*np;


p=space3d(nz2,ny2,nx2);


for(i=0; i<np; i++)
for(k=0;k<np;k++)
for(j=0; j<np; j++)
p[i][k][j]=init_model[0][0][0];
for(i=0; i<np; i++)
for(k=np;k<np+ny;k++)
for(j=0; j<np; j++)
p[i][k][j]=init_model[0][k-np][0];
for(i=0; i<np; i++)
for(k=np+ny;k<ny2;k++)
for(j=0; j<np; j++)
p[i][k][j]=init_model[0][ny-1][0];

for (i=np; i<nz+np; i++)
for(k=0;k<np;k++)
for (j=0; j<np; j++)
p[i][k][j]=init_model[i-np][0][0];
for (i=np; i<nz+np; i++)
for(k=np;k<np+ny;k++)
for (j=0; j<np; j++)
p[i][k][j]=init_model[i-np][k-np][0];
for (i=np; i<nz+np; i++)
for(k=ny+np;k<ny2;k++)
for (j=0; j<np; j++)
p[i][k][j]=init_model[i-np][ny-1][0];

for (i=nz+np; i<nz2; i++)
for(k=0;k<np;k++)
for (j=0; j<np; j++)
p[i][k][j]=init_model[nz-1][0][0];
for (i=nz+np; i<nz2; i++)
for(k=np;k<np+ny;k++)
for (j=0; j<np; j++)
p[i][k][j]=init_model[nz-1][k-np][0];
for (i=nz+np; i<nz2; i++)
for(k=ny+np;k<ny2;k++)
for (j=0; j<np; j++)
p[i][k][j]=init_model[nz-1][ny-1][0];



for(i=0; i<np; i++)
for(k=0;k<np;k++)
for(j=np; j<np+nx; j++)
p[i][k][j]=init_model[0][0][j-np];
for(i=0; i<np; i++)
for(k=np;k<np+ny;k++)
for(j=np; j<np+nx; j++)
p[i][k][j]=init_model[0][k-np][j-np];
for(i=0; i<np; i++)
for(k=ny+np;k<ny2;k++)
for(j=np; j<np+nx; j++)
p[i][k][j]=init_model[0][ny-1][j-np];

for (i=np; i<nz+np; i++)
for(k=0;k<np;k++)
for (j=np; j<np+nx; j++)
p[i][k][j]=init_model[i-np][0][j-np];
for (i=np; i<nz+np; i++)
for(k=np;k<np+ny;k++)
for (j=np; j<np+nx; j++)
p[i][k][j]=init_model[i-np][k-np][j-np];
for (i=np; i<nz+np; i++)
for(k=ny+np;k<ny2;k++)
for (j=np; j<np+nx; j++)
p[i][k][j]=init_model[i-np][ny-1][j-np];

for (i=nz+np; i<nz2; i++)
for(k=0;k<np;k++)
for (j=np; j<np+nx; j++)
p[i][k][j]=init_model[nz-1][0][j-np];
for (i=nz+np; i<nz2; i++)
for(k=np;k<np+ny;k++)
for (j=np; j<np+nx; j++)
p[i][k][j]=init_model[nz-1][k-np][j-np];
for (i=nz+np; i<nz2; i++)
for(k=ny+np;k<ny2;k++)
for (j=np; j<np+nx; j++)
p[i][k][j]=init_model[nz-1][ny-1][j-np];


for(i=0; i<np; i++)
for(k=0;k<np;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[0][0][nx-1];
for(i=0; i<np; i++)
for(k=np;k<np+ny;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[0][k-np][nx-1];
for(i=0; i<np; i++)
for(k=ny+np;k<ny2;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[0][ny-1][nx-1];

for (i=np; i<nz+np; i++)
for(k=0;k<np;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[i-np][0][nx-1];
for (i=np; i<nz+np; i++)
for(k=np;k<np+ny;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[i-np][k-np][nx-1];
for (i=np; i<nz+np; i++)
for(k=ny+np;k<ny2;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[i-np][ny-1][nx-1];

for (i=nz+np; i<nz2; i++)
for(k=0;k<np;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[nz-1][0][nx-1];
for (i=nz+np; i<nz2; i++)
for(k=np;k<np+ny;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[nz-1][k-np][nx-1];
for (i=nz+np; i<nz2; i++)
for(k=ny+np;k<ny2;k++)
for(j=np+nx; j<nx2; j++)
p[i][k][j]=init_model[nz-1][ny-1][nx-1];


return p;
}

[qq_32787919]

程序太长了，麻烦了！编译没错，运行出错！！！！

[赵4老师]

仅供参考：

#include <stdio.h>
#include <stdlib.h>
int **newarr2d(int rows,int cols) {
    int **p,i;

    p=new int *[rows];
    if (NULL==p) exit(1);
    for (i=0;i<rows;i++) {
        p[i]=new int[cols];
        if (NULL==p[i]) exit(1);
    }
    return p;
}
void deletearr2d(int **p,int rows) {
    int i;

    for (i=0;i<rows;i++) {
        delete[] p[i];
    }
    delete[] p;
}
int main() {
    int **arr2d,i,j,r,c;

    r=4;
    c=5;
    //在堆中开辟一个4×5的二维int数组
    arr2d=newarr2d(r,c);
    for (i=0;i<r;i++) {
        for (j=0;j<c;j++) {
            arr2d[i][j]=i*c+j;
        }
    }
    for (i=0;i<r;i++) {
        for (j=0;j<c;j++) {
            printf(" %2d",arr2d[i][j]);
        }
        printf("\n");
    }
    deletearr2d(arr2d,r);

    r=6;
    c=3;
    //在堆中开辟一个6×3的二维int数组
    arr2d=newarr2d(r,c);
    for (i=0;i<r;i++) {
        for (j=0;j<c;j++) {
            arr2d[i][j]=i*c+j;
        }
    }
    for (i=0;i<r;i++) {
        for (j=0;j<c;j++) {
            printf(" %2d",arr2d[i][j]);
        }
        printf("\n");
    }
    deletearr2d(arr2d,r);

    return 0;
}
//  0  1  2  3  4
//  5  6  7  8  9
// 10 11 12 13 14
// 15 16 17 18 19
//  0  1  2
//  3  4  5
//  6  7  8
//  9 10 11
// 12 13 14
// 15 16 17
//

[赵4老师]

崩溃的时候在弹出的对话框按相应按钮进入调试，按Alt+7键查看Call Stack即“调用堆栈”里面从上到下列出的对应从里层到外层的函数调用历史。双击某一行可将光标定位到此次调用的源代码或汇编指令处，看不懂时双击下一行，直到能看懂为止。

[棉猴]

在申请三维动态数组函数space3d中修改，估计是分配内存时出错了！

float ***space3d(int nr,int ny,int nc)
{
float ***b;
int i,j;
b=(float ***)malloc(sizeof(float **) *nr);
for(i=0;i<nr;i++)
{
b[i]=(float **)malloc(sizeof(float *) *ny);
}
for(i=0;i<nr;i++)
{
for(j=0;j<ny;j++)
{
b[i][j]=(float *)malloc(sizeof(float)*nc);//sizeof(float) nc改为sizeof(float)*nc
}}

[棉猴]

修改后，VS2015运行没问题

[qq_32787919]

这个问题 改过了，我在vc++6.0 里面还是运行不了，请问：在vs里面 运行 有it=0,it=100吗？

[棉猴]

引用 6 楼 qq_32787919 的回复:

这个问题 改过了，我在vc++6.0 里面还是运行不了，请问：在vs里面 运行 有it=0,it=100吗？

逗号表达式

[qq_32787919]

it=0
it=100
it=200
这个malloc 在vc中使用，是有问题吗？

[qq_32787919]

引用 7 楼 hou09tian 的回复:

[quote=引用 6 楼 qq_32787919 的回复:]
这个问题 改过了，我在vc++6.0 里面还是运行不了，请问：在vs里面 运行 有it=0,it=100吗？

逗号表达式[/quote]
it=0
it=100
it=200
这个malloc 在vc中使用，是有问题吗？

[qq_32787919]

引用 5 楼 hou09tian 的回复:

修改后，VS2015运行没问题

这个 换成new 还是有问题？
float ***space3d(int nr,int ny,int nc)
{
float ***b;

b = new float **[nr];

for(int i=0;i<nr;i++)
{
b[i] = new float *[ny];
}
for(i=0;i<nr;i++)
{
for(int j=0;j<ny;j++)
{
b[i][j]=new float [nc];
}}
return b;
}