f3.c

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

 sf1[g]+=fv*Sh1[n][fd][eo][g];  
 sf2[g]+=fv*Sh2[n][fd][eo][g];  
 }
                     }
}
/*********************************************************/
/*********************************************************/
void sl_Newt0(e,maxg) /* For planar flows */
ELEMENT e; int maxg;
{
NODE   *pn;
ROW    *pr;
NOF     nofi,nofj;
int     g,wn,sn,wf,i,j,eo,n,f;
double wt[MG],wt1[MG],wt2[MG];
double sh,sh1,sh2,res,coef;
extern double Re;

/*-------- For v1_weighting residuals-----------*/
                   nofi.o=0;
for(wn=0;wn<9;wn++)
{ 
                   nofi.n=e.n[wn];
  pn=Gnode+e.n[wn];
  for(wf=0,pr=pn->fh[0];wf<pn->fc[0];wf++,pr++) if(pr->typ==' ')
  { 
  eo=e.vo-1;
                   nofi.f=wf;
  for(res=0.,g=0;g<maxg;g++) 
  {
   wt[g]=Sh[wn][wf][eo][g];
   wt1[g]=Sh1[wn][wf][eo][g];
   wt2[g]=Sh2[wn][wf][eo][g];
   res+=Re*wt[g]*(v1[g]*v11[g]+v2[g]*v12[g])*WJ[g];
   res+=(2.*wt1[g]*v11[g]+wt2[g]*(v12[g]+v21[g])-p[g]*wt1[g])*WJ[g];
  }
   pr->load-=res;
   i=I(nofi);
                               /*--- v1_weight, v1_shape ---*/
                                 nofj.o=0;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[0];f++) if(Gnode[nofj.n].fh[0][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh=Sh[sn][f][eo][g];
      sh1=Sh1[sn][f][eo][g];
      sh2=Sh2[sn][f][eo][g];
      coef+=Re*wt[g]*(sh*v11[g]+v1[g]*sh1+v2[g]*sh2)*WJ[g];
      coef+=(2.*wt1[g]*sh1+wt2[g]*sh2)*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                 }}
                               /*--- v1_weight, v2_shape ---*/
                                 nofj.o=1;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[1];f++) if(Gnode[nofj.n].fh[1][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh=Sh[sn][f][eo][g];
      sh1=Sh1[sn][f][eo][g];
      coef+=Re*wt[g]*sh*v12[g]*WJ[g];
      coef+=wt2[g]*sh1*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                 }}
                               /*--- v1_weight, p_shape ---*/
   eo=e.po-1;
                                 nofj.o=2;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[2];f++) if(Gnode[nofj.n].fh[2][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh=Sh[sn][f][eo][g];
      coef-=(wt1[g]*sh)*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                  }}
  }
}

/*-------- For v2_weighting residuals-----------*/
                   nofi.o=1;
for(wn=0;wn<9;wn++)
{
                   nofi.n=e.n[wn];
  pn=Gnode+e.n[wn];
  for(wf=0,pr=pn->fh[1];wf<pn->fc[1];wf++,pr++) if(pr->typ==' ')
  {
  eo=e.vo-1;
                   nofi.f=wf;
  for(res=0.,g=0;g<maxg;g++) 
  {
   wt[g]=Sh[wn][wf][eo][g];	  
   wt1[g]=Sh1[wn][wf][eo][g];
   wt2[g]=Sh2[wn][wf][eo][g];
   res+=Re*wt[g]*(v1[g]*v21[g]+v2[g]*v22[g])*WJ[g];
   res+=(2.*wt2[g]*v22[g]+wt1[g]*(v21[g]+v12[g])-wt2[g]*p[g])*WJ[g];
  }
   pr->load-=res;
   i=I(nofi);
                               /*--- v2_weight, v1_shape ---*/
                                 nofj.o=0;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[0];f++) if(Gnode[nofj.n].fh[0][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh=Sh[sn][f][eo][g];
      sh2=Sh2[sn][f][eo][g];
      coef+=Re*wt[g]*sh*v21[g]*WJ[g];
      coef+=wt1[g]*sh2*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                 }}
                               /*--- v2_weight, v2_shape ---*/
                                 nofj.o=1;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[1];f++) if(Gnode[nofj.n].fh[1][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh=Sh[sn][f][eo][g];
      sh1=Sh1[sn][f][eo][g];
      sh2=Sh2[sn][f][eo][g];
      coef+=Re*wt[g]*(v1[g]*sh1+sh*v22[g]+v2[g]*sh2)*WJ[g];
      coef+=(wt1[g]*sh1+2.*wt2[g]*sh2)*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                 }}
                               /*--- v2_weight, p_shape ---*/
   eo=e.po-1;
                                 nofj.o=2;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[2];f++) if(Gnode[nofj.n].fh[2][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh=Sh[sn][f][eo][g];
      coef-=(wt2[g]*sh)*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                  }}
  }
}

