f3.c

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

/****************************** f3.c *************************/

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

extern NODE       Gnode[];
extern ELEMENT    Element[];
extern long       Egauss[][MO];
extern int        MAXNOD,MAXELE;
extern FILE       *Fprss,*Fg;
/**********************************************************/
extern double shape(), dsh_k(), dsh_a();
extern double X(),Y(),Xk(),Xa(),Yk(),Ya();
extern int  I(),J();  
extern void add_list();  
/*********************************************************/
double G2[2]={-0.577350269189626, 0.577350269189626};
double W2[2]={ 1.000000000000000, 1.000000000000000};
double G3[3]=
{-0.774596669241483,  0.000000000000000,  0.774596669241483 };
double W3[3]=
 {0.555555555555555,  0.888888888888889,  0.555555555555555 };
double G4[4]=
{-0.861136311594053, -0.339981043584856,  0.339981043584856,  0.861136311594053};
double W4[4]=
 {0.347854845137454,  0.652145154862546,  0.652145154862546,  0.347854845137454};
double G5[5]=
{-0.906179845938664, -0.538469310105683,  0.000000000000000,
  0.538469310105683,  0.906179845938664 };
double W5[5]=
 {0.236926885056189,  0.478628670499366,  0.568888888888889,
  0.478628670499366,  0.236926885056189 };
double G6[6]=
{-0.932469514203152, -0.661209386466265, -0.238619186083197,
  0.238619186083197,  0.661209386466265,  0.932469514203152 };
double W6[6]=
 {0.171324492379170,  0.360761573048139,  0.467913934572691,
  0.467913934572691,  0.360761573048139,  0.171324492379170 };
