为何文件生成了两个cell和一个double数组
本程序完全可以运行，只是奇怪为什么在calchi中的fk3文件会生成两个cell和一个double数组，其目的只是和其它文件一样只生成一个数组，其它文件输出正常。麻烦各位了。
程序代码：#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* FMG Solver of the EHL circular contact */
#define pi 3.1415926535897931
typedef struct
{
double hx;
int ii;
double *p, *f;
double *hrhs;
double *hfi;
double *w;
double *K,*K0;
double *pjac; /*old pressure for use in Jacobi relaxation*/
double rg;
double Hm;
} Level;
typedef struct
{
int nx0;
int maxlevel;
double xa,xb;
double h0; /*global constant and work unit*/
Level *Lk;
} Stack;
/********** GLOBAL VARIABLES **********************/
double MMoes, LMoes, H_0, rlambda, alphabar;
double p0r, alpha, eta0, zr, maxll;
double hfact,xi_l,urja,urgs;
void initialize(Stack *U, int nx0, int maxl, double xa, double xb, double h0)
{
double hx;
Level *L;
int i,ii;
U->xa=xa;
U->xb=xb;
U->maxlevel=maxl;
U->h0=h0;
U->Lk=(Level *)calloc(maxl+1,sizeof(Level));
hx=(xb-xa)/nx0;
ii=nx0;
for (i=1; i<=maxl; i++)
{
L=U->Lk+i;
L->hx=hx;
L->ii=ii;
L->p =(double *)calloc(L->ii+1,sizeof(double));
L->w =(double *)calloc(L->ii+1,sizeof(double));
L->f =(double *)calloc(L->ii+1,sizeof(double));
L->hrhs=(double *)calloc(L->ii+1,sizeof(double));
L->K =(double *)calloc(L->ii+1,sizeof(double));
L->K0 =(double *)calloc(L->ii+1,sizeof(double));
L->pjac =(double *)calloc(L->ii+1,sizeof(double));
L->hfi =(double *)calloc(L->ii+1,sizeof(double));
printf("\n level: %2d ii=%4d, hx=%f",i,ii,hx);
hx*=0.5; ii*=2;
}
}
/********* SINGLE GRID ROUTINES *****/
void init_f(Stack *U, int lev)
{
int i;
Level *L;
double x;
L=U->Lk+lev;
for (i=0; i<=L->ii; i++)
{
x=U->xa+i*L->hx;
L->f[i]=0.0;
L->hrhs[i]=0.5*x*x;
}
FILE *ff1,*fhrhs1;
ff1=fopen("f1.dat","w");
fhrhs1=fopen("hrhs1.dat","w");
for (i=0; i<=L->ii; i++)
{
x=U->xa+i*L->hx;
fprintf(ff1,"%f\n",L->f[i]);
fprintf(fhrhs1,"%f\n",L->hrhs[i]);
}
fprintf (ff1, "\n");
fprintf (fhrhs1, "\n");
fclose(ff1);
fclose(fhrhs1);
L->rg=-pi/2.0;
}
void init_log(Stack *U, int lev)
{
/* computes the kernel on level l */
int i;
Level *L;
double xp,xm;
L=U->Lk+lev;
for (i=0; i<=L->ii;i++)
{
xp=(i+0.5)*L->hx; xm=xp-L->hx;
L->K[i]=xp*(log(xp)-1)-xm*(log(xm)-1);
//L->K0[i]=(i+1/2)*(L->hx)*((log(fabs(i+1/2)*L->hx))-1)-(i-1/2)*(L->hx)*((log(fabs(i-1/2)*L->hx))-1);
}
}
void init_p(Stack *U, int lev)
{
int i;
Level *L;
double x;
void calchi(Stack *U, int lev);
L =U->Lk+lev;
for (i=0; i<=L->ii; i++)
{
x=U->xa+i*L->hx;
L->p[i]=0.0;
if (x*x<1.0) L->p[i]=sqrt(1.0-x*x);
}
FILE *fp1;
fp1=fopen("fp1.dat","w");
for (i=0; i<=L->ii; i++)
{
x=U->xa+i*L->hx;
fprintf(fp1,"%f\n",L->p[i]);
}
fprintf (fp1, "\n");
fclose(fp1);
calchi(U,lev);
}
/********** SPECIAL FUNCTIONS **********/
double reta(double p)
{
/* Barus */
return (exp(-alphabar*p));
/* Roelands */
/* return(exp(-alpha*p0r/zr*(-1.0+pow(1.0+(p/p0r)*(alphabar/alpha),zr)))) ; */
}
double rho(double p)
{
/* Incompressible */
return (1.0);
/* Compressible */
/* return( (5.9e8+1.34*(alphabar/alpha)*p)/(5.9e8+(alphabar/alpha)*p)); */
}
/********relax****************/
void relax(Stack *U, int lev)
/* relaxes the equation on level l */
{
int i;
Level *L;
double *P, *Pj, *Hfi, *KK;
double g,hx;
double r0,r0w,r0e,kesai_w,kesai_e,ri,pianli,pianlii;
void calchi(Stack *U, int l);
L =U->Lk+lev;
P =L->p;
Pj=L->pjac;
Hfi=L->hfi;
hx=L->hx;
KK=L->K;
calchi(U,lev); /*generate L->hfi*/
for (i=1; i<=L->ii; i++)
{
r0=rho(P[i]);
r0w=rho(P[i-1]);
r0e=rho(P[i+1]);
kesai_w=0.5*(r0*Hfi[i]*Hfi[i]*Hfi[i]*reta(P[i])*rlambda+r0w*Hfi[i-1]*Hfi[i-1]*Hfi[i-1]*reta(P[i-1])*rlambda);
kesai_e=0.5*(r0*Hfi[i]*Hfi[i]*Hfi[i]*reta(P[i])*rlambda+r0e*Hfi[i+1]*Hfi[i+1]*Hfi[i+1]*reta(P[i+1])*rlambda);
if ((kesai_w/hx/hx>1)||(kesai_e/hx/hx>1))
{
ri=-(kesai_w*P[i-1]-(kesai_w+kesai_e)*P[i]+kesai_e*P[i+1])/hx/hx+(r0*Hfi[i]-r0w*Hfi[i-1])/hx;
pianli=-(kesai_w+kesai_e)/hx/hx+(r0*L->K[0]-r0w*L->K[1])/pi/hx;
P[i]=P[i]+0.3*ri/pianli;
}
else
{
ri=-(kesai_w*Pj[i-1]-(kesai_w+kesai_e)*Pj[i]+kesai_e*Pj[i+1])/hx/hx+(r0*Hfi[i]-r0w*Hfi[i-1])/hx;
pianlii=-(2*kesai_w+kesai_e)/hx/hx+2*(r0*KK[0]-r0w*KK[1])/pi/hx;
P[i]=P[i]+0.2*ri/pianlii;
P[i-1]=P[i-1]-0.2*ri/pianlii;
if (P[i-1]<0) P[i-1]=0;
}
if (P[i]<0) P[i]=0;
}
g=0.0;
for (i=0; i<=L->ii; i++) g+=L->p[i];
g=g*L->hx+L->rg;
}
void relaxh0(Stack *U, int lev)
/* relaxes the force-balance equation */
/* hfact <0.05 gives stable convergence for W cycles */
{
int i;
Level *L;
double g;
L=U->Lk+lev;
g=0.0;
for (i=0; i<=L->ii; i++) g+=L->p[i];
g=g*L->hx+L->rg;
U->h0+=0.05*g;
}
void calchi(Stack *U, int lev)
{
int i,i1;
Level *L;
L=U->Lk+lev;
double help,x;
for (i1=0; i1<=L->ii; i1++)
{
help=0.0;
for (i=0; i<=L->ii; i++) help+=L->K[abs(i-i1)]*L->p[i];
L->w[i1]=help;
}
for (i=0; i<=L->ii; i++) L->hfi[i]=U->h0+2.0/pi/pi*L->w[i]+L->hrhs[i];
FILE *fk3,*fhrhs3,*fp3;
fk3=fopen("fk3.dat","w");
fhrhs3=fopen("fhrhs3.dat","w");
fp3=fopen("fp3.dat","w");
for (i=0; i<=L->ii; i++)
{
x=U->xa+i*L->hx;
fprintf(fk3, "%f\n",L->K[i]);
fprintf(fhrhs3,"%f\n",L->hrhs[i]);
fprintf(fp3,"%f\n",L->p[i]);
}
fprintf(fk3, "\n");
fprintf(fhrhs3, "\n");
fprintf(fp3, "\n");
fclose(fk3);
fclose(fhrhs3);
fclose(fp3);
}
/********** output **********/
void output(Stack *U)
/* writes an output file of p and h */
/* output of Hm and Hc to screen */
{
int i;
Level *L;
double x;
FILE *fp,*fh;
L=U->Lk+U->maxlevel;
fp=fopen("Pp.dat","w");
fh=fopen("Hh.dat","w");
//L->Hm=1e5; /* arbitrary large value */
for (i=0; i<=L->ii; i++)
{
x=U->xa+i*L->hx;
fprintf(fp,"%f %f\n",x,L->p[i]);
fprintf(fh,"%f %f\n",x,L->hfi[i]);
}
fprintf (fp, "\n");
fprintf (fh, "\n");
fclose(fp);
fclose(fh);
}
/********** MAIN PROGRAM **********/
void main()
{
Stack U;
printf("\ngive M ?"); scanf("%lf",&MMoes);
printf("\ngive L ?"); scanf("%lf",&LMoes);
/* conversion to parameters appearing in equations */
rlambda=1.0/(3*pi*pi/8/(MMoes*MMoes));
alphabar=LMoes*sqrt(MMoes/2/pi);
printf("\nrlambda=%8.5e alpha*p_h=%8.5e\n",rlambda,alphabar);
initialize(&U,128,1,-4.5,1.5,-0.5);//initialize(Stack *U, int nx0, int maxl, double xa, double xb, double h0)
init_log(&U,1);
init_f(&U,1);
init_p(&U,1);
//output(&U);
}
