fem3.cpp

单元网格划分 一个最基本的有限元计算程序 二维传热问题

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

int Gauss(double a[],double b[],int n);	//全选主元高斯消去法

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];					//单元有限元特征式系数矩阵
};

//--------------  全局变量  ---------------------------
	int	i,j,k,							//循环指标
		id;								//单元的循环指标
	const int nx=21,ny=21;              //x,y方向划分网格数,三角形单元个数=nx*ny*2 
	                                    //  (即24*24*2=1152个单元)
	const double Lx=1,Ly=1;				//矩形区域的两边长
	double	D;							//单元矩阵行列式值
	const int iNode=(nx+1)*(ny+1);		//结点个数=(nx+1)*(ny+1) (即25*25=625个节点)
	double* pMatrix;					//总体矩阵指针
	double* pMf;						//f向量指针
	ElementTriangle* pE;				//单元三角形结构体数组指针
	double ai,bi,ci;					//基函数的系数


//--------------------  主程序  -------------------------//
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
	}
	

	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++)
			{	
				pE[i].Aij[l][m]=( pE[i].b[l]*pE[i].b[m] + 
					pE[i].c[l]*pE[i].c[m] ) * pE[i].A/3;
			}
	}

	printf("OK!\n");
    
	//单元矩阵元素累加到总体矩阵相应的位置上	
	printf("单元矩阵元素累加到总体矩阵相应的位置上...\n");
	int idx=0;
	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 ]=0;
		}
	printf("OK!\n");

///////////////////////////////// respaired by zqy 2003,6,20

	double dBig=pow(10,20);				    //边界条件对角线扩大法处理所用的大数							
	for(i=0;i<nx+1;i++)
	{   double Ur=1.0;                      //强制边界条件1
		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;
	}
	double Ur=0.0;                          //强制边界条件2
    for(i=0;i<nx+1;i++)
	{   pMatrix[i*iNode+i]*=dBig;
	    pMf[i]=0;
	}
	for(j=0;j<nx;j++)
	{	i=(nx+1)*(j+1)-1;
		pMatrix[i*iNode+i]*=dBig;
		pMf[i]=0;
	}
	for(j=0;j<nx;j++)
    {   i=(nx+1)*j;
	    pMatrix[i*iNode+i]*=dBig;
		pMf[i]=0;
	}

//////////////////////////////// respaired by zqy 2003,6,20

	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]);
	}
	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);
}