| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 311 人关注过本帖
标题:为何文件生成了两个cell和一个double数组
只看楼主 加入收藏
huaijuliu
Rank: 1
等 级:新手上路
帖 子:20
专家分:0
注 册:2009-9-17
收藏
 问题点数:0 回复次数:2 
为何文件生成了两个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);
}

搜索更多相关主题的帖子: double 
2011-06-20 05:47
huaijuliu
Rank: 1
等 级:新手上路
帖 子:20
专家分:0
注 册:2009-9-17
收藏
得分:0 
而且这个fk3所生成的数组为128行,而其它正常文件为129行,感谢大家。
2011-06-20 05:50
huaijuliu
Rank: 1
等 级:新手上路
帖 子:20
专家分:0
注 册:2009-9-17
收藏
得分:0 
这个问题已经解决 因为xm越界了 添加fabs(xm)即可 感谢大家关注
2011-06-20 06:53
快速回复:为何文件生成了两个cell和一个double数组
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.020371 second(s), 7 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved