?? fem1.cpp
字號(hào):
// Finite Element Method solve the Poisson formula
// If this program works, it was written by hujinshan@2003.5.9
// If not, I don't know who wrote it :)
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int Gauss(double a[],double b[],int n); //全選主元高斯消去法
double GaussIntegral(int n,int js[], //計(jì)算多重積分的高斯方法
double(*fn)(int n,double x[]),
void(*fnGaussLimit)(int j, int n,double x[],double y[]));
struct ETNode{ //單元結(jié)點(diǎn)結(jié)構(gòu)體
double x,y; //單元結(jié)點(diǎn)坐標(biāo)
int number; //單元結(jié)點(diǎn)在總體區(qū)域劃分中的編號(hào)
};
struct ElementTriangle{ //三角形單元結(jié)構(gòu)體
ETNode nd[3]; //存儲(chǔ)相對(duì)應(yīng)的總體結(jié)點(diǎn)號(hào)
double a[3],b[3],c[3]; //基函數(shù)的系數(shù)
double A; //單元面積
double Aij[3][3]; //單元有限元特征式系數(shù)矩陣
double fi[3]; //單元上的積分值
};
//-------------- 全局變量 ---------------------------
int i,j,k, //循環(huán)指標(biāo)
id; //單元的循環(huán)指標(biāo)
const int nx=21,ny=21; //x,y方向劃分網(wǎng)格數(shù),三角形單元個(gè)數(shù)=nx*ny*2
const double Lx=1,Ly=1; //矩形區(qū)域的兩邊長(zhǎng)
double D; //單元矩陣行列式值
const int iNode=(nx+1)*(ny+1); //結(jié)點(diǎn)個(gè)數(shù)
double* pMatrix; //總體矩陣指針
double* pMf; //f向量指針
ElementTriangle* pE; //單元三角形結(jié)構(gòu)體數(shù)組指針
double ai,bi,ci; //基函數(shù)的系數(shù)
//-----------------------------------------------------
double fn(int n,double x[2]) //被積函數(shù),高斯積分函數(shù)的參數(shù)
{
return(ai+bi*x[0]+ci*x[1]);
}
void fnGaussLimit(int jFlag,int n, //積分上下限函數(shù),高斯積分函數(shù)的參數(shù)
double x[],double y[2])
{ //矩形中第一個(gè)三角形單元積分上下限計(jì)算
switch(jFlag)
{
case 0:{
y[0]=(Lx/nx)*i;
y[1]=(Lx/nx)*(i+1);
break;
}
case 1:{
y[0]=(Ly/ny)*j+ x[0]-(Lx/nx)*i;
y[1]=(Ly/ny)*(j+1);
break;
}
default:
break;
}
}
void fnGaussLimit2(int jFlag,int n, //積分上下限函數(shù),高斯積分函數(shù)的參數(shù)
double x[],double y[2])
{ //矩形中第二個(gè)三角形單元積分上下限計(jì)算
switch(jFlag)
{
case 0:{
y[0]=(Lx/nx)*i;
y[1]=(Lx/nx)*(i+1);
break;
}
case 1:{
y[0]=(Ly/ny)*j;
y[1]=(Ly/ny)*j+ x[0]-(Lx/nx)*i;
break;
}
default:
break;
}
}
//-------------------- 主程序 -------------------------
//有限元理論請(qǐng)參考章本照先生編著的<<流體力學(xué)中的有限元方法>>, PP.156-165
//機(jī)械工業(yè)出版社出版,1986
void main(void)
{
//為總體矩陣,三角形單元數(shù)組,f函數(shù)向量分配存儲(chǔ)內(nèi)存
pMatrix=(double*)malloc(iNode*iNode*sizeof(double));
pE=(ElementTriangle*)malloc(nx*ny*2*sizeof(ElementTriangle));
pMf=(double*)malloc(iNode*sizeof(double));
//初始化值為0,因?yàn)橄旅嬉奂涌傮w矩陣
for(i=0;i<iNode*iNode;i++)
pMatrix[i]=0;
for(i=0;i<iNode;i++)
pMf[i]=0;
try{
//------ 計(jì)算得到網(wǎng)格的信息 -----------
for(j=0;j<nx;j++)
for(i=0;i<ny;i++)
{
//for the first triangle in the rectangle
pE[i*2+j*ny*2].nd[0].x=(Lx/nx)*i;
pE[i*2+j*ny*2].nd[0].y=(Ly/ny)*j;
pE[i*2+j*ny*2].nd[0].number=i+j*(nx+1); //NO.0
pE[i*2+j*ny*2].nd[1].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2].nd[1].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2].nd[1].number=i+1+(nx+1)*(j+1); //NO.1
pE[i*2+j*ny*2].nd[2].x=(Lx/nx)*i;
pE[i*2+j*ny*2].nd[2].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2].nd[2].number=i+(nx+1)*(j+1); //NO.2
//for the second triangle in the rectangle
pE[i*2+j*ny*2+1].nd[0].x=(Lx/nx)*i;
pE[i*2+j*ny*2+1].nd[0].y=(Ly/ny)*j;
pE[i*2+j*ny*2+1].nd[0].number=i+j*(nx+1); //NO.0
pE[i*2+j*ny*2+1].nd[1].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2+1].nd[1].y=(Ly/ny)*j;
pE[i*2+j*ny*2+1].nd[1].number=i+j*(nx+1)+1; //NO.1
pE[i*2+j*ny*2+1].nd[2].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2+1].nd[2].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2+1].nd[2].number=i+1+(nx+1)*(j+1); //NO.2
}
//---------------------------------------------------------
//please turn to page 158 for more details
printf("計(jì)算基函數(shù)系數(shù)值...\n");
for(id=0;id<nx*ny*2;id++)
{
for(i=0;i<3;i++)
{
if(i==0) j=1,k=2;
else if(i==1) j=2,k=0;
else if(i==2) j=0,k=1;
pE[id].A=( (pE[id].nd[j].x-pE[id].nd[i].x)*(pE[id].nd[k].y-pE[id].nd[i].y)-
(pE[id].nd[j].y-pE[id].nd[i].y)*(pE[id].nd[k].x-pE[id].nd[i].x) )/2.0;
D=2.0*pE[id].A;
pE[id].a[i]=( pE[id].nd[j].x*pE[id].nd[k].y- pE[id].nd[k].x*pE[id].nd[j].y )/D;
pE[id].b[i]=( pE[id].nd[j].y-pE[id].nd[k].y )/D;
pE[id].c[i]=( pE[id].nd[k].x-pE[id].nd[j].x )/D;
}
}printf("OK!\n");
printf("計(jì)算單元有限元特征式系數(shù)矩陣...\n");
int l,m;
for(i=0;i<nx;i++) //計(jì)算單元有限元特征式系數(shù)矩陣
for(j=0;j<ny;j++)
{
for(l=0;l<3;l++) //for the first triangle in the rectangle
for(m=0;m<3;m++)
{
pE[i*2+j*ny*2].Aij[l][m]=( pE[i*2+j*ny*2].b[l]*pE[i*2+j*ny*2].b[m] +
pE[i*2+j*ny*2].c[l]*pE[i*2+j*ny*2].c[m] ) * pE[i*2+j*ny*2].A;
}
for(l=0;l<3;l++) //for the second triangle in the rectangle
for(m=0;m<3;m++)
{
pE[i*2+j*ny*2+1].Aij[l][m]=( pE[i*2+j*ny*2+1].b[l]*pE[i*2+j*ny*2+1].b[m] +
pE[i*2+j*ny*2+1].c[l]*pE[i*2+j*ny*2+1].c[m] ) * pE[i*2+j*ny*2+1].A;
}
}printf("OK!\n");
printf("計(jì)算積分值,填充到f函數(shù)向量數(shù)組...\n");
static int js[2]={4,4}; //每一層積分區(qū)間均分為4個(gè)子區(qū)間
int idx=0;
for(i=0;i<nx;i++) //計(jì)算積分值,填充到f函數(shù)向量數(shù)組
for(j=0;j<ny;j++)
{
for(idx=0;idx<3;idx++) //for the first triangle in the rectangle
{
ai=pE[i*2+j*ny*2].a[idx];
bi=pE[i*2+j*ny*2].b[idx];
ci=pE[i*2+j*ny*2].c[idx];
pE[i*2+j*ny*2].fi[idx]=GaussIntegral(2,js,fn,fnGaussLimit);
}
for(idx=0;idx<3;idx++) //for the second triangle in the rectangle
{
ai=pE[i*2+j*ny*2+1].a[idx];
bi=pE[i*2+j*ny*2+1].b[idx];
ci=pE[i*2+j*ny*2+1].c[idx];
pE[i*2+j*ny*2+1].fi[idx]=GaussIntegral(2,js,fn,fnGaussLimit2);
}
}printf("OK!\n");
//單元矩陣元素累加到總體矩陣相應(yīng)的位置上
printf("單元矩陣元素累加到總體矩陣相應(yīng)的位置上...\n");
for(idx=0;idx<nx*ny*2;idx++)
for(i=0;i<3;i++)
{
for(j=0;j<3;j++)
pMatrix[pE[idx].nd[i].number*iNode+pE[idx ].nd[j].number] +=pE[idx].Aij[i][j];
pMf[ pE[idx].nd[i].number ]+=pE[idx].fi[i];
}
printf("OK!\n");
double dBig=pow(10,20); //邊界條件對(duì)角線擴(kuò)大法處理所用的大數(shù)
double Ur=1.0; //邊界條件1(邊界條件2通過Galerkin弱解表達(dá)式自動(dòng)滿足)
for(i=0;i<nx+1;i++)
{ j=nx+1;
pMatrix[(j*nx+i)*iNode+(j*nx+i)]*=dBig;
pMf[(j*nx+i)]*=dBig*Ur;
}
for(j=0;j<nx+1;j++)
{ i=(nx+1)*(j+1)-1;
pMatrix[i*iNode+i]*=dBig;
pMf[i]*=dBig*Ur;
}
printf("調(diào)用全選主元高斯消去法函數(shù)解方程組...\n");
Gauss(pMatrix,pMf,iNode); //調(diào)用全選主元高斯消去法函數(shù)解方程組
printf("OK!\n");
printf("寫計(jì)算結(jié)果數(shù)據(jù)到文件...\n");
FILE *wfp; //文件操作
if((wfp=fopen("dat.txt","w+"))==NULL)
printf("Cann't open the file... ");
//fprintf(wfp,"計(jì)算得各結(jié)點(diǎn)上的溫度值為:\n");
//for(i=0;i<iNode;i++)fprintf(wfp,"%d %f\n",i,pMf[i]);
//--------輸出為Matlab數(shù)組格式------------------
for(i=0;i<iNode;i++)fprintf(wfp,"%f\n",pMf[i]);
/* fprintf(wfp,"[");
for(i=0;i<nx+1;i++)
{
for(j=0;j<nx+1;j++)
fprintf(wfp,"%f ",pMf[i]);
fprintf(wfp,";\n");
}
fprintf(wfp,"]");*/
//-------------------------------------------
printf("OK!\n");
fclose(wfp);
}
catch(...)
{
printf("Error occured...\n");
}
//釋放總體矩陣和三角形單元數(shù)組占用內(nèi)存
free(pMf); free(pE); free(pMatrix);
printf("Please press Enter to exit...");
getchar();
}
//---------- 全選主元高斯消去法 ------------------------------
// a 體積為n*n的雙精度實(shí)型二維數(shù)組，方程組系數(shù)矩陣，返回時(shí)將被破壞
// b 長(zhǎng)度為n的雙精度實(shí)型一維數(shù)組，方程組右端的常數(shù)向量，返回方程組的解向量
// n 整型變量，方程組的階數(shù)
//--------------------------------------------------------------
int Gauss(double a[],double b[],int n)
{
int *js,l,k,i,j,is,p,q;
double d,t;
js=(int*)malloc(n*sizeof(int));
l=1;
for(k=0;k<=n-2;k++)
{
d=0.0;
for(i=k;i<=n-1;i++)
for(j=k;j<=n-1;j++)
{
t=fabs(a[i*n+j]);
if(t>d) { d=t; js[k]=j; is=i;}
}
if(d+1.0==1.0) l=0;
else
{
if(js[k]!=k)
for(i=0;i<=n-1;i++)
{
p=i*n+k; q=i*n+js[k];
t=a[p]; a[p]=a[q]; a[q]=t;
}
if(is!=k)
{
for(j=k;j<=n-1;j++)
{
p=k*n+j; q=is*n+j;
t=a[p]; a[p]=a[q]; a[q]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
}
if(l==0)
{
free(js);
printf("Gauss funtion failed 1...\n");
return(0);
}
d=a[k*n+k];
for(j=k+1;j<=n-1;j++)
{
p=k*n+j; a[p]=a[p]/d;
}
b[k]=b[k]/d;
for(i=k+1;i<=n-1;i++)
{
for(j=k+1;j<=n-1;j++)
{
p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if(fabs(d)+1.0==1.0)
{
free(js);
printf("Gauss funtion failed 2...\n");
return(0);
}
b[n-1]=b[n-1]/d;
for(i=n-2;i>=0;i--)
{
t=0.0;
for(j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for(k=n-1;k>=0;k--)
if(js[k]!=k)
{
t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;
}
free(js);
return(1);
}
//------------ 計(jì)算多重積分的高斯方法 ------------------------
// n 整型變量，積分重?cái)?shù)
// js 整型一維數(shù)組，長(zhǎng)度為n，
// 其js(i) (i=0,1,...,n-1)表示第i層積分區(qū)間所劃分的子區(qū)間個(gè)數(shù)
// fn() 被積函數(shù)(函數(shù)指針)
// fnGaussLimit() 函數(shù)(函數(shù)指針)計(jì)算各層積分上下限值(上限>下限)
//--------------------------------------------------------------
double GaussIntegral(int n,int js[],double (*fn)(int n,double x[]),
void (*fnGaussLimit)(int j,int n,double x[],double y[]))
{
int m,j,k,q,*is;
double y[2],p,s,*x,*a,*b;
static double t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
static double c[5]={0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851};
is=(int*)malloc(2*(n+1)*sizeof(int));
x=(double*)malloc(n*sizeof(double));
a=(double*)malloc(2*(n+1)*sizeof(double));
b=(double*)malloc((n+1)*sizeof(double));
m=1;a[n]=1.0; a[2*n+1]=1.0;
while(true)
{
for(j=m;j<=n;j++)
{
fnGaussLimit(j-1,n,x,y);
a[j-1]=0.5*(y[1]-y[0])/js[j-1];
b[j-1]=a[j-1]+y[0];
x[j-1]=a[j-1]*t[0]+b[j-1];
a[n+j]=0.0; is[j-1]=1; is[n+j]=1;
}
j=n; q=1;
while(q==1)
{
k=is[j-1];
if(j==n) p=fn(n,x);
else p=1.0;
a[n+j]=a[n+j+1]*a[j]*p*c[k-1]+a[n+j];
is[j-1]=is[j-1]+1;
if(is[j-1]>5)
{
if(is[n+j]>=js[j-1])
{
j=j-1; q=1;
if(j==0)
{
s=a[n+1]*a[0];
free(is); free(x);free(a); free(b);
return(s);
}
}
else
{
is[n+j]=is[n+j]+1;
b[j-1]=b[j-1]+a[j-1]*2.0;
is[j-1]=1; k=is[j-1];
x[j-1]=a[j-1]*t[k-1]+b[j-1];
if(j==n)
q=1;
else
q=0;
}
}
else
{
k=is[j-1];
x[j-1]=a[j-1]*t[k-1]+b[j-1];
if(j==n)
q=1;
else
q=0;
}
}
m=j+1;
}
return 0;
}
//--------------------------------------------------------------
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -