finiteem.cpp

求电磁场边值

#include<iostream.h>
#include<iomanip.h>
#include<fstream.h>
#include"FiniteEM.h"
#include"Vector.h"

const double Width=0.2;
const double High=0.2;
const int numOfPoint=41;
const int numOfArea=58;  //int numOfPoint=((width/2)/h+1)*((high/2)/h+1+1);

Point pt;
double triArea(int);

FiniteEM::FiniteEM()
{
	int i,j;
	pArray=new Array(3,numOfArea);	

	//构造边界联系数组
	pBonArr2D=new Array(9,2);
	pBoundArr2D=pBonArr2D->getArr2D();
	for(i=0;i<4;i++)
		for(j=0;j<2;j++)
		{
			int pnC=pt.getPointCode(0.02*i+0.02*j,0.1,0.02);
			pBoundArr2D[i][j]=pnC;
		}
	for(i=4;i<9;i++)
		for(j=0;j<2;j++)
		{
			int pnC=pt.getPointCode(0.02*j,0.1,0.02);
			pBoundArr2D[i][j]=pnC;
		}
}

void FiniteEM::areaDisper(double h)        //h:单位长度
{                   
	double wid=Width/2;
	double hig=High/2;
	point=new Point[numOfPoint];
	//构造点并初始化
	int i=0;
	int j,numTri,pnC,tag;  //pnC:全局编码  tag:区别点的类型 (0:内部点 1:0V的点 2:10V的点 3:自然边界的点)	
	for(;i<5;i++)
	{
		for(j=0;j<6;j++)
		{
			pnC=(int)(wid/0.02+1)*i+j+1;
			int n=pnC-1;
			if(i==0||j==0) tag=1;             //0V的点
			else if((i!=0)&&(i!=5)&&(j==5)) tag=3;             //自然边界的点
			else tag=0;             //内部点
			point[n].setPoint(double(j*h),double(i*h),pnC,tag);
		}
	}
	for(j=0;j<6;j++)
	{
		i=5;
		pnC=(int)(wid/0.02+1)*i+j+1;
		int n=pnC-1;
		if(j==0) tag=1;
		else if((i==5&&j==4)||(i==5&&j==5)) tag=2;
		else if((i==5&&j==0)) tag=3;
		else tag=0;
		point[n].setPoint(double(j*h),double(4*h+(i-4)*0.01),pnC,tag);
	}
	for(j=0;j<5;j++)
	{
		i=6;
		pnC=(int)(wid/0.02+1)*i+j+1;              
		int n=pnC-1;
		if(j==0) tag=1;
		else if(i==6&&j==4)  tag=2;
		else  tag=3;
		point[n].setPoint(double(j*h),double(4*h+(i-4)*0.01),pnC,tag);
	}

	pArr2D=pArray->getArr2D();
	//设置联系数组第一行
	for(numTri=1;numTri<=40;numTri++)
	{
		int tg=1;
		if(numTri%10==0)
		{	
			pnC=pt.getPointCode(0.1,(numTri/10)*0.02,tg);		
			pArr2D[1-1][numTri-1]=pnC;
		}
		else if(numTri%10!=0)
		{
			pnC=pt.getPointCode((numTri%10/2)*0.02,(numTri/10+1)*0.02,tg); 
			pArr2D[1-1][numTri-1]=pnC;
		}
	}
	for(numTri=41;numTri<=58;numTri++)
	{
		int tg=2;
		if(numTri%10==0)
		{
			pnC=pt.getPointCode(0.1,(numTri/10)*0.02-0.01,tg);
			pArr2D[1-1][numTri-1]=pnC;
		}
		else if(numTri%10!=0)
		{
			pnC=pt.getPointCode((numTri%10/2)*0.02,(4*0.02)+((numTri-40)/10+1)*0.01,tg); 
			pArr2D[1-1][numTri-1]=pnC;
		}
	}	
	//设置第二行
	for(numTri=1;numTri<=40;numTri++)
	{
		int tg=1;
		if(numTri%2!=0)          //设置奇数元的局部编码
		{
			pnC=pt.getPointCode((numTri%10/2)*0.02,(numTri/10)*0.02,tg);
			pArr2D[2-1][numTri-1]=pnC;
		}
		else if(numTri%2==0&&numTri%10!=0)          //设置偶元的局部编码
		{
			pnC=pt.getPointCode((numTri%10/2-1)*0.02,(numTri/10+1)*0.02,tg); 
			pArr2D[2-1][numTri-1]=pnC;
		}
		else if(numTri%2==0&&numTri%10==0)
		{
			pnC=pt.getPointCode(0.1-0.02,(numTri/10)*0.02,tg);
			pArr2D[2-1][numTri-1]=pnC;
		}
	}
	for(numTri=41;numTri<=58;numTri++)
	{
		int tg=2;
		if(numTri%2!=0)          //设置奇数元的局部编码
		{
			pnC=pt.getPointCode((numTri%10/2)*0.02,4*0.02+((numTri-40)/10)*0.01,tg);
			pArr2D[2-1][numTri-1]=pnC;
		}
		else if(numTri%2==0&&numTri%10!=0)          //设置偶元的局部编码
		{
			pnC=pt.getPointCode((numTri%10/2-1)*0.02,4*0.02+((numTri-40)/10+1)*0.01,tg); 
			pArr2D[2-1][numTri-1]=pnC;
		}
		else if(numTri%2==0&&numTri%10==0)
		{
			pnC=pt.getPointCode(0.1-0.02,4*0.02+((numTri-40)/10)*0.01,tg);
			pArr2D[2-1][numTri-1]=pnC;
		}
	}
	//设置第三行
	for(numTri=1;numTri<=40;numTri++)
	{
		int tg=1;
		if(numTri%2!=0)
		{
			pnC=pt.getPointCode((numTri%10/2)*0.02+0.02,(numTri/10)*0.02,tg);
			pArr2D[3-1][numTri-1]=pnC;
			pArr2D[3-1][numTri+1-1]=pnC;
		}
		else continue;
	}
	for(numTri=41;numTri<=58;numTri++)
	{
		int tg=2;
		if(numTri%2!=0)
		{
			pnC=pt.getPointCode((numTri%10/2)*0.02+0.02,4*0.02+((numTri-40)/10)*0.01,tg);
			pArr2D[3-1][numTri-1]=pnC;
			pArr2D[3-1][numTri+1-1]=pnC;
		}
		else continue;
	}
}

//单元插值
void FiniteEM::insertValueOfArea()
{
//	pA=A_D2.setDoubArr2D(numOfArea,3);
	pB=B_D2.setDoubArr2D(numOfArea,3);      //存b1 b2 b3
	pC=C_D2.setDoubArr2D(numOfArea,3);      //存c1 c2 c3
	for(int i=0;i<58;i++)		
	{			
		int index_1=pArr2D[0][i];	//全局编码		
		int index_2=pArr2D[1][i];			
		int index_3=pArr2D[2][i];		
		double y_1=point[index_1-1].getYcoord();			
		double y_2=point[index_2-1].getYcoord();			
		double y_3=point[index_3-1].getYcoord();			
		double x_1=point[index_1-1].getXcoord();			
		double x_2=point[index_2-1].getXcoord();			
		double x_3=point[index_3-1].getXcoord();					
		pB[i][0]=y_2-y_3;					      //b1
		pB[i][1]=y_3-y_1;					//b2
		pB[i][2]=y_1-y_2;			
		pC[i][0]=x_3-x_2;					
		pC[i][1]=x_1-x_3;				
		pC[i][2]=x_2-x_1;	
	}
}

//构造单元K矩阵
void FiniteEM::consUnitMatrix()
{
	pMatrix_D3=new double[58][3][3];
	for(int n=0;n<58;n++)
	{
		for(int i=0;i<3;i++)
		{
			for(int j=0;j<3;j++)
			{
				pMatrix_D3[n][i][j]=(pB[n][i]*pB[n][j]+pC[n][i]*pC[n][j])/(4*triArea(n));
			}
		}
	}
}

//由单元矩阵生成整体矩阵
double** FiniteEM::consK_Matrix()
{
	pK_matrix=K_matrix.setDoubArr2D(numOfPoint,numOfPoint);
	for(int n=0;n<58;n++)
	{
		for(int i=0;i<3;i++)
		{
			for(int j=0;j<3;j++)
			{
				int pCi=pArr2D[i][n]-1;
				int pCj=pArr2D[j][n]-1;
				pK_matrix[pCi][pCj]=pK_matrix[pCi][pCj]+pMatrix_D3[n][i][j];
			}
		}
	}
	return pK_matrix;
			
}

//加入自然边界条件
Vector FiniteEM::addDirBonValue()          //1:0V     2:10V
{
	
	Vector b(41,0);
	for(int i=0;i<41;i++)
	{
		if(point[i].getTag()==1)
		{
			b[i]=0;                             //电位值为零 i=j-1
			for(int j=0;j<41;j++)
			{
				if(j!=i)  b[j]=b[j]-pK_matrix[j][i]*0;
			}

			for(int r=0;r<41;r++)
			{
				if(r!=i) pK_matrix[r][i]=0;
			}
			for(int c=0;c<41;c++)
			{
				if(c!=i) pK_matrix[i][c]=0;
			}
			pK_matrix[i][i]=1;
		}
		if(point[i].getTag()==2)
		{
			b[i]=10;
			for(int j=0;j<41;j++)
			{
				if(j!=i)  b[j]=b[j]-pK_matrix[j][i]*10;
			}

			for(int r=0;r<41;r++) pK_matrix[r][i]=0;
			for(int c=0;c<41;c++) pK_matrix[i][c]=0;
			pK_matrix[i][i]=1;
		}	
	}
	return b;
}

//矩阵压缩
Vector FiniteEM::condense_equation(Vector& b)
{
	int n=0;
	for(int i=0;i<41;i++)		
	{		
		if(point[i].getTag()==1||point[i].getTag()==2)		
		{
			K_matrix.delRow(i-n);
			K_matrix.delCol(i-n);
			b.del_elemt(i-n);
			n++;
		}
	}
	return b;
}

double FiniteEM::CG(Vector& x,Vector& b,double eps,int& iter)
{
	if (K_matrix.getRow2D()!=b.size())
		cout<<"matrix and vector sizes do not match"<<endl;
	int i=0;
	const int max_iter = iter;
	Vector r = b - K_matrix*x;   //初始残差
	Vector p = r;              //p:搜索方向
	Vector v_d;
	double zr = v_d.dot(r,r);   //r和r的内积
	const double stp = eps*b.max_norm(); //制定停止运算的标准
	if (r.max_norm()==0)      //如果初始估计是真实值
	{
		eps = 0.0;             //返回，否则发生除零错误  
		iter = 0;
		x.disp_x_vector();
	}
	for (i = 0;i < max_iter;i++) //主循环
	{
		Vector mp =K_matrix*p;     //矩阵-向量乘
		double alpha = zr/v_d.dot(mp,p);  //当且仅当r=0，除数才为零
		x += alpha*p;         //更新迭代解
		r -= alpha*mp;        //更新残差
		if(r.max_norm() <= stp)	break; //如果已收敛，停机
		double zrold = zr;
		zr = v_d.dot(r,r);          //点积
		p = r + (zr/zrold)*p;   //当且仅当r=0，zrold=0
	}	
	eps = r.max_norm();
	x.disp_x_vector();
	return r.max_norm();
}


void FiniteEM::dispK_matrix()
{
	fstream output;
	output.open("f:\\K_matrix.txt",ios::out);
	for(int i=0;i<41;i++)
	{
		for(int j=0;j<41;j++)
		{
			output<<setw(6)<<pK_matrix[i][j]<<setw(6);
		}
		output<<endl<<endl;
	}
}

void FiniteEM::dispconK_matrix()            //压缩后的K矩阵
{
	fstream output;
	output.open("f:\\condenseK_matrix.txt",ios::out);
	for(int i=0;i<26;i++)
	{
		for(int j=0;j<26;j++)
		{
			output<<setw(6)<<pK_matrix[i][j]<<setw(6);
		}
		output<<endl<<endl;
	}
}

void FiniteEM::dispCon_matrix()              //联系数组
{
	fstream outCon_matrix;
	outCon_matrix.open("f:\\con_matrix.txt",ios::out);
	for(int i=0;i<3;i++)
	{
		for(int j=0;j<58;j++)
		{
			outCon_matrix<<setw(5)<<pArr2D[i][j]<<setw(5);
		}
		outCon_matrix<<endl<<endl<<endl;
	}
}

void FiniteEM::disparray_D3()
{
	fstream outCon_a_D3;
	outCon_a_D3.open("f:\\array_D3.txt",ios::out);
	for(int i=0;i<58;i++)
	{
		outCon_a_D3<<i+1<<":";
		for(int j=0;j<3;j++)	
		{		
			for(int k=0;k<3;k++)		
			{
			outCon_a_D3<<setw(20)<<pMatrix_D3[i][j][k];		
			}	
		}
		outCon_a_D3<<endl<<endl<<endl;
	}
	outCon_a_D3<<endl<<endl<<endl;
}

double triArea(int numTri)
{
	if(numTri<40)
	{
		return 0.0002;
	}
	else 
		return 0.0001;
}