为何文件生成了两个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); }