f2.c

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

/******************************* f2.c ***************************/

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

/*********************************************************/
double Legend(n,x)  int n; double x;
{
if(n<0) { printf("Error in MT(n,x)\n"); exit(1);}
if(n==0) return 1.;
if(n==1) return x ;
return ((2.*n-1)*x*Legend(n-1,x)-(n-1)*Legend(n-2,x))/n; 
}
/*********************************************************/
double MT(n,x) int n; double x;
{
if(n<0) { printf("Error in MT(n,x)\n"); exit(1);}
if(n==0) return((1.-x)/2.);
if(n==1) return((1.+x)/2.);
return (Legend(n,x)-Legend(n-2,x))/sqrt(4.*n-2.);
}
/*********************************************************/
double dMT(n,x)  /*derivative of MT(n,x) resp. to x */
int n;
double x;
{
if(n<0) { printf("Error in MT(n,x)\n"); exit(1);}
if(n==0) return -0.5;
if(n==1) return  0.5;
return   Legend(n-1,x)*sqrt(n-0.5);
}
/*********************************************************/
double L2(n,x) int   n; double x;
{
switch(n) {
case 0: return (1.-x)/2.;
case 1: return (1.+x)/2.;
default: printf(" L2_ERROR!!"); exit(1);
          }
} 
/*********************************************************/
double dL2(n,x) int   n; double x;
{
switch(n) {
case 0: return -0.5;
case 1: return  0.5;
default: printf(" dL2_ERROR!!"); exit(1);
          }
} 
/*********************************************************/
double L4(i, gp) int   i; PAIR gp;
{
double k,a;
k=gp.k;  a=gp.a;
switch(i) {
case 0: return   L2(0,k)*L2(0,a);
case 1: return   L2(1,k)*L2(0,a);
case 2: return   L2(1,k)*L2(1,a);
case 3: return   L2(0,k)*L2(1,a);
default: printf(" L4_ERROR!!"); exit(1); }
} 
/*********************************************************/
double L4k(i, gp) int   i; PAIR gp;
{
double k,a;
k=gp.k;  a=gp.a;
switch(i) {
case 0: return   dL2(0,k)*L2(0,a);
case 1: return   dL2(1,k)*L2(0,a);
case 2: return   dL2(1,k)*L2(1,a);
case 3: return   dL2(0,k)*L2(1,a);
default: printf(" L4_ERROR!!"); exit(1); }
} 
/*********************************************************/
double L4a(i, gp) int   i; PAIR gp;
{
double k,a;
k=gp.k;  a=gp.a;
switch(i) {
case 0: return   L2(0,k)*dL2(0,a);
case 1: return   L2(1,k)*dL2(0,a);
case 2: return   L2(1,k)*dL2(1,a);
case 3: return   L2(0,k)*dL2(1,a);
default: printf(" L4_ERROR!!"); exit(1); }
} 
/*********************************************************/
double L3(n,x) int   n; double x;
{
switch(n) {
case 0: return x*(x-1.)/2.;
case 1: return x*(x+1.)/2.;
case 2: return -(x-1.)*(x+1.);
default: printf(" L3_ERROR!!"); exit(1); }
} 
/*********************************************************/
double dL3(n,x) int   n; double x;
{
switch(n) {
case 0: return (2.*x-1.)/2.;
case 1: return (2.*x+1.)/2.;
case 2: return -2.*x;
default: printf(" dL3_ERROR!!"); exit(1); }
} 
/*********************************************************/
double L9(i, gp) int   i; PAIR gp;
{
double k,a;
k=gp.k;  a=gp.a;
switch(i) {
case 0: return   L3(0,k)*L3(0,a);
case 1: return   L3(1,k)*L3(0,a);
case 2: return   L3(1,k)*L3(1,a);
case 3: return   L3(0,k)*L3(1,a);
case 4: return   L3(2,k)*L3(0,a);
case 5: return   L3(1,k)*L3(2,a);
case 6: return   L3(2,k)*L3(1,a);
case 7: return   L3(0,k)*L3(2,a);
case 8: return   L3(2,k)*L3(2,a);
default: printf(" L9_ERROR!!"); exit(1); }
} 
/*********************************************************/
double L9k(i, gp) int   i; PAIR gp;
{
double k,a;
k=gp.k;  a=gp.a;
switch(i) {
case 0: return   dL3(0,k)*L3(0,a);
case 1: return   dL3(1,k)*L3(0,a);
case 2: return   dL3(1,k)*L3(1,a);
case 3: return   dL3(0,k)*L3(1,a);
case 4: return   dL3(2,k)*L3(0,a);
case 5: return   dL3(1,k)*L3(2,a);
case 6: return   dL3(2,k)*L3(1,a);
case 7: return   dL3(0,k)*L3(2,a);
case 8: return   dL3(2,k)*L3(2,a);
default: printf(" L9k_ERROR!!"); exit(1); }
} 
/*********************************************************/
double L9a(i, gp) int   i; PAIR gp;
{
double k,a;
k=gp.k;  a=gp.a;
switch(i) {
case 0: return   L3(0,k)*dL3(0,a);
case 1: return   L3(1,k)*dL3(0,a);
case 2: return   L3(1,k)*dL3(1,a);
case 3: return   L3(0,k)*dL3(1,a);
case 4: return   L3(2,k)*dL3(0,a);
case 5: return   L3(1,k)*dL3(2,a);
case 6: return   L3(2,k)*dL3(1,a);
case 7: return   L3(0,k)*dL3(2,a);
case 8: return   L3(2,k)*dL3(2,a);
default: printf(" L9a_ERROR!!"); exit(1); }
} 
/*********************************************************/
double shapeL(nod,of,eo,gp) int  nod,of,eo; PAIR    gp;
{
if(eo==1) return L4(nod,gp);
else      return L9(nod,gp);
}
/*********************************************************/
/*********************************************************/
double X(nodx,gp) double nodx[9]; PAIR gp;
{
int i;
double x;
for(x=0.,i=0;i<9;i++) x+=nodx[i]*L9(i,gp);
return(x);
}
/*********************************************************/
double Y(nody,gp) double nody[9]; PAIR gp;
{
int   i;
double y;
for(y=0.,i=0;i<9;i++) y+=nody[i]*L9(i,gp);
return(y);
}
/*********************************************************/
double Xk(nodx,gp) double nodx[9]; PAIR gp;
{
int i;
double xk;
for(xk=0.,i=0;i<9;i++) xk+=nodx[i]*L9k(i,gp);
return(xk);
}
/*********************************************************/
double Xa(nodx,gp) double nodx[9]; PAIR gp;
{
int i;
double xa;
for(xa=0.,i=0;i<9;i++) xa+=nodx[i]*L9a(i,gp);
return(xa);
}
/*********************************************************/
double Yk(nody,gp) double nody[9]; PAIR gp;
{
int i;
double yk;
for(yk=0.,i=0;i<9;i++) yk+=nody[i]*L9k(i,gp);
return(yk);
}
/*********************************************************/
double Ya(nody,gp) double nody[9]; PAIR gp;
{
int i;
double ya;
for(ya=0.,i=0;i<9;i++) ya+=nody[i]*L9a(i,gp);
return(ya);
}
/*********************************************************/
double Jacob(nodx,nody,gp)  /* Jacob=(xk*ya-xa*yk) */
double nodx[9],nody[9];
PAIR gp;
{ return(Xk(nodx,gp)*Ya(nody,gp)-Xa(nodx,gp)*Yk(nody,gp)); }
/*********************************************************/
void MN(nod,f,eo,ma,na) int     nod,f,eo,*ma,*na;
{
int   m,n;

if((nod<4)&&(f>1)) {printf(" nod and f missmatched!!\n");exit(1);}
if((eo==1)&&(nod==8)){printf(" nod and eo missmatched!!\n");exit(1);}
switch(nod) 
{
case 0: m=0; n=0; break;
case 1: m=1; n=0; break;
case 2: m=1; n=1; break;
case 3: m=0; n=1; break;
case 4: m=f+2; n=0; break;
case 5: m=1; n=f+2; break;
case 6: m=f+2; n=1; break;
case 7: m=0; n=f+2; break;
case 8: m=f/(eo-1)+2; n=f%(eo-1)+2; break;
default: printf(" nod ERROR!!\n"); exit(1);
}
*ma=m;  *na=n;
}
/*********************************************************/
double shape(nod,f,eo,gp) int  nod,f,eo; PAIR    gp;
{
int   m,n;
MN(nod,f,eo,&m,&n);
return MT(m,gp.k)*MT(n,gp.a);
}
/*********************************************************/
double dsh_k(nod,f,e,gp) int  nod,f,e; PAIR  gp;
{
int   m,n;
MN(nod,f,e,&m,&n);
return dMT(m,gp.k)*MT(n,gp.a);
}
/*********************************************************/
double dsh_a(nod,f,e,gp) int  nod,f,e; PAIR  gp;
{
int   m,n;
MN(nod,f,e,&m,&n);
return MT(m,gp.k)*dMT(n,gp.a);
}
/*********************************************************/