f6.c

本源代码实现了通过网格插值的方法来进行有限元流线处理的结果

#include <stdio.h>
#include <math.h>
#include "com.h"

int               MAXO;
extern char       Eq,C_sol[],Btyp[];
extern NODE       Gnode[];
extern double     Re;
extern ELEMENT    Element[];
extern int        MAXNOD,MAXELE,MAXVERT,MAXBD,Bd[];
extern FILE       *Fprss;

/**********************************************************/
extern double shape(),dsh_k(),dsh_a(),X(),Y(),Xk(),Xa(),Yk(),Ya();
/**********************************************************/
double nod_value(e,eo,nod,o) ELEMENT e; int eo,nod,o;
{
NODE    *pn;
ROW     *pr;
int     n,fd;
double   val;
PAIR    gp;

switch(nod)
{
case 0:   gp.k=-1.0;  gp.a=-1.0;  break;
case 1:   gp.k= 1.0;  gp.a=-1.0;  break;
case 2:   gp.k= 1.0;  gp.a= 1.0;  break;
case 3:   gp.k=-1.0;  gp.a= 1.0;  break;
case 4:   gp.k= 0.0;  gp.a=-1.0;  break;
case 5:   gp.k= 1.0;  gp.a= 0.0;  break;
case 6:   gp.k= 0.0;  gp.a= 1.0;  break;
case 7:   gp.k=-1.0;  gp.a= 0.0;  break;
case 8:   gp.k= 0.0;  gp.a= 0.0;  break;
}

for(val=0.,n=0;n<9;n++)
{
 pn=Gnode+e.n[n];
 for(fd=0,pr=pn->fh[o];fd<pn->fc[o];fd++,pr++) 
 val+=pr->fv*shape(n,fd,eo,gp);
}
return(val);
}
/**********************************************************/
void tecplot()
{
ELEMENT e;
double v1[MMAXNOD],v2[MMAXNOD],p[MMAXNOD],sf[MMAXNOD],x,y;
double length,l,m;
int  i,j,n,en;
char str[35];
FILE *fp,*ofp;

for(i=0;i<MAXELE;i++) for(n=0;n<9;n++)
{
e=Element[i];
v1[e.n[n]]=nod_value(e,e.vo,n,0);
v2[e.n[n]]=nod_value(e,e.vo,n,1);
p[e.n[n]]=nod_value(e,e.po,n,2);
sf[e.n[n]]=nod_value(e,2,n,MB);
}

fp=fopen(C_sol,"w");
fprintf(fp,"title=\"De=%f\"\n",Re);
fprintf(fp,"variables=\"x\",\"y\",\"Vx\",\"Vy\",\"p\",\"sf\"\n");
fprintf(fp,"zone N=%d, E=%d, F=FEPOINT, ET=QUADRILATERAL\n",MAXNOD,MAXELE*4);
for(n=0;n<MAXNOD;n++)
{
 x=Gnode[n].x; y=Gnode[n].y;
 fprintf(fp,"%e %e %e %e %e %e\n",x,y,v1[n],v2[n],p[n],sf[n]);
}
fprintf(fp,"#####################################\n");
for(en=0;en<MAXELE;en++)
{
e=Element[en];
fprintf(fp,"%d %d %d %d\n",e.n[0]+1,e.n[4]+1,e.n[8]+1,e.n[7]+1);
fprintf(fp,"%d %d %d %d\n",e.n[4]+1,e.n[1]+1,e.n[5]+1,e.n[8]+1);
fprintf(fp,"%d %d %d %d\n",e.n[8]+1,e.n[5]+1,e.n[2]+1,e.n[6]+1);
fprintf(fp,"%d %d %d %d\n",e.n[7]+1,e.n[8]+1,e.n[6]+1,e.n[3]+1);
}
fclose(fp);
}
/**********************************************************/
void BD_sf_0()
{

int i,n,r,o;
ROW  *pr;
NODE *pn;
double v0,Q;

v0=1.5;  Q=1.0;
                                                                                
for(i=0;i<MAXBD;i++)
{

switch(Btyp[i]){ 

 case 'i':           /*---inflow boundary , fully developed flow ---*/
 break;
 case 'w':           /*-------wall boundary-------*/
 
 {
  pn=Gnode+Bd[i];
  for(r=0,pr=pn->fh[MB];r<pn->fc[MB];r++,pr++) pr->typ='e';
  if(Bd[i]<MAXVERT) pn->fh[MB]->fv=Q; else pn->fh[MB]->fv=0.;
 }
 break;
 case 'o':           /*-------outflow boundary-------*/
 break;
 case 's':           /*-------symmetric boundary-------*/

 {
  pn=Gnode+Bd[i];
  for(r=0,pr=pn->fh[MB];r<pn->fc[MB];r++,pr++) { pr->typ='e'; pr->fv=0.;}
 }

 break;
 case 'b':           /*------- block boundary-------*/

 {
  pn=Gnode+Bd[i];
  for(r=0,pr=pn->fh[MB];r<pn->fc[MB];r++,pr++) pr->typ='e';
  if(Bd[i]<MAXVERT) pn->fh[MB]->fv=Q; else pn->fh[MB]->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_sf_1()
{

int i,n,r,o;
ROW  *pr;
NODE *pn;
double v0,Q;

v0=2.0;  Q=PI;
                                                                                
for(i=0;i<MAXBD;i++)
{
switch(Btyp[i]){ 

 case 'i':           /*---inflow boundary , fully developed flow ---*/
 break;
 case 'w':           /*-------wall boundary-------*/
 
 {
  pn=Gnode+Bd[i];
  for(r=0,pr=pn->fh[MB];r<pn->fc[MB];r++,pr++) pr->typ='e';
  if(Bd[i]<MAXVERT) pn->fh[MB]->fv=Q; else pn->fh[MB]->fv=0.;
 }
 break;
 case 'o':           /*-------outflow boundary-------*/
 break;
 case 's':           /*-------symmetric boundary-------*/
 
 {
  pn=Gnode+Bd[i];
  for(r=0,pr=pn->fh[MB];r<pn->fc[MB];r++,pr++) { pr->typ='e'; pr->fv=0.;}
 }

 break;
 case 'b':           /*------- block boundary-------*/
 
 {
  pn=Gnode+Bd[i];
  for(r=0,pr=pn->fh[MB];r<pn->fc[MB];r++,pr++) pr->typ='e';
  if(Bd[i]<MAXVERT) pn->fh[MB]->fv=Q; else pn->fh[MB]->fv=0.;
 }

 break;
 default: printf("boundary type error!! Btyp=%c\n",Btyp[i]); exit(1);
                  }  /** end of switch(btyp) **/
} /*end of s-loop*/

}
/**********************************************************/
void post(out) char *out;
{
extern void boundary0(),boundary1(),write_data(),solv(),md_fdm();
extern double Re,supg;
extern int coor;

fprintf(Fprss,"------ Post Processing ------\n");
fprintf(Fprss,"--- Re=%f SUPG=%f ---\n",Re,supg);
if(coor==0) BD_sf_0();
if(coor==1) BD_sf_1();
solv('S');
md_fdm(3,4);
tecplot();
write_data(out);
}
/**********************************************************/
void only_post(in) char *in;
{
extern void read_data();
extern double Re,supg;

read_data(in);
fprintf(Fprss,"------ Processing: %s  Re=%f SUPG=%f ----\n",in,Re,supg);
tecplot();
}
/**********************************************************/