/*-------- For p_weighting residuals-----------*/
                   nofi.o=2;
for(wn=0;wn<9;wn++)
{
                   nofi.n=e.n[wn];
  pn=Gnode+e.n[wn];
  for(wf=0,pr=pn->fh[2];wf<pn->fc[2];wf++,pr++) if(pr->typ==' ')  
  { 
  eo=e.po-1;
                   nofi.f=wf;
  for(res=0.,g=0;g<maxg;g++) 
  {
   wt[g]=Sh[wn][wf][eo][g];
   res-=wt[g]*(v11[g]+v22[g])*WJ[g];
  }
   pr->load-=res;
   i=I(nofi);
                               /*--- p_weight, v1_shape ---*/
   eo=e.vo-1;
                                 nofj.o=0;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[0];f++) if(Gnode[nofj.n].fh[0][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh=Sh[sn][f][eo][g];
      sh1=Sh1[sn][f][eo][g];
      coef-=wt[g]*sh1*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                 }}
                               /*--- p_weight, v2_shape ---*/
                                 nofj.o=1;
   for(sn=0;sn<9;sn++)         { nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[1];f++) if(Gnode[nofj.n].fh[1][f].typ==' ')
   { nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh2=Sh2[sn][f][eo][g];
      coef-=wt[g]*sh2*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                  }}
  }
}
}
/*********************************************************/
void sl_sf0(e,maxg)
ELEMENT e; int maxg;
{
NOF    nofi,nofj;
NODE   *pn;
ROW    *pr;
int g,wn,sn,eo,wf,f,i,j;
double wt[MG],wt1[MG],wt2[MG];
double sh1,sh2;
double res,coef;


    /*---- For sf_source terms ----*/
                      nofi.o=MB;
for(wn=0;wn<9;wn++)
{ 
                      nofi.n=e.n[wn];
  pn=Gnode+e.n[wn];
  for(wf=0,pr=pn->fh[MB];wf<pn->fc[MB];wf++,pr++) if(pr->typ==' ')
  { 
  eo=2-1;           nofi.f=wf;
  for(res=0.,g=0;g<maxg;g++) 
  {
   wt[g]=Sh[wn][wf][eo][g];
   wt1[g]=Sh1[wn][wf][eo][g];
   wt2[g]=Sh2[wn][wf][eo][g];
   res+=((sf1[g]*wt1[g]+sf2[g]*wt2[g])+wt[g]*(v21[g]-v12[g]))*WJ[g];
  }
  pr->load-=res;
  i=I(nofi);
       /*--- sf_stiff matrix ---*/
                  nofj.o=MB;
   for(sn=0;sn<9;sn++)        {  nofj.n=e.n[sn];
   for(f=0;f<Gnode[nofj.n].fc[MB];f++) if(Gnode[nofj.n].fh[MB][f].typ==' ')
   {      nofj.f=f; 
   for(coef=0.,g=0;g<maxg;g++) 
   {
      sh1=Sh1[sn][f][eo][g];
      sh2=Sh2[sn][f][eo][g];
      coef+=(sh1*wt1[g]+sh2*wt2[g])*WJ[g];
   }
      j=J(nofj);
      add_list(i,j,coef);
                                        }}
  }
}

}
/*********************************************************/
/*********************************************************/
void asm_Newt(en)  int en;
{
ELEMENT e;
int ng;

e=Element[en];
ng=Gnum(e);
Gauss_read(en,ng);
veloc_press(e,ng*ng);
sl_Newt0(e,ng*ng);
}
/*********************************************************/
void asm_sf(en)  int en;
{
ELEMENT e;
int ng;

e=Element[en];
ng=Gnum(e);
Gauss_read(en,ng);
veloc_sf(e,ng*ng);
sl_sf0(e,ng*ng);
}
/*********************************************************/