shenliu.cpp

有限元计算基本程序

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream.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];					//单元有限元特征式系数矩阵
};

struct  Modulus_penetr{                //地层的渗透系数
	int  material;                     //地层编号
	char name_mat[15];                   //地层名称
	double B[2][2];                    //渗透系数
};

struct Char_infile{                    //输入文件中的注释字符
	char title[30];
};

struct Border{                         //边界条件
	int  jd;                           //节点编号
	double dx;                         //水头大小
};

//struct ETNoded{							//迭代增加的单元结点结构体
//	double x,y;							//单元结点坐标
//	int number;							//单元结点在总体区域划分中的编号
//};

struct Border_ele{                       //边界单元结构体
	int a,b,c,d;
};


 

//--------------  全局变量  ---------------------------
	int	i,j,k,l,m,							//循环指标
		id;	                                //单元的循环指标
//    int iN[13];                          //注意改变数组的大小(节点总数)
	int Ele[19];                         //注意改变数组的大小(单元总数)
	double zhi[16];                       //注意改变数组的大小(节点总数)
	double  tem[16];                       //注意改变数组的大小(节点总数)
//	int  tem1[4];                          //注意改变数组的大小(可能的出水节点总数)
//	int  sum_border_ele;                 //边界单元的各数
//	int  jinshui_point;                 //进水点节点编号

	int num_border;                     //已知水头节点的总数
	int num_waterout;                   //可能出水点的节点总数
	int outpoint[200];                  //可能出水节点
	int sum_material;                   //地层渗透系数种类
	int place_mat[40][3];               //不同地层分布位置
	double	D;							//单元矩阵行列式值
	int iNode,element;                  //节点数和单元数
	
//	ElementTriangle* pE;				//单元三角形结构体数组指针

	double ai,bi,ci;					//基函数的系数
	double water_error;                 //自由水面误差结束标准
	



//--------------------  主程序  -------------------------//
void main()
{
    struct Char_infile *chf=new Char_infile[13];                     //输入文件注释结构数组指针及分配内存
	struct Border *bl=new Border[400];                                //边界条件，结构数组指针及分配内存
	FILE *cfptr;
    if ((cfptr=fopen("shu.dat","r"))==NULL)
		printf("无法打开文件.\n");
	else
		fscanf(cfptr,"%s",chf[0].title);
	    fscanf(cfptr,"%s",chf[1].title);
		fscanf(cfptr,"%d",&iNode);
   		fscanf(cfptr,"%s",chf[2].title);
		fscanf(cfptr,"%d",&element);

     //为总体矩阵,三角形单元数组,f函数向量等分配存储内存

//	pE=(ElementTriangle*)malloc(element*sizeof(ElementTriangle));
//	cnode=(ETNode*)malloc(iNode*sizeof(ETNode));
	double *pMatrix=new double[iNode*iNode*sizeof(double)];				//总体矩阵指针
	double *pMf=new double[iNode*sizeof(double)];						//f向量指针
	struct ElementTriangle *pE=new ElementTriangle[element];            //单元
	struct ETNode *cnode=new ETNode[iNode];                     //存坐标

	////"50"为认为规定迭代次数，要求计算次数较多的，根据情况自行改变此值//

	//初始化值为0,因为下面要累加总体矩阵
	for(i=0;i<iNode*iNode;i++)
		pMatrix[i]=0;
	for(i=0;i<iNode;i++)
		pMf[i]=0;

		fscanf(cfptr,"%s",chf[3].title);
		fscanf(cfptr,"%lf",&water_error);
		fscanf(cfptr,"%s",chf[4].title);
        for(i=0;i<element;i++)
			fscanf(cfptr,"%d    %d     %d    %d",&Ele[i],&pE[i].nd[0].number,&pE[i].nd[1].number,&pE[i].nd[2].number);
		fscanf(cfptr,"%s",chf[5].title);
	    for(j=0;j<iNode;j++)
			fscanf(cfptr,"%d    %lf    %lf",&cnode[j].number,&cnode[j].x,&cnode[j].y);
		 for(i=0;i<element;i++)
		{  
			pE[i].nd[0].x=cnode[pE[i].nd[0].number].x;
			pE[i].nd[0].y=cnode[pE[i].nd[0].number].y;
			pE[i].nd[1].x=cnode[pE[i].nd[1].number].x;
		    pE[i].nd[1].y=cnode[pE[i].nd[1].number].y;
			pE[i].nd[2].x=cnode[pE[i].nd[2].number].x;
			pE[i].nd[2].y=cnode[pE[i].nd[2].number].y;
		}
		fscanf(cfptr,"%s",chf[6].title);
        fscanf(cfptr,"%d",&num_border);
		fscanf(cfptr,"%s",chf[7].title);
		for(k=0;k<num_border;k++)                      
			fscanf(cfptr,"%d    %lf",&bl[k].jd,&bl[k].dx);
        fscanf(cfptr,"%s",chf[8].title);
        fscanf(cfptr,"%d",&num_waterout);
		fscanf(cfptr,"%s",chf[9].title);
		for(m=0;m<num_waterout;m++)                           
			fscanf(cfptr,"%d",&outpoint[m]);
		fscanf(cfptr,"%s",chf[10].title);
		fscanf(cfptr,"%d",&sum_material);
		fscanf(cfptr,"%s",chf[11].title);

    struct  Modulus_penetr *mp=new Modulus_penetr[sum_material];     //不同地层的渗透系数，结构数组指针及分配内存

		for(id=0;id<sum_material;id++)
		{
			fscanf(cfptr,"%d    %s",&mp[id].material,mp[id].name_mat);
			fscanf(cfptr,"%lf   %lf",&mp[id].B[0][0],&mp[id].B[0][1]);
			fscanf(cfptr,"%lf   %lf",&mp[id].B[1][0],&mp[id].B[1][1]);
		}
		fscanf(cfptr,"%s",chf[12].title);
		for(l=0;l<sum_material;l++)
			fscanf(cfptr,"%d    %d    %d",&place_mat[l][0],&place_mat[l][1],&place_mat[l][2]);

///	输出读入的文件数据
		printf("%s\n",chf[0].title);
	    printf("%s\n",chf[1].title);
		printf("%d\n",iNode);
		printf("%s\n",chf[2].title);
		printf("%d\n",element);
		printf("%s\n",chf[3].title);
		printf("%E\n",water_error);
		printf("%s\n",chf[4].title);
		for(i=0;i<element;i++)
			printf("%d    %d     %d    %d\n",Ele[i],pE[i].nd[0].number,pE[i].nd[1].number,pE[i].nd[2].number);
		printf("%s\n",chf[5].title);
		for(j=0;j<iNode;j++)
			printf("%d    %lf    %lf\n",cnode[j].number,cnode[j].x,cnode[j].y);
		printf("%s\n",chf[6].title);
        printf("%d\n",num_border);
	    printf("%s\n",chf[7].title);
		for(k=0;k<num_border;k++)                              
	         printf("%d    %lf\n",bl[k].jd,bl[k].dx);
		printf("%s\n",chf[8].title);
        printf("%d\n",num_waterout);
		printf("%s\n",chf[9].title);
		for(m=0;m<num_waterout;m++)                              
			printf("%d\n",outpoint[m]);
		printf("%s\n",chf[10].title);
		printf("%d\n",sum_material);
		printf("%s\n",chf[11].title);
		for(id=0;id<sum_material;id++)
		{
			printf("%d    %s\n",mp[id].material,mp[id].name_mat);
			printf("%lf   %lf\n",mp[id].B[0][0],mp[id].B[0][1]);
			printf("%lf   %lf\n",mp[id].B[1][0],mp[id].B[1][1]);
		}
		printf("%s\n",chf[12].title);
		for(l=0;l<sum_material;l++)
			printf("%d    %d    %d\n",place_mat[l][0],place_mat[l][1],place_mat[l][2]);
	
		printf("节点信息.\n");
		for(i=0;i<element;i++)
		{  
			printf("%d   %lf    %lf\n",pE[i].nd[0].number,pE[i].nd[0].x,pE[i].nd[0].y);
		
			printf("%d   %lf    %lf\n",pE[i].nd[1].number,pE[i].nd[1].x,pE[i].nd[1].y);
			
			printf("%d   %lf    %lf\n\n",pE[i].nd[2].number,pE[i].nd[2].x,pE[i].nd[2].y);
		
		}


try{
	
	printf("计算基函数系数值...\n");
	for(id=0;id<element;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");

    for(j=0;j<sum_material;j++)
	{
		for(i=place_mat[j][0];i<=place_mat[j][1];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]*mp[j].B[0][0] + 
					pE[i].c[l]*pE[i].c[m]*mp[j].B[1][1])*pE[i].A;
				}
		}
	}

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

    ///////////边界条件的处理////////////
    printf("边界条件的处理...\n");

	double dBig=pow(10,25);				    //边界条件对角线扩大法处理所用的大数
	for(k=0;k<num_border;k++)                          
	{   
		pMatrix[bl[k].jd*iNode+bl[k].jd] *= dBig;
	    pMf[bl[k].jd]=pMatrix[bl[k].jd*iNode+bl[k].jd]*bl[k].dx;
	}
	
	printf("OK!\n");


/////////////////////////////////////////////////////////////////////////////////

	printf("调用全选主元高斯消去法函数解方程组...\n");
	Gauss(pMatrix,pMf,iNode);			         //调用全选主元高斯消去法函数解方程组
	for(i=0;i<iNode;i++)