f1.c

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

/******************************** f1.c ************************/

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

/**********************************************************/
extern NODE       Gnode[MMAXNOD];
extern ELEMENT    Element[];
extern int        MAXNOD,MAXELE,MAXVERT;
extern FILE       *Fprss,*Fg;
/**********************************************************/
#define MAXWD 300
double  Fmat[MAXWD][MAXWD],MAXpiv,MINpiv;
NOF     Fi[MAXWD],Fj[MAXWD];
int     Lasti,Lastj;
int     Last_e[MMAXNOD],Maxwd;
char    ntyp[9];

PIV   *Pivot;
FILE  *Fpiv,*fpiv;
long       Mem;
/**********************************************************/
/**********************************************************/
ROW *newrow(count)
int count;
{
  ROW *pr;

  if(count==0) return NULL;
  pr=(ROW*)calloc(count,SROW);
  if(pr!=NULL)  { Mem+=SROW*count; return(pr); }
  else {printf("newrow failed\n"); exit(1);}
}
/**********************************************************/
PIV *newpiv()
{
  PIV *pp;
  pp=(PIV*)malloc(SPIV);
  if(pp!=NULL)  { Mem+=SPIV;  return(pp);}
  else {printf("newpiv failed\n"); exit(1);}

}
/**********************************************************/
PIV *nextpiv(current) PIV *current;
{
PIV *piv;

piv=newpiv(); piv->next=current;
return piv;
}
/**********************************************************/
int nof_mat(f1,f2) NOF f1,f2;
{
if((f1.n==f2.n)&&(f1.o==f2.o)&&(f1.f==f2.f))
      return(1);
else  return(0);
}
/**********************************************************/
int I(nof) NOF nof;
{
int i;

for(i=0;i<Lasti;i++) if(nof_mat(Fi[i],nof)) return i;

Fi[Lasti]=nof; Lasti++;
if(Lasti>Maxwd) Maxwd=Lasti;
if(Lasti>=MAXWD) {printf("MAXWD too small!!!\n"); exit(1);}
return i;
}
/**********************************************************/
int J(nof) NOF nof;
{
int j;

for(j=0;j<Lastj;j++)
if(nof_mat(Fj[j],nof)) return j;

Fj[Lastj]=nof; Lastj++;
if(Lastj>Maxwd) Maxwd=Lastj;
if(Lastj>=MAXWD) {printf("MAXWD too small!!!\n"); exit(1);}
return j;
}
/**********************************************************/
void add_list(i,j,cf) int i,j; double cf;
{
if(ABS(cf)<ZERO) return;
Fmat[i][j]+=cf;
return;
}
/**********************************************************/
void adjust_front(ip,jp)
int ip,jp;
{
int i,j,L1;

L1=Lasti-1;
Fi[ip]=Fi[L1];
for(j=0;j<Lastj;j++)
{ Fmat[ip][j]=Fmat[L1][j]; Fmat[L1][j]=0.; }

L1=Lastj-1;
Fj[jp]=Fj[L1];
for(i=0;i<Lasti;i++)
{ Fmat[i][jp]=Fmat[i][L1]; Fmat[i][L1]=0.; }

Fi[Lasti-1].n=-2; Fj[Lastj-1].n=-2;
Lasti--;          Lastj--;
}
/*******************************************************/
void write_piv(ip,jp) 
int  ip,jp;
{
int j; long  m;

Pivot->begin=ftell(Fpiv);
m=0;
m+=fwrite(&(Fj[jp]),SNOF,1,Fpiv);
m+=fwrite(&(Fmat[ip][jp]),Sdouble,1,Fpiv);
for(j=0;j<Lastj;j++) 
if((j!=jp)&&(ABS(Fmat[ip][j])>ZERO))
{
m+=fwrite(&(Fj[j]),SNOF,1,Fpiv);
m+=fwrite(&(Fmat[ip][j]),Sdouble,1,Fpiv);
} 
Pivot->end=ftell(Fpiv);
}
/**********************************************************/
void eli_piv(ip,jp)
int ip,jp;
{
ROW   *pr;
double scale;
int i,j;

for(i=0;i<Lasti;i++) 
if((i!=ip)&&(ABS(Fmat[i][jp])>ZERO))
{
  scale=Fmat[i][jp];
  for(j=0;j<Lastj;j++)
     if((j!=jp)&&(ABS(Fmat[ip][j])>ZERO))
        add_list(i,j,scale*Fmat[ip][j]);

  pr=Gnode[Fi[i].n].fh[Fi[i].o]+Fi[i].f;
  pr->load-=scale*Fmat[ip][jp];
}

}
/*********************************************************/
void mk_piv(ip,jp)
int ip,jp;
{
ROW *pr;
double val;
int j;

val=Fmat[ip][jp];
for(j=0;j<Lastj;j++)
   if(ABS(Fmat[ip][j])>ZERO) Fmat[ip][j]/=(-val);

pr=Gnode[Fi[ip].n].fh[Fi[ip].o]+Fi[ip].f;
Fmat[ip][jp]=pr->load/val;
}

