?? fem2.cpp
字號:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int Gauss(double a[],double b[],int n); //全選主元高斯消去法
/*double GaussIntegral(int n,int js[], //計算多重積分的高斯方法
double(*fn)(int n,double x[]),
void(*fnGaussLimit)(int j, int n,double x[],double y[])); */
struct ETNode{ //單元結點結構體
double x,y; //單元結點坐標
int number; //單元結點在總體區域劃分中的編號
};
struct ElementTriangle{ //三角形單元結構體
ETNode nd[3]; //存儲相對應的總體結點號
double a[3],b[3],c[3]; //基函數的系數
double A; //單元面積
double Aij[3][3]; //單元有限元特征式系數矩陣
double fi[3]; //單元上的積分值
};
//-------------- 全局變量 ---------------------------
int i,j,k, //循環指標
id; //單元的循環指標
const int nx=21,ny=21; //x,y方向劃分網格數,三角形單元個數=nx*ny*2
const double Lx=1,Ly=1; //矩形區域的兩邊長
double D; //單元矩陣行列式值
const int iNode=(nx+1)*(nx+1); //結點個數
double* pMatrix; //總體矩陣指針
double* pMf; //f向量指針
ElementTriangle* pE; //單元三角形結構體數組指針
double ai,bi,ci; //基函數的系數
//-----------------------------------------------------
double fn(int n,double x[2]) //被積函數,高斯積分函數的參數
{
return(ai+bi*x[0]+ci*x[1]);
}
/*
void fnGaussLimit(int jFlag,int n, //積分上下限函數,高斯積分函數的參數
double x[],double y[2])
{ //矩形中第一個三角形單元積分上下限計算
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, //積分上下限函數,高斯積分函數的參數
double x[],double y[2])
{ //矩形中第二個三角形單元積分上下限計算
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;
}
}
*/
//-------------------- 主程序 -------------------------
//有限元理論請參考章本照先生編著的<<流體力學中的有限元方法>>, PP.156-165
//機械工業出版社出版,1986
void main(void)
{
//為總體矩陣,三角形單元數組,f函數向量分配存儲內存
pMatrix=(double*)malloc(iNode*iNode*sizeof(double));
pE=(ElementTriangle*)malloc(nx*ny*2*sizeof(ElementTriangle));
pMf=(double*)malloc(iNode*sizeof(double));
//初始化值為0,因為下面要累加總體矩陣
for(i=0;i<iNode*iNode;i++)
pMatrix[i]=0;
for(i=0;i<iNode;i++)
pMf[i]=0;
try{
//------ 計算得到網格的信息 -----------
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("計算基函數系數值...\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("計算單元有限元特征式系數矩陣...\n");
//--------------------------------------計算單元有限元特征式系數矩陣可以不再分兩個三角形循環
int l,m;
for(i=0;i<nx*ny*2;i++)
{ for(l=0;l<3;l++)
for(m=0;m<3;m++) // Respaired
{
pE[i].Aij[l][m]=( pE[i].b[l]*pE[i].b[m] +
pE[i].c[l]*pE[i].c[m] ) * pE[i].A;
}
}
//------------------------------------------------------------------------------------------
/* for(i=0;i<nx;i++) //計算單元有限元特征式系數矩陣
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("計算積分值,填充到f函數向量數組...\n");
//---------------------------計算三角形單元f函數向量不用多重積分公式,直接代用單元三角形面積的三分之一
int idx=0;
for(i=0;i<2*nx*ny;i++)
{ for(idx=0;idx<3;idx++)
{
pE[i].fi[idx]=(1.0/3.0)*(0.5)*(Lx/nx)*(Ly/ny); // Respaired
}
}
/* printf("pE[0].fi[0]=%f\n",pE[0].fi[0]);
printf("pE[0].fi[1]=%f\n",pE[0].fi[1]);
printf("pE[0].fi[2]=%f\n",pE[0].fi[2]); */
//------------------------------------------------------------------------------------------
/* static int js[2]={4,4}; //每一層積分區間均分為4個子區間
int idx=0;
for(i=0;i<nx;i++) //計算積分值,填充到f函數向量數組
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");
//單元矩陣元素累加到總體矩陣相應的位置上
printf("單元矩陣元素累加到總體矩陣相應的位置上...\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); //邊界條件對角線擴大法處理所用的大數
double Ur=0.0; //邊界條件1(邊界條件2通過Calerkin弱解表達式自動滿足)
for(i=0;i<nx+1;i++)
{ j=nx+1;
pMatrix[(j*nx+i)*iNode+(j*nx+i)]*=dBig;
pMf[(j*nx+i)]=pMatrix[(j*nx+i)*iNode+(j*nx+i)]*Ur;
}
for(j=0;j<nx+1;j++)
{ i=(nx+1)*(j+1)-1;
pMatrix[i*iNode+i]*=dBig;
pMf[i]=pMatrix[i*iNode+i]*Ur;
}
printf("調用全選主元高斯消去法函數解方程組...\n");
Gauss(pMatrix,pMf,iNode); //調用全選主元高斯消去法函數解方程組
printf("OK!\n");
printf("寫計算結果數據到文件...\n");
FILE *wfp; //文件操作
if((wfp=fopen("dat.txt","w+"))==NULL)
printf("Cann't open the file... ");
//fprintf(wfp,"計算得各結點上的溫度值為:\n");
for(i=0;i<iNode;i++)
fprintf(wfp,"%f\n", pMf[i]);
//fprintf(wfp,"%d %f\n",i+1,pMf[i]);
printf("OK!\n");
fclose(wfp);
}
catch(...)
{
printf("Error occured...\n");
}
//釋放總體矩陣和三角形單元數組占用內存
free(pMf); free(pE); free(pMatrix);
printf("Please press Enter to exit...");
getchar();
}
//---------- 全選主元高斯消去法 ------------------------------
// a 體積為n*n的雙精度實型二維數組，方程組系數矩陣，返回時將被破壞
// b 長度為n的雙精度實型一維數組，方程組右端的常數向量，返回方程組的解向量
// n 整型變量，方程組的階數
//--------------------------------------------------------------
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);
}
/*
//------------ 計算多重積分的高斯方法 ------------------------
// n 整型變量，積分重數
// js 整型一維數組，長度為n，
// 其js(i) (i=0,1,...,n-1)表示第i層積分區間所劃分的子區間個數
// fn() 被積函數(函數指針)
// fnGaussLimit() 函數(函數指針)計算各層積分上下限值(上限>下限)
//--------------------------------------------------------------
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;
}
//-------------------------------------------------------------*/
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -