shenliu.cpp

有限元计算基本程序

	{
		printf("%d   %lf    ",i,pMf[i]);
		if(i%4==0)
			printf("\n");
	}


    printf("OK!\n");
	
//////////////////////////////////////////////////////////////////////////////////////////////
	int  cishu=0;
	double max;
	

    for(m=0;m<iNode;m++)                                    
	    zhi[m]=cnode[m].y-pMf[m];

	max=fabs(zhi[0]);
	for(i=1;i<num_waterout;i++)                                      
	{
		if(fabs(zhi[i])>max)
			max=fabs(zhi[i]);
	}
	int out_point=num_waterout;
	int sumbory=num_border;
	


	printf("迭代使用高斯公式计算...\n");
	
  while(max > water_error)
  {     
        cishu++;
		
		for(i=0;i<out_point;i++)
		{
			if(zhi[outpoint[i]]<=-1.0E-10)
			{
			             
			    bl[num_border].jd=outpoint[i];
				bl[num_border].dx=cnode[outpoint[i]].y;
				num_border++;
			}
			
			
		}

        out_point=num_waterout-(num_border-sumbory);


		

	//处理潜润面上下部分单元的渗透系数问题///////
     for(j=0;j<sum_material;j++)
	{
		for(i=place_mat[j][0];i<=place_mat[j][1];i++)
		{
		     if (zhi[i]>=1.0E-10)
			 {
                   mp[j].B[0][0]=mp[j].B[0][0]/1000;
			       mp[j].B[0][1]=mp[j].B[0][1]/1000;
		           mp[j].B[1][0]=mp[j].B[1][0]/1000;
			       mp[j].B[1][1]=mp[j].B[1][1]/1000;
			 }
		}
	}
	
///	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");

////初始化值为0,因为下面要累加总体矩阵
	for(k=0;k<iNode*iNode;k++)
        pMatrix[k]=0;
	for(i=0;i<iNode;i++)
		pMf[i]=0;
	
//单元矩阵元素累加到总体矩阵相应的位置上////////////////////////////////////////////////////	
	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++)
	{
		printf("%d   %lf    ",i,pMf[i]);
		if(i%4==0)
			printf("\n");
	}
	printf("OK\n");
	printf("迭代进行次数：%d\n",cishu);
	


////////////////二次计算完成////////////////////////////////////////////////////
	
	for(m=0;m<iNode;m++)                                    
	    zhi[m]=pMf[m]-cnode[m].y;

	max=fabs(zhi[0]);
	for(i=1;i<num_waterout;i++)                                      
	{
		if(fabs(zhi[i])>max)
			max=fabs(zhi[i]);
	}
		
	}
	
	printf("OK!\n");

	printf("写计算结果数据到文件...\n");
	FILE *wfp;							         //文件操作
	if((wfp=fopen("dat.txt","w+"))==NULL)
		printf("Cann't open the file... ");

 //   fprintf(wfp,"重复计算的次数为:\n");
//	fprintf(wfp,"%d\n",cishu);
	fprintf(wfp,"计算得各结点上的水头值为:\n");
	for(i=0;i<iNode;i++)                        //-cishu*num_waterout
	{
		fprintf(wfp,"%d   %f\n",i,pMf[i]);
	}
	printf("OK!\n");

	fclose(wfp);

}

catch(...)
{	
	printf("Error occured...\n");
}
	//释放总体矩阵和三角形单元数组占用内存
	delete [] pMf;	delete []  pMatrix;	delete []  chf;  delete []  bl;	
	delete []  mp;
	delete [] pE;  delete [] cnode;

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