C程序执行有错，求教~
我写的程序编译连接没有错，在执行时出现****.exe遇到问题需要关闭，我们对此引起不便表示抱歉。请各位帮我分析一下，找找错 ~不胜感激源程序：
#include<stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.141592653589793
//#define N 6
void get_pet_matrix(double *p,int *e,int *t,int *tp,double *onxy,int np,int ne, int nt); // 读取三角元的信息
void show_p_e_t(double *p,int *e,int *t,int *tp,double *onxy,int np,int ne, int nt);// 显示三角元的信息
double CAB_lk(int Sn_t,int nt,double *p,int *e,int *t,int *tp,double *onxy,double *v_t,double detJs,int et[],double Kk[6][6],double Ki[6][6],double A_lk[6][6],double B_lk[6][6],double AB_lk[6][6],double K_lk[6][6],double dispit_s[6],double *A_lk_t,double *B_lk_t,double *AB_lk_t,double *K_lk_t);//计算相关系数
void show_A_B_lk(double A_lk[6][6],double B_lk[6][6],double AB_lk[6][6],double K_lk[6][6],double dispit_s[6]);//显示计算得到的系数
double m_inv(double AM[6][6]);//计算逆矩阵
void Ctime(int tmax,double dt,double tdelay,double gz[],double gw[],double f0,double M0,int nt,int Sn_t,int *t,int *e,double detJs,double Mkl[6][6],double Kk[6][6],double Ki[6][6],double *K_lk_t,double *A_P_t,double *B_P_t,double dispit_s[],double dispit[]);//时间循环
void main()
{
int n;
int tmax;//时间迭代次数
int Sn_t;//震源所在三角单元编号
int i ;
//int i,j,k,l;
double dt;//时间采样间隔
double f0;//Ricker子波主频
double tdelay;//Ricker子波时间延迟
double M0;
double Ksai,Ita;//标准坐标系下坐标
int np;//节点总数
int ne;//边总数
int nt;//三角单元总数
double *p;//节点矩阵P，大小为nt*3，记录所有节点的编号及实际的xy坐标值
int *e; //边矩阵e，大小为nt*5，记录每条边的编号，边的两个节点，边所在的三角单元 （两个三角单元共用一条边）
int *t;//单元矩阵t，大小为nt*5，记录所有单元的编号，组成单元的三条边，单元所属的计算区域
int *tp;//矩阵tp，大小为nt*5，记录所有单元的编号，组成单元的三个节点，单元所属的计算区域
double *onxy; //矩阵onxy，大小为nt*3*2，记录组成单元的三条边的法向单位矢量
double *v_t; //所有单元的介质速度参数
double *A_lk_t,*B_lk_t,*AB_lk_t,*K_lk_t;
//double *A_P,*B_P,*C_P;
double *A_P_t,*B_P_t,*C_P_t;
double detJs=0.0;//所有单元的雅克比矩阵值
double A_lk[6][6],B_lk[6][6],AB_lk[6][6],K_lk[6][6];
double dispit[6];//基函数
//double PHIs[6]={1.0, -1.0/4, -0.1875, -1.0/4, -0.0625, -0.375};
//矩阵Mkl
double Mkl[6][6]={
{1.0/2, 0.0, 0.0, 0.0, 0.0, 0.0},
{0.0, 1.0/12, 0.0, 0.0, 0.0, 0.0},
{0.0, 0.0, 1.0/30, 0.0, 0.0, 0.0},
{0.0, 0.0, 0.0, 1.0/4, 0.0, 0.0},
{0.0, 0.0, 0.0, 0.0, 1.0/18, 0.0},
{0.0, 0.0, 0.0, 0.0, 0.0, 1.0/6}};
//矩阵Kk
double Kk[6][6]={ // phi_l/ksai*phi_k/ksai
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
{0.0, 2.0, 0.0, 0.0, 4.0/3, 0.0},
{0.0, 0.0, 3.0, 0.0, 0.0, 0.0},
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
{0.0, 4.0/3.0, 0.0, 0.0, 11.0/3.0, 0.0},
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
double Ki[6][6]={
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // phi_l/yita*phi_k/yita
{0.0, 1.0/2, 1.0/3, 3.0/2.0, 1.0/3, -2.0/3},
{0.0, 1.0/3, 1.0, 1.0, -2.0, 7.0/3},
{0.0, 3.0/2.0, 1.0, 9.0/2, 1.0, -2.0},
{0.0, 1.0/3, -2.0, 1.0, 3.0, 7.0/3},
{0.0, -2.0/3, 7.0/3, -2.0, 7.0/3, 12.0}};
int gn=10; //Lengdre多项式求解震源系数中高斯数值积分节点数
double gz[10]={-0.9739065285,-0.8650633677,-0.6794095683,-0.4333953941,-0.1488743990,0.1488743990,0.4333953941,
0.6794095683,0.8650633677,0.9739065285}; //高斯数值积分节点
double gw[10]={0.0666713443,0.1494513492,0.2190863625,0.2692667193,0.2955242247,0.2955242247,0.2692667193,
0.2190863625,0.1494513492,0.0666713443};//高斯数值积分系数
double dispit_s[6]={0.0};
int et[3]={0};
n=2;
dt=0.0001;
tmax=1000;
f0=30.0;
tdelay=0.050;
M0=1000000.0;//震源幅度
Sn_t=5172;//震源所在三角单元
Ksai=1.0/4; //标准单元下显示点的横坐标
Ita=1.0/4; //标准单元下显示点的纵坐标
np=2916;//节点总数
ne=8545;//边总数
nt=5630;//单元总数
p=(double *)malloc(sizeof(double)*np*3);
e=(int *)malloc(sizeof(int)*ne*5);
t=(int *)malloc(sizeof(int)*nt*5);
tp=(int *)malloc(sizeof(int)*nt*5);
onxy=(double *)malloc(sizeof(double)*nt*3*2);
v_t=(double *)malloc(sizeof(double)*nt);
A_lk_t=(double *)malloc(sizeof(double)*nt*6*6);
B_lk_t=(double *)malloc(sizeof(double)*nt*6*6);
AB_lk_t=(double *)malloc(sizeof(double)*nt*6*6);
K_lk_t=(double *)malloc(sizeof(double)*nt*6*6);
A_P_t=(double *)malloc(sizeof(double)*nt*6*6);
B_P_t=(double *)malloc(sizeof(double)*nt*6*6);
C_P_t=(double *)malloc(sizeof(double)*nt*6*6);
// 介质速度参数设置
for (i=0;i<nt;i++)
{
if (t[5*i+4]==48)
{
v_t[i]=4000.0;
}
else
{
v_t[i]=4000.0;
}
}
//基函数表达式
dispit[0]=1.0;
dispit[1]=Ita-1.0+2.0*Ksai;
dispit[2]=Ita*Ita-2.0*Ita+6.0*Ksai*Ita+6.0*Ksai*Ksai-6.0*Ksai+1.0;
dispit[3]=-1.0+3.0*Ita;
dispit[4]=5.0*Ita*Ita+10.0*Ksai*Ita-6.0*Ita-2.0*Ksai+1.0;
dispit[5]=1.0-8.0*Ita+10.0*Ita*Ita;
//读取矩阵p,e,tp等
get_pet_matrix(p,e,t,tp,onxy,np,ne,nt);
//显示矩阵 p,e,tp等
show_p_e_t(p,e,t,tp,onxy,np,ne,nt);
CAB_lk( Sn_t, nt, p, e, t, tp, onxy, v_t, detJs, et, Kk, Ki, A_lk, B_lk, AB_lk, K_lk, dispit_s, A_lk_t, B_lk_t,AB_lk_t, K_lk_t);
show_A_B_lk( A_lk, B_lk, AB_lk, K_lk, dispit_s);
m_inv(Mkl);
//时间迭代更新
Ctime( tmax, dt, tdelay, gz, gw, f0, M0, nt, Sn_t, t, e, detJs, Mkl, Kk, Ki,K_lk_t,A_P_t,B_P_t, dispit_s, dispit);
free(p);
free(e);
free(t);
free(tp);
free(onxy);
free(v_t);
free(A_lk_t);
free(B_lk_t);
free(AB_lk_t);
free(K_lk_t);
free(A_P_t);
free(B_P_t);
free(C_P_t);
}
//读取与单元有关的矩阵p，e，t等矩阵
void get_pet_matrix(double *p,int *e,int *t,int *tp,double *onxy,int np,int ne, int nt)
{
FILE *fp;
if((fp=fopen("p","rb"))==NULL)
{
printf("cannot open file\n");
return;
}
fread(p,sizeof(double),np*3,fp);
fclose(fp);
if((fp=fopen("e","rb"))==NULL)
{
printf("cannot open file\n");
return;
}
fread(e,sizeof(int),ne*5,fp);
fclose(fp);
if((fp=fopen("t","rb"))==NULL)
{
printf("cannot open file\n");
return;
}
fread(t,sizeof(int),nt*5,fp);
fclose(fp);
if((fp=fopen("tp","rb"))==NULL)
{
printf("cannot open file\n");
return;
}
fread(tp,sizeof(int),nt*5,fp);
fclose(fp);
if((fp=fopen("onxy","rb"))==NULL)
{
printf("cannot open file\n");
return;
}
fread(onxy,sizeof(double),nt*3*2,fp);
fclose(fp);
/*if((fp=fopen("v_t","rb"))==NULL)
{
printf("cannot open file\n");
return;
}
fread(v_t,sizeof(double),nt*2,fp);
fclose(fp); */
}
//显示p，e，t等与网格信息有关的矩阵
void show_p_e_t(double *p,int *e,int *t,int *tp,double *onxy,int np,int ne, int nt)
{
int i,j;
FILE *fp0,*fp1,*fp2,*fp3,*fp4,*fp5;
fp0=fopen("p.txt","w");
for (i=0;i<np;i++)
{
for (j=0;j<3;j++)
fprintf(fp0,"%f ",p[i*3+j]);
fprintf(fp0,"\n");
}
fclose(fp0);
fp1=fopen("e.txt","w");
for (i=0;i<ne;i++)
{
for (j=0;j<5;j++)
fprintf(fp1,"%d ",e[i*5+j]);
fprintf(fp1,"\n");
}
fclose(fp1);
fp2=fopen("t.txt","w");
for (i=0;i<nt;i++)
{
for (j=0;j<5;j++)
fprintf(fp2,"%d ",t[i*5+j]);
fprintf(fp2,"\n");
}
fclose(fp2);
fp3=fopen("tp.txt","w");
for (i=0;i<nt;i++)
{
for (j=0;j<5;j++)
fprintf(fp3,"%d ",tp[i*5+j]);
fprintf(fp3,"\n");
}
fclose(fp3);
fp4=fopen("onxy.txt","w");
for (i=0;i<3*nt;i++)
{
for (j=0;j<2;j++)
fprintf(fp4,"%f ",onxy[i*2+j]);
fprintf(fp4,"\n");
}
fclose(fp4);
/*fp5=fopen("v_t.txt","w");
for (i=0;i<nt;i++)
{
for (j=0;j<2;j++)
fprintf(fp5,"%f ",v_t[i+j]);
fprintf(fp5,"\n");
}
fclose(fp5);*/
}
// 计算A_lk=Kk_lk*(dksi/dx)^2 and B_lk=Ki_lk*(dyita/dy)^2
double CAB_lk(int Sn_t,int nt,double *p,int *e,int *t,int *tp,double *onxy,double *v_t,double detJs,int et[],double Kk[6][6],double Ki[6][6],double A_lk[6][6],double B_lk[6][6],double AB_lk[6][6],double K_lk[6][6],double dispit_s[6],double *A_lk_t,double *B_lk_t,double *AB_lk_t,double *K_lk_t)
{
int p1,p2,p3;//所计算单元的三个节点
double x1,y1,x2,y2,x3,y3;//所有节点的实际坐标
//double dispit_s[6];//基函数在震源点处的值
double v,x_s,y_s;
double Ksai_s,Ita_s;
// double A_lk[6][6]={{0.0}};
//double B_lk[6][6]={{0.0}}; // 所有单元的有限元方程组的系数
//double K_lk[6][6]={{0.0}};
//double AB_lk[6][6]={{0.0}};
int i,j;
int in;
// double v;
in=Sn_t-1;
//读取组成单元的三条边
et[0]=t[(in)*5+1];
et[1]=t[(in)*5+2];
et[2]=t[(in)*5+3];
//读取组成单元的三个节点
p1=tp[in*5+1];
p2=tp[in*5+2];
p3=tp[in*5+3];
//读取所有节点的实际坐标值
x1=p[(p1-1)*3+1];
x2=p[(p2-1)*3+1];
x3=p[(p3-1)*3+1];
y1=p[(p1-1)*3+2];
y2=p[(p2-1)*3+2];
y3=p[(p3-1)*3+2];
detJs=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1);//该计算单元的雅克比行列式值
x_s=250.0; // 震源位置
y_s=350.0;
Ksai_s=((x3*y1-x1*y3)+x_s*(y3-y1)+y_s*(x1-x3))/detJs;
Ita_s=((x1*y2-x2*y1)+x_s*(y1-y2)+y_s*(x2-x1))/detJs; //震源位置在标准单元坐标系中的位置
// 基函数在震源位置处的值
dispit_s[0]=1.0;
dispit_s[1]=Ita_s-1.0+2.0*Ksai_s;
dispit_s[2]=Ita_s*Ita_s-2.0*Ita_s+6.0*Ksai_s*Ita_s+6.0*Ksai_s*Ksai_s-6.0*Ksai_s+1.0;
dispit_s[3]=-1.0+3.0*Ita_s;
dispit_s[4]=5.0*Ita_s*Ita_s+10.0*Ksai_s*Ita_s-6.0*Ita_s-2.0*Ksai_s+1.0;
dispit_s[5]=1.0-8.0*Ita_s+10.0*Ita_s*Ita_s;
v=v_t[in]; //存储计算单元的速度值
for (i=0;i<6;i++)
for (j=0;j<6;j++)
{
A_lk[i][j]=Kk[i][j]*(y3-y1)*(y3-y1)/detJs;
B_lk[i][j]=Ki[i][j]*(x2-x1)*(x2-x1)/detJs;
AB_lk[i][j]=A_lk[i][j]+B_lk[i][j];
K_lk[i][j]=v*v*AB_lk[i][j];
}
for (in=1;in<nt+1;in++)
for (i=0;i<6;i++)
for (j=0;j<6;j++)
{
A_lk_t[(in-1)*6*6+i*6+j]=A_lk[i][j];
B_lk_t[(in-1)*6*6+i*6+j]=B_lk[i][j];
AB_lk_t[(in-1)*6*6+i*6+j]=AB_lk[i][j];
K_lk_t[(in-1)*6*6+i*6+j]=K_lk[i][j];
}
return 0.0;
}
void show_A_B_lk(double A_lk[6][6],double B_lk[6][6],double AB_lk[6][6],double K_lk[6][6],double dispit_s[6])
{
int i,j;
FILE *fp0,*fp1,*fp2,*fp3,*fp4;
fp0=fopen("A_lk.txt","w");
for (i=0;i<6;i++)
{
for(j=0;j<6;j++)
fprintf(fp0,"%f",A_lk[i][j]);
fprintf(fp0,"\n");
}
fclose(fp0);
fp1=fopen("B_lk.txt","w");
for (i=0;i<6;i++)
{
for(j=0;j<6;j++)
fprintf(fp1,"%f",B_lk[i][j]);
fprintf(fp1,"\n");
}
fclose(fp1);
fp2=fopen("AB_lk.txt","w");
for (i=0;i<6;i++)
{
for(j=0;j<6;j++)
fprintf(fp2,"%f",AB_lk[i][j]);
fprintf(fp2,"\n");
}
fclose(fp2);
fp3=fopen("K_lk.txt","w");
for (i=0;i<6;i++)
{
for(j=0;j<6;j++)
fprintf(fp3,"%f",K_lk[i][j]);
fprintf(fp3,"\n");
}
fclose(fp3);
fp4=fopen("A_lk.txt","w");
for (i=0;i<6;i++)
{
fprintf(fp4,"%f",dispit_s[i]);
fprintf(fp4,"\n");
}
fclose(fp4);
}
//矩阵求逆
double m_inv(double AM[6][6])
{
int nM=14;
int iM,jM,kM;
double dM;
int JSM[14];//error
int ISM[14];//error
for (kM=0;kM<nM;kM++)
{
dM=0;
for (iM=kM;iM<nM;iM++)
for (jM=kM;jM<nM;jM++)
{
if (fabs(AM[iM][jM])>dM)
{
dM=fabs(AM[iM][jM]);
ISM[kM]=iM;
JSM[kM]=jM;
}
}
if (dM+1.0==1.0)
return 0;
if (ISM[kM]!=kM)
for (jM=0;jM<nM;jM++)
{
double c=AM[kM][jM];
AM[kM][jM]=AM[ISM[kM]][jM];
AM[ISM[kM]][jM]=c;
}
if (JSM[kM]!=kM)
for (iM=0;iM<nM;iM++)
{
double c=AM[iM][kM];
AM[iM][kM]=AM[iM][JSM[kM]];
AM[iM][JSM[kM]]=c;
}
AM[kM][kM]=1/AM[kM][kM];
for (jM=0;jM<nM;jM++)
if (jM!=kM)
AM[kM][jM]=AM[kM][jM]*AM[kM][kM];
for (iM=0;iM<nM;iM++)
if (iM!=kM)
for (jM=0;jM<nM;jM++)
if (jM!=kM)
AM[iM][jM]=AM[iM][jM]-AM[iM][kM]*AM[kM][jM];
for (iM=0;iM<nM;iM++)
if (iM!=kM)
AM[iM][kM]=-AM[iM][kM]*AM[kM][kM];
}
for (kM=nM-1;kM>=0;kM--)
{
for (jM=0;jM<nM;jM++)
if (JSM[kM]!=kM)
{
double c=AM[kM][jM];
AM[kM][jM]=AM[JSM[kM]][jM];
AM[JSM[kM]][jM]=c;
}
for (iM=0;iM<nM;iM++)
if (ISM[kM]!=kM)
{
double c=AM[iM][kM];
AM[iM][kM]=AM[iM][ISM[kM]];
AM[iM][ISM[kM]]=c;
}
}
return 0.0;
}
//时间层迭代更新
void Ctime(int tmax,double dt,double tdelay,double gz[],double gw[],double f0,double M0,int nt,int Sn_t,int *t,int *e,double detJs,double Mkl[6][6],double Kk[6][6],double Ki[6][6],double *K_lk_t,double *A_P_t,double *B_P_t, double dispit_s[],double dispit[])
{
int it,inp;//时间迭代变量
int i,j;
int in;
int k,l,o;
//double *u,*ua,*up,*ul;
double u[6]={0.0};
double uap[6]={0.0};
double ua[6]={0.0};
double *dwti1,*ul;
double *uaip,*dwtp1,*up,*uip;
double psai[3];
int gn=10;
ul=(double *)malloc(sizeof(double)*nt*6*tmax);//所有时刻所有单元里u的系数
dwtp1=(double *)malloc(sizeof(double)*nt); // 存储所有单元里u的近似值
//所有进程需要存储的的系数
up=(double *)malloc(sizeof(double)*nt*6); //存储所有单元ul(t+dt)时刻的值
uaip=(double *)malloc(sizeof(double)*nt*6);//存储所有单元ul(t)时刻的值
uip=(double *)malloc(sizeof(double)*nt*6);//存储所有单元ul(t-dt)时刻的值
//up=(double *)malloc(sizeof(double)*nt*6);
//dwt1=(double *)malloc(sizeof(double)*nt);
dwti1=(double *)malloc(sizeof(double)*nt*tmax); // 存储所有时刻所有单元里的u的近似值
//初始化u，ua
for (i=0;i<nt*6;i++)
{
u[i]=0.0;
ua[i]=0.0;
uap[i]=0.0;
}
printf("the time is starting!...\n");
for (it=0;it<tmax;it++)
{
double S=0.0;
double So[3]={0.0};
double gg[3]={0.0};
double ga,gb;
double gy,gt,gaa,gricker;
double S_t[6]={0.0};
//计算震源的系数Spo
S=1.0;
ga=it*dt-tdelay;
gb=ga+dt;
if(it%20==0)
printf("the time is %d\n",it);
for (i=0;i<gn;i++)
{
gy=(gz[i]*(gb-ga)+ga+gb)/2.0;//将【ga-gb】间隔转化为Legendre多项式的定义域【-1，1】
gt=gy;
gaa=PI*f0*gt;
gricker=-M0*(1-2*gaa*gaa)*exp(-gaa*gaa);
gg[0]=gg[0]+gw[i]*gricker;
gg[1]=gg[1]+gw[i]*gricker*gz[i];
gg[2]=gg[2]+gw[i]*gricker*(3.0*gz[i]*gz[i]-1.0)/2.0;
}
gg[0]=gg[0]*(gb-ga)/2;
gg[1]=gg[1]*(gb-ga)/2; //求解数值积分
gg[2]=gg[2]*(gb-ga)/2;
// Legendre 多项式
psai[0]=1;
psai[1]=ga;
psai[2]=(3.0*ga*ga-1.0)/2.0;
So[0]=S*gg[0]/2.0/(dt/2.0);
So[1]=S*gg[1]/(1/3.0)/(dt/2.0);
So[2]=S*gg[2]/(2.0/5.0)/(dt/2.0);//求解震源时间上用Legendre多项式分解的分解系数
for(k=0;k<6;k++)
{
for(o=0;o<3;o++)
{
S_t[k]=S_t[k]+So[o]*psai[o];
}
S_t[k]=S_t[k]*dispit_s[k];
}
for(k=0;k<6;k++)
{
S_t[k]=S_t[k]*detJs;
}
//
for(inp=0;inp<nt;inp++) // 所有三角单元里计算
{
in=inp;
if (in>nt)
{
in=nt-1;
}
m_inv(Mkl);
// 计算有限元离散方程组的系数
for(i=0;i<6;i++)
for(j=0;j<6;j++)
{
A_P_t[(in-1)*6*6+i*6+j]=Mkl[i][j]/dt/dt;
B_P_t[(in-1)*6*6+i*6+j]=-A_P_t[(in-1)*6*6+i*6+j]*K_lk_t[(in-1)*6*6+i*6+j];
//C_P_t[(in-1)*6*6+i*6+j]=-A_P_t[(in-1)*6*6+i*6+j];
}
for (k=0;k<6;k++)
{
for (l=0;l<6;l++)
{
dt=pow(dt,2);
//A_P_t[(in-1)*6*6+k*6+l]*uap[l]=B_P_t[(in-1)*6*6+k*6+l]*u[l]+C_P_t[(in-1)*6*6+k*6+l]*ua[l]+S_t[k]; //
uap[k]=uap[k]+2.0/dt*u[l]+B_P_t[(in-1)*6*6+k*6+l]*u[l]-ua[l]+A_P_t[(in-1)*6*6+k*6+l]*S_t[k];
}
}
dwtp1[inp]=0.0;
for (i=0;i<6;i++)
{ up[inp*6+i]=uap[i]; //存储所有计算单元(t+dt)时刻的离散方程组的系数值
uip[inp*6+i]=u[i]; //存储所有计算单元(t)时刻的离散方程组的系数值
uaip[inp*6+i]=ua[i]; //存储所有计算单元(t-dt)时刻的离散方程组的系数值
dwtp1[inp]=dwtp1[inp]+up[inp*6+i]*dispit[i]; //计算所有单元里u的近似值
}
}//in end
//% update
for (i=0;i<nt;i++)//
{
for (k=0;k<6;k++)
{
/*u[i*6+k]=uap[i*6+k];
ua[i*6+k]=u[i*6+k];
*/
uip[i*6+k]=up[i*6+k];
uaip[i*6+k]=uip[i*6+k];
}
ul[it*6*tmax+i*6+k]=up[i*6+k];//存储所有时刻所有单元里的离散方程组的系数
}
for (i=0;i<nt;i++)
{
dwti1[i+it*nt]=dwtp1[i]; //存储所有时间所有单元里u的近似值
}
}//it end
//存储计算结果
int count2,count1;
FILE *fp1,*fp2;
if((fp1=fopen("dwtc1_FE_npml.data","wb"))==NULL)
{
printf("cannot open file1 \n");
//return 0;
}
count1=nt*tmax;
fwrite(dwti1,sizeof(double),count1,fp1);
//if(!=count2)printf("file write error");
printf("file dwtc1 write %d (double) sucessfully!\n",count1);
fclose(fp1);
if((fp2=fopen("ul_FE_npml.data","wb"))==NULL)
{
printf("cannot open file upl\n");
//return 0;
}
count2=nt*tmax*6;
fwrite(ul,sizeof(double),count2,fp2);
//if(!=count2)printf("file write error");
printf("file ul write %d (double) sucessfully!\n",count2);
fclose(fp2);
printf("Mission Complete!\n");
//system("PAUSE");
free(dwtp1);
free(dwti1);
free(ul);
free(uip);
free(uaip);
free(up);
}