/*********************************************************/
/*********************************************************/
double x[MG],y[MG],xk[MG],xa[MG],yk[MG],ya[MG],Jac[MG],WJ[MG];
double v1[MG],v2[MG],v11[MG],v12[MG],v21[MG],v22[MG],p[MG];
double sf1[MG],sf2[MG];
double Sh[9][MF][MO][MG],Sh1[9][MF][MO][MG],Sh2[9][MF][MO][MG];
/*********************************************************/
/*********************************************************/
double *G(maxo) int maxo;
{
switch(maxo) {
case 0: printf("Element order Error !!\n"); exit(1);
case 2: return G2;
case 3: return G3;
case 4: return G4;
case 5: return G5;
case 6: return G6;
default: printf("Element order too high!!\n"); exit(1);
             }
}
/*********************************************************/
double *W(maxo) int maxo;
{
switch(maxo) {
case 0: printf("Element order Error !!\n"); exit(1);
case 2: return W2;
case 3: return W3;
case 4: return W4;
case 5: return W5;
case 6: return W6;
default: printf("Element order too high!!\n"); exit(1);
             }
}
/*********************************************************/
/*********************************************************/
int Gnum(e) ELEMENT e;
{
int n,i,maxo;
NODE *pn;

for(maxo=0,n=4;n<8;n++) 
{
pn=Gnode+e.n[n];
for(i=0;i<2;i++) if(maxo<pn->fc[i]) maxo=pn->fc[i];
}
return(maxo+2);
}
/*********************************************************/
/*********************************************************/
double dsh_x(n,f,eo,gp,g) int n,f,eo,g; PAIR gp;
{
 extern double dsh_k(),dsh_a(),ya[],yk[],Jac[];
 return((dsh_k(n,f,eo,gp)*ya[g]-dsh_a(n,f,eo,gp)*yk[g])/Jac[g]);
}
/*********************************************************/
double dsh_y(n,f,eo,gp,g) int n,f,eo,g; PAIR gp;
{
 extern double dsh_k(),dsh_a(),xa[],xk[],Jac[];
 return((dsh_a(n,f,eo,gp)*xk[g]-dsh_k(n,f,eo,gp)*xa[g])/Jac[g]);
}
/*********************************************************/
void Gauss_write(en,maxg,GS,WT)
     int en,maxg; double GS[],WT[];
{
ELEMENT e;
PAIR   gp[MG];
int i,n,g,maxg2,eo,maxo,f;
double nodx[9],nody[9];
long m=0;

e=Element[en];
for(i=0;i<9;i++)
{ n=e.n[i]; nodx[i]=Gnode[n].x; nody[i]=Gnode[n].y; }

maxg2=maxg*maxg;
for(g=0;g<maxg2;g++)
{
gp[g].k=GS[g/maxg]; gp[g].a=GS[g%maxg];
x[g]=X(nodx,gp[g]); y[g]=Y(nody,gp[g]);
xk[g]=Xk(nodx,gp[g]);    xa[g]=Xa(nodx,gp[g]);
yk[g]=Yk(nody,gp[g]);    ya[g]=Ya(nody,gp[g]);
Jac[g]=xk[g]*ya[g]-xa[g]*yk[g];
if(Jac[g]<0.) {printf("Negative Jac=%f\n",Jac[g]);
for(i=0;i<9;i++)
printf("%d\t",Element[en].n[i]);
printf("\n");
exit(1);}
}

for(g=0;g<maxg2;g++) { WJ[g]=WT[g/maxg]*WT[g%maxg]*Jac[g]; }

maxo=maxg-1;
for(n=0;n<4;n++) for(eo=0;eo<maxo;eo++)
for(g=0;g<maxg2;g++) 
 { 
 Sh[n][0][eo][g]=shape(n,0,eo+1,gp[g]);
 Sh1[n][0][eo][g]=dsh_x(n,0,eo+1,gp[g],g);
 Sh2[n][0][eo][g]=dsh_y(n,0,eo+1,gp[g],g);
 }

for(n=4;n<8;n++) for(f=0;f<maxg-2;f++) for(eo=1;eo<maxo;eo++)
for(g=0;g<maxg2;g++) 
 { 
 Sh[n][f][eo][g]=shape(n,f,eo+1,gp[g]);
 Sh1[n][f][eo][g]=dsh_x(n,f,eo+1,gp[g],g);
 Sh2[n][f][eo][g]=dsh_y(n,f,eo+1,gp[g],g);
 }

for(eo=1;eo<maxo;eo++) for(f=0;f<eo*eo;f++)
for(g=0;g<maxg2;g++) 
 { 
 Sh[8][f][eo][g]=shape(8,f,eo+1,gp[g]);
 Sh1[8][f][eo][g]=dsh_x(8,f,eo+1,gp[g],g);
 Sh2[8][f][eo][g]=dsh_y(8,f,eo+1,gp[g],g);
 }

Egauss[en][maxg-3]=ftell(Fg);
m+=fwrite(x,Sdouble,maxg2,Fg);
m+=fwrite(y,Sdouble,maxg2,Fg);
m+=fwrite(Jac,Sdouble,maxg2,Fg);
m+=fwrite(WJ,Sdouble,maxg2,Fg);
for(n=0;n<4;n++) for(eo=0;eo<maxo;eo++)
 { 
 m+=fwrite(Sh[n][0][eo],Sdouble,maxg2,Fg);
 m+=fwrite(Sh1[n][0][eo],Sdouble,maxg2,Fg);
 m+=fwrite(Sh2[n][0][eo],Sdouble,maxg2,Fg);
 }
for(n=4;n<8;n++) for(f=0;f<maxg-2;f++) for(eo=1;eo<maxo;eo++)
 { 
 m+=fwrite(Sh[n][f][eo],Sdouble,maxg2,Fg);
 m+=fwrite(Sh1[n][f][eo],Sdouble,maxg2,Fg);
 m+=fwrite(Sh2[n][f][eo],Sdouble,maxg2,Fg);
 }
for(eo=1;eo<maxo;eo++) for(f=0;f<eo*eo;f++)
 { 
 m+=fwrite(Sh[8][f][eo],Sdouble,maxg2,Fg);
 m+=fwrite(Sh1[8][f][eo],Sdouble,maxg2,Fg);
 m+=fwrite(Sh2[8][f][eo],Sdouble,maxg2,Fg);
 }
}
/*********************************************************/
void Gauss_read(en,maxg)
int en,maxg;
{
int maxg2,n,f,eo,maxo;  long m=0;

maxg2=maxg*maxg;
maxo=maxg-1;

fseek(Fg,Egauss[en][maxg-3],0);
m+=fread(x,Sdouble,maxg2,Fg);
m+=fread(y,Sdouble,maxg2,Fg);
m+=fread(Jac,Sdouble,maxg2,Fg);
m+=fread(WJ,Sdouble,maxg2,Fg);
for(n=0;n<4;n++) for(eo=0;eo<maxo;eo++)
 { 
 m+=fread(Sh[n][0][eo],Sdouble,maxg2,Fg);
 m+=fread(Sh1[n][0][eo],Sdouble,maxg2,Fg);
 m+=fread(Sh2[n][0][eo],Sdouble,maxg2,Fg);
 }
for(n=4;n<8;n++) for(f=0;f<maxg-2;f++) for(eo=1;eo<maxo;eo++)
 { 
 m+=fread(Sh[n][f][eo],Sdouble,maxg2,Fg);
 m+=fread(Sh1[n][f][eo],Sdouble,maxg2,Fg);
 m+=fread(Sh2[n][f][eo],Sdouble,maxg2,Fg);
 }
for(eo=1;eo<maxo;eo++) for(f=0;f<eo*eo;f++)
 { 
 m+=fread(Sh[8][f][eo],Sdouble,maxg2,Fg);
 m+=fread(Sh1[8][f][eo],Sdouble,maxg2,Fg);
 m+=fread(Sh2[8][f][eo],Sdouble,maxg2,Fg);
 }
}
/*********************************************************/
/*********************************************************/
void veloc_press(e,maxg)
ELEMENT e;
int maxg;
{
int g,n,fd,eo;
double fv;

for(g=0;g<maxg;g++) {
eo=e.vo-1;
for(v1[g]=v11[g]=v12[g]=0.,n=0;n<9;n++)
 for(fd=0;fd<Gnode[e.n[n]].fc[0];fd++)
 { 
 fv=Gnode[e.n[n]].fh[0][fd].fv;
 v1[g]+=fv*Sh[n][fd][eo][g];  
 v11[g]+=fv*Sh1[n][fd][eo][g];
 v12[g]+=fv*Sh2[n][fd][eo][g];
 }
for(v2[g]=v21[g]=v22[g]=0.,n=0;n<9;n++)
 for(fd=0;fd<Gnode[e.n[n]].fc[1];fd++)
 { 
 fv=Gnode[e.n[n]].fh[1][fd].fv;
 v2[g]+=fv*Sh[n][fd][eo][g];  
 v21[g]+=fv*Sh1[n][fd][eo][g];  
 v22[g]+=fv*Sh2[n][fd][eo][g];  
 }
eo=e.po-1;
for(p[g]=0.,n=0;n<9;n++)
 for(fd=0;fd<Gnode[e.n[n]].fc[2];fd++)
 { 
 fv=Gnode[e.n[n]].fh[2][fd].fv;
 p[g]+=fv*Sh[n][fd][eo][g];  
 }
                     }
}
/*********************************************************/
void veloc_sf(e,maxg)
ELEMENT e; int maxg;
{
int g,n,fd,eo;
double fv;

for(g=0;g<maxg;g++) {

eo=e.vo-1;
for(v1[g]=v11[g]=v12[g]=0.,n=0;n<9;n++)
 for(fd=0;fd<Gnode[e.n[n]].fc[0];fd++)
 { 
 fv=Gnode[e.n[n]].fh[0][fd].fv;
 v1[g]+=fv*Sh[n][fd][eo][g];  
 v11[g]+=fv*Sh1[n][fd][eo][g];
 v12[g]+=fv*Sh2[n][fd][eo][g];
 }
for(v2[g]=v21[g]=v22[g]=0.,n=0;n<9;n++)
 for(fd=0;fd<Gnode[e.n[n]].fc[1];fd++)
 { 
 fv=Gnode[e.n[n]].fh[1][fd].fv;
 v2[g]+=fv*Sh[n][fd][eo][g];  
 v21[g]+=fv*Sh1[n][fd][eo][g];  
 v22[g]+=fv*Sh2[n][fd][eo][g];  
 }
eo=2-1;
for(sf1[g]=sf2[g]=0.,n=0;n<9;n++) 
for(fd=0;fd<Gnode[e.n[n]].fc[MB];fd++)
 { 
 fv=Gnode[e.n[n]].fh[MB][fd].fv;