📄 f4.c
字号:
#include <stdio.h>
#include <math.h>
#include "com.h"
NODE Gnode[MMAXNOD];
ELEMENT Element[MMAXELE];
long Egauss[MMAXELE][MO];
int Bd[MMAXBD];
char Btyp[MMAXBD];
int MAXNOD,MAXELE,MAXBD,MAXVERT,coor;
char Eq,D_out[30],D_in[30],Mesh[30],JOB[30],C_sol[30];
FILE *Fprss,*Fg;
double Maxerr,FE_ertol,Re,Re_max,Re_step,Re0,supg;
int N_FE;
/**********************************************************/
extern int Gnum();
extern ROW *newrow();
extern void solv(),post(),Gauss_write();
extern double *G(),*W();
/*********************************************************/
void read_mesh()
{
FILE *filep=NULL;
FILE *fileout=NULL;
int num=0;
int EleNum=0;
int Point=0;
filep=fopen("meshfile.dat","r");
fileout=fopen("checkfile","w+");
if (fileout==0||filep==0)
{
printf("This file was not open");
}
else
{
fscanf(filep,"%d\n",&MAXNOD);
fprintf(fileout,"%d\n",MAXNOD);
fscanf(filep,"%d\n",&MAXELE);
fprintf(fileout,"%d\n",MAXELE);
fscanf(filep,"%d\n",&MAXBD);
fprintf(fileout,"%d\n",MAXBD);
fscanf(filep,"%d\n",&MAXVERT);
fprintf(fileout,"%d\n",MAXVERT);
for(num=0;num<MAXNOD;num++)
{
fscanf(filep,"%lf\t",&(Gnode[num].x));
fscanf(filep,"%lf\n",&(Gnode[num].y));
fprintf(fileout,"%lf\t",Gnode[num].x);
fprintf(fileout,"%lf\n",Gnode[num].y);
}
for(EleNum=0;EleNum<MAXELE;EleNum++)
{
for(Point=0;Point<8;Point++)
{
fscanf(filep,"%d\t",&(Element[EleNum].n[Point]));
fprintf(fileout,"%d\t",Element[EleNum].n[Point]);
}
fscanf(filep,"%d\n",&(Element[EleNum].n[8]));
fprintf(fileout,"%d\n",Element[EleNum].n[8]);
}
for(num=0;num<MAXBD;num++)
{
fscanf(filep,"%d\t",Bd+num);
fscanf(filep,"%c\n",Btyp+num);
fprintf(fileout,"%d\t",Bd[num]);
fprintf(fileout,"%c\n",Btyp[num]);
}
fclose(filep);
fclose(fileout);
}
/********************************************
Open an mesh data file;
Read Gnode[n].x and Gnode[n].y for n=0,1,……MAXNOD-1 (Global)
Read Bd[i] and Btyp[i] for i=0,1,……MAXBD-1; (Global)
Read Element[en].n[i] for en=0,1…...MAXELE-1, for i=0,1,…8 (Local)
*********************************************/
}
/*********************************************************/
static void initial()
{
ELEMENT e;
int uo,vo,po,count[9];
int n,i,j,o,f;
ROW *pr; NODE *pn;
for(i=0;i<MAXELE;i++)
{
Element[i].vo=2; Element[i].po=1; Element[i].fo=2;
}
for(n=0;n<MAXNOD;n++) for(i=0;i<MB+1;i++) Gnode[n].fc[i]=0;
for(j=0;j<MAXELE;j++)
{
e=Element[j];
count[0]=count[1]=count[2]=count[3]=1;
count[4]=count[5]=count[6]=count[7]=e.vo-1; count[8]=(e.vo-1)*(e.vo-1);
for(o=0;o<2;o++) for(i=0;i<9;i++)
if(Gnode[e.n[i]].fc[o]<count[i]) Gnode[e.n[i]].fc[o]=count[i];
count[4]=count[5]=count[6]=count[7]=e.po-1; count[8]=(e.po-1)*(e.po-1);
o=2; for(i=0;i<9;i++)
if(Gnode[e.n[i]].fc[o]<count[i]) Gnode[e.n[i]].fc[o]=count[i];
o=MB; for(i=0;i<9;i++) Gnode[e.n[i]].fc[o]=1;
}
for(n=0;n<MAXNOD;n++) for(i=0;i<MB+1;i++)
{
Gnode[n].fh[i]=newrow(Gnode[n].fc[i]);
for(f=0;f<Gnode[n].fc[i];f++)
{
Gnode[n].fh[i][f].fv=0.; Gnode[n].fh[i][f].dfv=0.;
Gnode[n].fh[i][f].load=0.; Gnode[n].fh[i][f].typ=' ';
}
}
}
/**********************************************************/
void ready()
{
int n,i,r;
ROW *pr; NODE *pn;
for(n=0,pn=Gnode;n<MAXNOD;n++,pn++)
for(i=0;i<MB+1;i++)
for(r=0,pr=pn->fh[i];r<pn->fc[i];r++,pr++)
{ pr->load=0.; pr->dfv=0.; }
}
/**********************************************************/
void BD_vp_0() /* 2D Cartesian coor., plane channel*/
{
int i,n,r,o;
ROW *pr;
NODE *pn;
double x,x0,x1,x2,u0,u1,u2,sh2_coef,v0;
sh2_coef=3./sqrt(6.); v0=1.5;
for(i=0;i<MAXBD;i++)
{
switch(Btyp[i]){
case 'i': /*---inflow boundary , fully developed flow ---*/
if(Bd[i]<MAXVERT) /*assign the values */
{
x=Gnode[Bd[i]].x; pn=Gnode+Bd[i];
pr=pn->fh[0]; pr->typ='e'; pr->fv=0.;
pr=pn->fh[1]; pr->typ='e'; pr->fv=v0*(1.-x*x/4.);
}
else /* assign 2nd tangential derivatives */
{
x0=Gnode[Bd[i]].x; x1=Gnode[Bd[i-1]].x; x2=Gnode[Bd[i+1]].x;
pr=Gnode[Bd[i]].fh[0]; pr->typ='e'; pr->fv=0.;
u0=v0*(1.-x0*x0/4.);
u1=v0*(1.-x1*x1/4.);
u2=v0*(1.-x2*x2/4.);
pr=Gnode[Bd[i]].fh[1]; pr->typ='e'; pr->fv=(u1+u2-2.*u0)/sh2_coef;
}
break;
case 'w': /*-------wall boundary-------*/
{
pn=Gnode+Bd[i];
for(r=0,pr=pn->fh[0];r<pn->fc[0];r++,pr++) { pr->typ='e'; pr->fv=0.;}
for(r=0,pr=pn->fh[1];r<pn->fc[1];r++,pr++) { pr->typ='e'; pr->fv=0.;}
}
break;
case 'o': /*-------outflow boundary-------*/
if(Bd[i]<MAXVERT) /*assign the values */
{
x=Gnode[Bd[i]].x; pn=Gnode+Bd[i];
pr=pn->fh[0]; pr->typ='e'; pr->fv=0.;
pr=pn->fh[1]; pr->typ='e'; pr->fv=v0*(1.-x*x/4.);
}
else /* assign 2nd tangential derivatives */
{
x0=Gnode[Bd[i]].x; x1=Gnode[Bd[i-1]].x; x2=Gnode[Bd[i+1]].x;
pr=Gnode[Bd[i]].fh[0]; pr->typ='e'; pr->fv=0.;
u0=v0*(1.-x0*x0/4.);
u1=v0*(1.-x1*x1/4.);
u2=v0*(1.-x2*x2/4.);
pr=Gnode[Bd[i]].fh[1]; pr->typ='e'; pr->fv=(u1+u2-2*u0)/sh2_coef;
}
pr=pn->fh[2]; pr->typ='e'; pr->fv=0.; /*--- pressure boundary one datum */
break;
case 's': /*-------symmetric boundary-------*/
{
pn=Gnode+Bd[i];
for(r=0,pr=pn->fh[0];r<pn->fc[0];r++,pr++) { pr->typ='e'; pr->fv=0.;}
}
break;
case 'b': /*------- block boundary-------*/
{
pn=Gnode+Bd[i];
for(r=0,pr=pn->fh[0];r<pn->fc[0];r++,pr++) { pr->typ='e'; pr->fv=0.;}
for(r=0,pr=pn->fh[1];r<pn->fc[1];r++,pr++) { pr->typ='e'; pr->fv=0.;}
}
break;
default: printf("boundary type error!! Btyp=%c\n",Btyp[i]); exit(1);
} /** end of switch(btyp) **/
} /*end of s-loop*/
}
/**********************************************************/
void BD_vp_1() /* 2D cylindrical coor., tube problem */
{
int i,n,r,o;
ROW *pr;
NODE *pn;
double x,x0,x1,x2,u0,u1,u2,sh2_coef,v0,Q;
sh2_coef=3./sqrt(6.); v0=2.0; Q=PI;
for(i=0;i<MAXBD;i++)
{
switch(Btyp[i]){
case 'i': /*---inflow boundary , fully developed flow ---*/
if(Bd[i]<MAXVERT) /*assign the values */
{
x=Gnode[Bd[i]].x; pn=Gnode+Bd[i];
pr=pn->fh[0]; pr->typ='e'; pr->fv=0.;
pr=pn->fh[1]; pr->typ='e'; pr->fv=v0*(1.-x*x);
}
else /* assign 2nd tangential derivatives */
{
x0=Gnode[Bd[i]].x; x1=Gnode[Bd[i-1]].x; x2=Gnode[Bd[i+1]].x;
pr=Gnode[Bd[i]].fh[0]; pr->typ='e'; pr->fv=0.;
u0=v0*(1.-x0*x0);
u1=v0*(1.-x1*x1);
u2=v0*(1.-x2*x2);
pr=Gnode[Bd[i]].fh[1]; pr->typ='e'; pr->fv=(u1+u2-2*u0)/sh2_coef;
}
break;
case 'w': /*-------wall boundary-------*/
{
pn=Gnode+Bd[i];
for(r=0,pr=pn->fh[0];r<pn->fc[0];r++,pr++) { pr->typ='e'; pr->fv=0.;}
for(r=0,pr=pn->fh[1];r<pn->fc[1];r++,pr++) { pr->typ='e'; pr->fv=0.;}
}
break;
case 'o': /*-------outflow boundary-------*/
if(Bd[i]<MAXVERT) /*assign the values */
{
x=Gnode[Bd[i]].x; pn=Gnode+Bd[i];
pr=pn->fh[0]; pr->typ='e'; pr->fv=0.;
pr=pn->fh[1]; pr->typ='e'; pr->fv=v0*(1.-x*x);
}
else /* assign 2nd tangential derivatives */
{
x0=Gnode[Bd[i]].x; x1=Gnode[Bd[i-1]].x; x2=Gnode[Bd[i+1]].x;
pr=Gnode[Bd[i]].fh[0]; pr->typ='e'; pr->fv=0.;
u0=v0*(1.-x0*x0);
u1=v0*(1.-x1*x1);
u2=v0*(1.-x2*x2);
pr=Gnode[Bd[i]].fh[1]; pr->typ='e'; pr->fv=(u1+u2-2*u0)/sh2_coef;
}
pr=pn->fh[2]; pr->typ='e'; pr->fv=0.; /*--- pressure boundary one datum */
break;
case 's': /*-------symmetric boundary-------*/
{
pn=Gnode+Bd[i];
for(r=0,pr=pn->fh[0];r<pn->fc[0];r++,pr++) { pr->typ='e'; pr->fv=0.;}
}
break;
case 'b': /*------- block boundary-------*/
{
pn=Gnode+Bd[i];
for(r=0,pr=pn->fh[0];r<pn->fc[0];r++,pr++) { pr->typ='e'; pr->fv=0.;}
for(r=0,pr=pn->fh[1];r<pn->fc[1];r++,pr++) { pr->typ='e'; pr->fv=0.;}
}
break;
default: printf("boundary type error!! Btyp=%c\n",Btyp[i]); exit(1);
} /** end of switch(btyp) **/
} /*end of s-loop*/
}
/**********************************************************/
void read_data(in) char *in;
{
NODE *pn;
ROW *pr;
FILE *fp;
int m,n,r,i;
fp=fopen(in,"rb"); m=0;
m+=fread(&MAXNOD,Sint,1,fp);
m+=fread(&MAXELE,Sint,1,fp);
m+=fread(&MAXBD,Sint,1,fp);
m+=fread(&MAXVERT,Sint,1,fp);
m+=fread(Bd,Sint,MAXBD,fp);
m+=fread(Btyp,Schar,MAXBD,fp);
m+=fread(&Re0,Sdouble,1,fp);
for(n=0;n<MAXNOD;n++)
{
pn=Gnode+n;
m+=fread(pn,SNODE,1,fp);
for(i=0;i<MB+1;i++)
{
pr=pn->fh[i]=newrow(pn->fc[i]);
m+=fread(pr,SROW,pn->fc[i],fp);
}
}
m+=fread(Element,SELEM,MAXELE,fp);
fclose(fp);
fp=fopen("GAUSSL","rb");
m+=fread(Egauss,Slong,MAXELE*MO,fp);
fclose(fp);
}
/**********************************************************/
void write_data(out) char *out;
{
NODE *pn;
FILE *fp;
int m,n,i;
fp=fopen(out,"wb"); m=0;
m+=fwrite(&MAXNOD,Sint,1,fp);
m+=fwrite(&MAXELE,Sint,1,fp);
m+=fwrite(&MAXBD,Sint,1,fp);
m+=fwrite(&MAXVERT,Sint,1,fp);
m+=fwrite(Bd,Sint,MAXBD,fp);
m+=fwrite(Btyp,Schar,MAXBD,fp);
m+=fwrite(&Re,Sdouble,1,fp);
for(n=0;n<MAXNOD;n++)
{
pn=Gnode+n;
m+=fwrite(pn,SNODE,1,fp);
for(i=0;i<MB+1;i++) m+=fwrite(pn->fh[i],SROW,pn->fc[i],fp);
}
m+=fwrite(Element,SELEM,MAXELE,fp);
fclose(fp);
}
/**********************************************************/
void renew(out) char *out;
{
int en,ng,maxg;
ELEMENT e;
FILE *fp;
read_mesh();
initial();
Fg=fopen("GAUSS","wb");
for(en=0;en<MAXELE;en++)
{
for(ng=0;ng<MO;ng++) Egauss[en][ng]=-1;
e=Element[en]; maxg=Gnum(e,0,MB);
for(ng=3;ng<maxg+1;ng++) Gauss_write(en,ng,G(ng),W(ng));
}
fclose(Fg);
fp=fopen("GAUSSL","wb");
fwrite(Egauss,Slong,MAXELE*MO,fp);
fclose(fp);
write_data(out);
}
/**********************************************************/
void check_error(OB,OE) int OB,OE;
{
int nd,fd,i,omax;
NODE *pn;
ROW *pr;
for(pn=Gnode,nd=0;nd<MAXNOD;nd++,pn++)
for(i=OB;i<OE;i++)
for(fd=0;fd<pn->fc[i];fd++)
{
pr=pn->fh[i]+fd;
if(fabs(Maxerr)<fabs(pr->dfv)) { Maxerr=pr->dfv; omax=i; }
}
fprintf(Fprss,"max_err=%e o=%d ",Maxerr,omax);
printf("max_err=%e o=%d\n",Maxerr,omax);
}
/***********************************************************/
void md_fdm(OB,OE) int OB,OE;
{
int nd,fd,i;
NODE *pn;
for(nd=0;nd<MAXNOD;nd++)
{
pn=Gnode+nd;
for(i=OB;i<OE;i++)
for(fd=0;fd<pn->fc[i];fd++)
pn->fh[i][fd].fv+=pn->fh[i][fd].dfv;
}
}
/**********************************************************/
void Re_continuation(in,out) char *in,*out;
{
int itnum;
read_data(in);
for(;Re<Re_max;Re+=Re_step)
{
fprintf(Fprss,"---------------------------------------------\n");
fprintf(Fprss,"----- Reynolds Number=%f SUPG=%f ------\n",Re,supg);
fprintf(Fprss,"---------------------------------------------\n");
if(coor==0) BD_vp_0();
if(coor==1) BD_vp_1();
itnum=0;
do /*------ Newton iteration ------*/
{
ready(); Maxerr=0.;
solv(Eq);
md_fdm(0,3);
check_error(0,3);
itnum++;
fprintf(Fprss,"iteration:%d\n",itnum);
} while((itnum<N_FE)&&(fabs(Maxerr)>FE_ertol));
if(fabs(Maxerr)<FE_ertol) post(out);
else { printf("flow fails to converge!\n"); exit(1); }
}
}
/**********************************************************/
void assign_job(job) char *job;
{
FILE *ifp;
char s[31];
extern void post();
s[30]='\0';
ifp=fopen(job,"r");
fscanf(ifp,"%30c%s\n",s,Mesh); printf("%s%s\n",s,Mesh);
fscanf(ifp,"%30c%d\n",s,&coor); printf("%s%d\n",s,coor);
fscanf(ifp,"%30c%s\n",s,JOB); printf("%s%s\n",s,JOB);
fscanf(ifp,"%30c%s\n",s,D_in); printf("%s%s\n",s,D_in);
fscanf(ifp,"%30c%s\n",s,D_out); printf("%s%s\n",s,D_out);
fscanf(ifp,"%30c%s\n",s,C_sol); printf("%s%s\n",s,C_sol);
fscanf(ifp,"%30c%c\n",s,&Eq); printf("%s%c\n",s,Eq);
fscanf(ifp,"%30c%lf\n",s,&Re); printf("%s %e\n",s,Re);
fscanf(ifp,"%30c%lf\n",s,&Re_step); printf("%s%e\n",s,Re_step);
fscanf(ifp,"%30c%lf\n",s,&Re_max); printf("%s%e\n",s,Re_max);
fscanf(ifp,"%30c%lf\n",s,&FE_ertol); printf("%s%e\n",s,FE_ertol);
fscanf(ifp,"%30c%d\n",s,&N_FE); printf("%s%d\n",s,N_FE);
fscanf(ifp,"-------------------------------------------\n");
printf("-------------------------------------------\n");
fclose(ifp);
fprintf(Fprss,"\n---------------------------------------------\n");
fprintf(Fprss,"------------ JOB: %s -----------\n",JOB);
fprintf(Fprss,"---------------------------------------------\n");
if(strcmp(JOB,"renew")==0) { renew(D_in); exit(1); }
if(strcmp(JOB,"only_post")==0) { only_post(D_in); exit(1); }
if(strcmp(JOB,"Re_continuation")==0)
{ Re_continuation(D_in,D_out); exit(1); }
}
/**********************************************************/
main()
{
int itnum,n,i;
Fprss=fopen("prss","a");
assign_job("job");
fclose(Fprss);
}
/**********************************************************/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -