twodimension.h

偏微分方程数值解 有限元法 面向对象 变分问题 剖分问题 边值处理 误差分析 椭圆方程

    Guass4(F82,i,Ncol,T[i][Ncol-1],T[i+1][Ncol-1],T[i+1][Ncol],T[i][Ncol]);
	stiffness[i][Ncol][i+1][Ncol-1]=
		Guass4(F9,i,Ncol,T[i][Ncol-1],T[i+1][Ncol-1],T[i+1][Ncol],T[i][Ncol]);
 	}
	//处理左边界点
	for(j=Ncol-1;j>0;j--){
	stiffness[0][j][0][j]=Guass4(F11,0,j,T[0][j],T[0+1][j],T[0+1][j+1],T[0][j+1])+
	Guass4(F14,0,j,T[0][j-1],T[0+1][j-1],T[0+1][j],T[0][j]);
	stiffness[0][j][0+1][j]=Guass4(F21,0,j,T[0][j],T[0+1][j],T[0+1][j+1],T[0][j+1])+
	Guass4(F22,0,j,T[0][j-1],T[0+1][j-1],T[0+1][j],T[0][j]);
	stiffness[0][j][0+1][j+1]=Guass4(F3,0,j,T[0][j],T[0+1][j],T[0+1][j+1],T[0][j+1]);
   	stiffness[0][j][0][j+1]=Guass4(F41,0,j,T[0][j],T[0+1][j],T[0+1][j+1],T[0][j+1]);
	stiffness[0][j][0][j-1]=
    Guass4(F82,0,j,T[0][j-1],T[0+1][j-1],T[0+1][j],T[0][j]);
	stiffness[0][j][0+1][j-1]=Guass4(F9,0,j,T[0][j-1],T[0+1][j-1],T[0+1][j],T[0][j]);
	}
}
double nonf(Point p,int r,int l){//非齐次项函数
	//return 2;
	//return (pi*pi*sin(pi*p.x/2))/4+2;
	return (pi*cos(pi*(1-p.x)/4))/4+2;
}
double nf1(Point p,int r,int l){//下面函数代表一系列小公共操作
	return (1-((p.x-T[r][l].x)/H))*(1-((p.y-T[r][l].y)/H));
}
double nf2(Point p,int r,int l){
	return (1+((p.x-T[r][l].x)/H))*(1-((p.y-T[r][l].y)/H));
}
double nf3(Point p,int r,int l){
	return (1+((p.x-T[r][l].x)/H))*(1+((p.y-T[r][l].y)/H));
}
double nf4(Point p,int r,int l){
	return (1-((p.x-T[r][l].x)/H))*(1+((p.y-T[r][l].y)/H));
}
double NF1(Point p,int r,int l){//下面函数为计算非齐次项系数所调用
	return nf1( p, r, l)*nonf( p, r, l);
}
double NF2(Point p,int r,int l){
	return nf2( p, r, l)*nonf( p, r, l);
}
double NF3(Point p,int r,int l){
	return nf3( p, r, l)*nonf( p, r, l);
}
double NF4(Point p,int r,int l){
	return nf4( p, r, l)*nonf( p, r, l);
}
void Twodimension::calnoncoef(){//计算非齐次项系数
	int i,j;
	//处理所有中央点
	for(i=1;i<Nrow;i++)
	{for(j=1;j<Ncol;j++)
	{
	noncoef[i][j]=Guass4(NF1,i,j,T[i][j],T[i+1][j],T[i+1][j+1],T[i][j+1])+
	Guass4(NF2,i,j,T[i-1][j],T[i][j],T[i][j+1],T[i-1][j+1])+
	Guass4(NF3,i,j,T[i-1][j-1],T[i][j-1],T[i][j],T[i-1][j])+
	Guass4(NF4,i,j,T[i][j-1],T[i+1][j-1],T[i+1][j],T[i][j]);
	}
	}
	//处理第一第二边值交界点，有四个[0][0],[Nrow][0],[Nrow][Ncol],[0][Ncol];
    noncoef[0][0]=
		Guass4(NF1,0,0,T[0][0],T[0+1][0],T[0+1][0+1],T[0][0+1]);
	//左下点结束
    noncoef[Nrow][0]=
		Guass4(NF2,Nrow,0,T[Nrow-1][0],T[Nrow][0],T[Nrow][0+1],T[Nrow-1][0+1]);
    //右下点结束
    noncoef[Nrow][Ncol]=
		Guass4(NF3,Nrow,Ncol,T[Nrow-1][Ncol-1],T[Nrow][Ncol-1],T[Nrow][Ncol],T[Nrow-1][Ncol]);
     //右上点结束
	noncoef[0][Ncol]=
		Guass4(NF4,0,Ncol,T[0][Ncol-1],T[0+1][Ncol-1],T[0+1][Ncol],T[0][Ncol]);
    //左上点结束/**/
	//处理下边界点
	for(i=1;i<Nrow;i++){
	noncoef[i][0]=Guass4(NF1,i,0,T[i][0],T[i+1][0],T[i+1][0+1],T[i][0+1])+
	Guass4(NF2,i,0,T[i-1][0],T[i][0],T[i][0+1],T[i-1][0+1]);
	}
	//处理右边界点
	for(j=1;j<Ncol;j++){
	noncoef[Nrow][j]=
	Guass4(NF2,Nrow,j,T[Nrow-1][j],T[Nrow][j],T[Nrow][j+1],T[Nrow-1][j+1])+
	Guass4(NF3,Nrow,j,T[Nrow-1][j-1],T[Nrow][j-1],T[Nrow][j],T[Nrow-1][j]);
	}
	//处理上边界点
	for(i=Nrow-1;i>0;i--){
	noncoef[i][Ncol]=
	Guass4(NF3,i,Ncol,T[i-1][Ncol-1],T[i][Ncol-1],T[i][Ncol],T[i-1][Ncol])+
	Guass4(NF4,i,Ncol,T[i][Ncol-1],T[i+1][Ncol-1],T[i+1][Ncol],T[i][Ncol]);
	}
	//处理左边界点
	for(j=Ncol-1;j>0;j--){
	noncoef[0][j]=Guass4(NF1,0,j,T[0][j],T[0+1][j],T[0+1][j+1],T[0][j+1])+
	Guass4(NF4,0,j,T[0][j-1],T[0+1][j-1],T[0+1][j],T[0][j]);
	}
}
double Twodimension::upedge1(Point p){//第一上边值条件
	//return sin((p.x*pi)/2);
	return (4*(cos(pi*(1-p.x)/4)))/pi;
}
double Twodimension::downedge1(Point p){//第一下边值条件
	//return sin((p.x*pi)/2);
	return (4*(cos(pi*(1-p.x)/4)))/pi;
}
double leftedge2(Point p){//第二左边值条件du/dn的值
	return -1;
}
double Twodimension::rightedge2(Point p){//第二右边值条件
	return 0;
}
double leftdir1(double y,int j){//第一象限基函数
	return 1-(y-T[0][j].y)/H;
}
double Leftintetral1(double y,int j){//左边值法向方向导数与基函数积分
	return leftedge2(T[0][j])*leftdir1(y,j);
}
double leftdir4(double y,int j){//第四象限基函数
	return 1+(y-T[0][j].y)/H;
}
double Leftintetral4(double y,int j){//左边值法向方向导数与基函数积分
	return leftedge2(T[0][j])*leftdir4(y,j);
}
void Twodimension::egdeprocess(){//边值处理
	int i,j;
	for(i=0;i<Nrow+1;i++)
	{sol[i][0]=downedge1(t[i][0]);
	sol[i][Ncol]=upedge1(t[i][Ncol]);
	}
	for(i=1;i<Nrow;i++){//上下边界点
		noncoef[i][1]=noncoef[i][1]-
			sol[i][0]*stiffness[i][0][i][1]
			-sol[i-1][0]*stiffness[i-1][0][i][1]
			-sol[i+1][0]*stiffness[i+1][0][i][1];
		noncoef[i][Ncol-1]=noncoef[i][Ncol-1]-
			sol[i][Ncol]*stiffness[i][Ncol][i][Ncol-1]
			-sol[i-1][Ncol]*stiffness[i-1][Ncol][i][Ncol-1]
			-sol[i+1][Ncol]*stiffness[i+1][Ncol][i][Ncol-1];
	}
	noncoef[0][1]=noncoef[0][1]-//方程系数中左下边值点
		sol[0][0]*stiffness[0][0][0][1]
			-sol[0+1][0]*stiffness[0+1][0][0][1];
	noncoef[0][Ncol-1]=noncoef[0][Ncol-1]-//方程系数中左上边值点
		sol[0][Ncol]*stiffness[0][Ncol][0][Ncol-1]
			-sol[0+1][Ncol]*stiffness[0+1][Ncol][0][Ncol-1];
	noncoef[Nrow][1]=noncoef[Nrow][1]-//方程系数中左下边值点
		sol[Nrow][0]*stiffness[Nrow][0][Nrow][1]
			-sol[Nrow-1][0]*stiffness[Nrow-1][0][Nrow][1];
	noncoef[Nrow][Ncol-1]=noncoef[Nrow][Ncol-1]-//方程系数中左上边值点
		sol[Nrow][Ncol]*stiffness[Nrow][Ncol][Nrow][Ncol-1]
			-sol[Nrow-1][Ncol]*stiffness[Nrow-1][Ncol][Nrow][Ncol-1];
    for (j=1;j<Ncol;j++){
		noncoef[0][j]=noncoef[0][j]+LBG2(Leftintetral1,j,T[0][j].y,T[0][j+1].y)
			+LBG2(Leftintetral4,j,T[0][j-1].y,T[0][j].y);
	}
}
void Twodimension::calequation(){//计算有限元方程组
	int i,j,k,l;
	int N=(Nrow+1)*(Ncol+1);
    int NM=(Nrow+1)*(Ncol+1)-2*(Nrow+1);
	double *b=new double[N];
	double *x=new double[N];
	double a[125][125];
	for(i=0;i<125;i++)
	for(j=0;j<125;j++)
	a[i][j]=0;
	for(i=0;i<Nrow+1;i++)
		for(j=0;j<Ncol+1;j++)
		{b[i+j*(Nrow+1)]=noncoef[i][j];
		}
	for(i=0;i<Nrow+1;i++)
	for(j=0;j<Ncol+1;j++)
	for(k=0;k<Nrow+1;k++)
	for(l=0;l<Ncol+1;l++)
	{a[i+j*(Nrow+1)][k+l*(Nrow+1)]=stiffness[i][j][k][l];
    }
	double *B=new double[NM];
	double *X=new double[NM];
	double A[100][100];
	for(i=0;i<100;i++)
	for(j=0;j<100;j++)
	A[i][j]=0;
	for(i=0;i<NM;i++)
	B[i]=b[i+Nrow+1];//这里可能还有非齐次边界条件的扩充
	for(i=0;i<NM;i++)
	for(j=0;j<NM;j++)
	A[i][j]=a[j+Nrow+1][i+Nrow+1];
	for(i=0;i<NM;i++)
	X[i]=0;
    /*cout<<endl<<"以下输出刚度矩阵及其非齐次项系数："<<endl;
    for(i=0;i<NM;i++){
		for(j=0;j<NM;j++)
			cout<<A[i][j]<<" ";cout<<B[i]<<endl;
	}*/
    X=LU(A,B,NM);//
    //GS(A,B,X,NM,15);
    //for(i=0;i<Nrow+1;i++)
		//x[i]=0;
	for(i=Nrow+1,j=0;j<NM;i++,j++)
    x[i]=X[j];
	//for(;i<N;i++)
	//x[i]=0; 
	for(i=0;i<Nrow+1;i++)
		for(j=1;j<Ncol;j++)//0,+1
		{sol[i][j]=x[i+j*(Nrow+1)];
		}
	delete []b;
	delete []x;
	delete []B;
	delete []X;
}
void Twodimension::Calculate(){
    calstiffness();
	calnoncoef();
    egdeprocess();
	calequation();
}
double Twodimension::exactf(Point p){
	//return p.x-(p.y+1)*(p.y-1);
	//return sin((pi*p.x)/2)-(p.y+1)*(p.y-1);
	return  (4*(cos(pi*(1-p.x)/4)))/pi-(p.y+1)*(p.y-1);
}
void Twodimension::Exact(){
	int i,j;
	for(i=0;i<Nrow+1;i++){
		for(j=0;j<Ncol+1;j++)
		{exactsol[i][j]=exactf(t[i][j]);}
	} 
}
void Twodimension::Error(){
	int i,j;
	double x=0;
	cout<<"以下输出误差分析："<<endl;
    cout<<setw(11)<<"自变量(x,y)";
	cout<<setw(14)<<"有限元解"<<setw(14)<<"精确解"
		<<setw(14)<<"误差"<<endl;
	for(j=0;j<Ncol+1;j++)
	for(i=0;i<Nrow+1;i++)
	{
		x=fabs(sol[i][j]-exactsol[i][j]);
		cout<<t[i][j];
 		cout<<setiosflags(ios::fixed);
        cout.precision(8);
		cout<<setw(14)<<sol[i][j]<<setw(14)<<exactsol[i][j];
        cout<<resetiosflags(ios::fixed);
		cout<<setw(17)<<x<<endl;
	}

}
#endif