/**********************************************************/
int cond(nod,o) int nod,o;
{
NOF nof;
int   f,i,j,ip,jp,count;
double coef;

count=0;
nof.n=nod;
nof.o=o;
for(f=0;f<Gnode[nod].fc[o];f++) if(Gnode[nod].fh[o][f].typ-'e')
 {
   nof.f=f;
   ip=jp=-1;
   for(i=0;i<Lasti;i++) if(nof_mat(Fi[i],nof)) {ip=i; break;}
   for(j=0;j<Lastj;j++) if(nof_mat(Fj[j],nof)) {jp=j; break;}

  if((ABS(Fmat[ip][jp])>ZERO)&&(ip>=0)&&(jp>=0))
  { 
   coef=Fmat[ip][jp];
   mk_piv(ip,jp);
   eli_piv(ip,jp);
   write_piv(ip,jp);
   adjust_front(ip,jp);
   Pivot=nextpiv(Pivot);
   if(ABS(coef)>ABS(MAXpiv)) MAXpiv=coef;
   if(ABS(coef)<ABS(MINpiv)) MINpiv=coef;
   count++;
  }
  else {
        printf("n=%d  o=%d  f=%d \n",nof.n,nof.o,nof.f);
        printf("ip=%d  jp=%d  Lasti=%d  Lastj=%d\n",ip,jp,Lasti,Lastj);
        printf(" Fails to condense f.!\n"); exit(1);
       }
 }
return count;
}
/**********************************************************/
int condense(e,Eq) ELEMENT e; char Eq;
{
int   o,nod,count,n,j;

count=0;
if(Eq=='S')  /* solve the stream function */
{
 for(n=8;n>-1;n--) if(ntyp[n]=='f')
 { nod=e.n[n]; count+=cond(nod,3);}
 return(count);
}

if(Eq=='N')  /* solve Newtonian problems */
{
 for(n=8;n>-1;n--) if(ntyp[n]=='f')
 { nod=e.n[n]; for(o=0;o<3;o++) count+=cond(nod,o); }
 return(count);
}
}
/**********************************************************/
void backsub()
{
double cf;
NOF col,*pc;
ROW   *pr1,*pr2;
FILE  *fp;
PIV *head;

fp=fopen("piv_dat","rb");
do
{
 fseek(fp,Pivot->begin,0);
 fread(&(col),SNOF,1,fp); fread(&(cf),Sdouble,1,fp);
 pr1=Gnode[col.n].fh[col.o]+col.f;  pr1->dfv=cf;

 while(ftell(fp)<Pivot->end)
 {
  fread(&(col),SNOF,1,fp); fread(&(cf),Sdouble,1,fp);
  pr2=Gnode[col.n].fh[col.o]+col.f;
  pr1->dfv+=pr2->dfv*cf;
 }
 head=Pivot;
 Pivot=Pivot->next; 
 free(head);
} while (Pivot);
fclose(fp);
}
/***********************************************************/
void nod_type(ne) int ne;
{
int loc,n,o,f,i,j;
NODE *pn;
ROW  *pr;
NOF nof;

for(loc=0;loc<9;loc++)
{
n=Element[ne].n[loc];
if(Last_e[n]==ne) ntyp[loc]='f';
else              ntyp[loc]='u';
}
}
/**********************************************************/
static void initial()
{
int i,j,n;

for(i=0;i<MAXWD;i++)
{
Fi[i].n=-2; Fj[i].n=-2;
for(j=0;j<MAXWD;j++) Fmat[i][j]=0.;
}
Lasti=0; Lastj=0; Maxwd=0;

Pivot=newpiv(); Pivot->next=NULL; 
for(i=0;i<MAXELE;i++) for(n=0;n<9;n++) Last_e[Element[i].n[n]]=i;
MAXpiv=0.; MINpiv=10000.;
}
/***********************************************************/
void solv(Eq) char Eq;
{
extern void  asm_Newt(),asm_sf();
ELEMENT e;
int ne,n,count;

Fg=fopen("GAUSS","rb");
Fpiv=fopen("piv_dat","wb");
initial(); count=0;
printf("************* solv %c\n",Eq);
for(ne=0;ne<MAXELE;ne++)
{
  e=Element[ne];
  if(Eq=='N') asm_Newt(ne);
  if(Eq=='S') asm_sf(ne);
  nod_type(ne);
  count+=condense(e,Eq);
}
fclose(Fpiv);
fclose(Fg);

printf("Number of freedom solved =%d\n",count);
printf("Lasti=%d   Lastj=%d\n",Lasti,Lastj);
fprintf(Fprss,"Number of freedom solved =%d\n",count);
fprintf(Fprss,"Max. width =%d\n",Maxwd);
fprintf(Fprss,"Max. pivot =%e\n",MAXpiv);
fprintf(Fprss,"Min. pivot =%e\n",MINpiv);

Pivot=Pivot->next;
backsub();
}
/***********************************************************/