📄 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 + -