pp.hpp

微分方程数值解法实验--二维有限元（用C++实现）

double Pp::f_j(Pos j)
{
	double ret=0.0;
	int k1,k2,l1,l2,m1,m2;
	k2=(j.getN()-1)%ord;
	k1=(j.getN()-1-k2)/ord;

	int u=0;
	while(u++<6){
		if(u==1){l1=k1;l2=k2-1;m1=k1+1;m2=k2-1;}
		else {l1=m1;l2=m2;}
		if(u==2){m1=k1+1;m2=k2;}
		if(u==3){m1=k1;m2=k2+1;}
		if(u==4){m1=k1-1;m2=k2+1;}
		if(u==5){m1=k1-1;m2=k2;}
		if(u==6){m1=k1;m2=k2-1;}
	    ret+=T6(A[k1][k2].getX(),A[k1][k2].getY(),A[l1][l2].getX(),A[l1][l2].getY(),A[m1][m2].getX(),A[m1][m2].getY())/s(A[k1][k2],A[l1][l2],A[m1][m2]);
	}

	return ret/2.0;

}

void Pp::formMat()
{
	Mat=new double*[xN];          //申请空间
	int i,j;
	for(i=0;i<xN;i++)
		Mat[i]=new double[xN];
	b=new double[xN];
	X=new double[xN];

	int count=0;
	Pos *temp;                 //内点数组
    Pos *e;                    //边界节点数组

            
	temp=new Pos[xN]; 
	e=new Pos[eN];
	INN=new Pos[xN];

	for(i=1,count=0;i<ord-1;i++){
		for(j=1;(j<ord-1)&&(i*ord+j<nodN-1);j++,count++)
		{
			INN[count]=A[i][j];
			temp[count]=A[i][j];
		}
	}

	//放边节点坐标
	for(j=0;j<ord;j++){e[j]=A[0][j];e[end[0]+j]=A[ord-1][j];}
	for(j=1;j<ord-1;j++){e[end[1]+j-1]=A[j][0];e[end[2]+j-1]=A[j][ord-1];}
    //边值
	for(i=0;i<4;i++){
		if(flag[i]!=0){
			if(i<2){
				for(j=0;j<ord;j++){
					if(i==0)eV[i*ord+j]=Ef1(e[i*ord+j].getX(),e[i*ord+j].getY());//input>>eV[i*ord+j]>>c;
					else if(i==1)eV[i*ord+j]=Ef2(e[i*ord+j].getX(),e[i*ord+j].getY());
				}
			}
			else{
				for(j=0;j<ord-2;j++){
					if(i==2)eV[2*ord+(i-2)*(ord-2)+j]=Ef3(e[2*ord+(i-2)*(ord-2)+j].getX(),e[2*ord+(i-2)*(ord-2)+j].getY());
					if(i==3)eV[2*ord+(i-2)*(ord-2)+j]=Ef4(e[2*ord+(i-2)*(ord-2)+j].getX(),e[2*ord+(i-2)*(ord-2)+j].getY());
				}
			}
		}
	}

	//for(i=0;i<eN;i++)cout<<"eV:"<<eV[i]<<endl;cout<<endl;
	for(i=0;i<count;i++){               //内节点
		for(j=0;j<count;j++){
			Mat[i][j]=a(temp[j],temp[i]);
		}
		double r=0.0;              //处理第一边值条件
		for(int l=0;l<4;l++){
			if(flag[l]==1){
				int k=l-1>=0?end[l-1]:0;
				for(;k<end[l];k++)
					r+=a(e[k],temp[i])*eV[k];
			}                      //-----------------
		}
		b[i]=f_j(temp[i])-r;
	}

	delete []temp;
	delete []e;
}

void Pp::Matsolve()
{
	int m,i,j,k;
	int n=xN;
	double **aa;

	aa=new double*[n];              //中间矩阵
	for(i=0;i<n;i++)
		aa[i]=new double[n];
	X=new double[n];

	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			aa[i][j]=Mat[i][j];
		X[i]=b[i];
	}

	for(m=0;m<n;m++)
		{
			for(i=m;i<n;i++)
				for(k=0;k<m;k++)
					aa[i][m]-=aa[i][k]*aa[k][m];
			for(j=m+1;j<n;j++)
			{
				for(k=0;k<m;k++)
					aa[m][j]-=aa[m][k]*aa[k][j];
				aa[m][j]/=aa[m][m];
			}
		}
		for(i=0;i<n;i++)
		{
			for(k=0;k<i;k++)
				X[i]-=aa[i][k]*X[k];
			X[i]/=aa[i][i];
		}
		for(i=n-1;i>=0;i--)
		{
			for(k=i+1;k<n;k++)
				X[i]-=aa[i][k]*X[k];
		}
		for(i=0;i<n;i++)           //回收
			delete[] aa[i];
		delete[] aa;
		return;
}
/*******************设置边值条件*********************/
void Pp::setE1(double (*func)(double,double),int flag1)
{
	flag[0]=flag1;
	Ef1=func;
	return;
}

void Pp::setE2(double (*func)(double ,double), int flag2)
{
	flag[1]=flag2;
	Ef2=func;
	return;
}
void Pp::setE3( double (*func)(double ,double),int flag3)
{
	flag[2]=flag3;
	Ef3=func;
	return;
}
void Pp::setE4( double (*func)(double ,double),int flag4)
{
	flag[3]=flag4;
	Ef4=func;
	return;
}
/***********************private functions********************************/
int Pp::next(int i,int j)
{
	int ret=-1;
	if((i<=0||i>nodN)||(j<=0||j>nodN))ret=-1;
	if(i==j)ret=0;
	else {
	int k1,k2,l1,l2;
	k2=(i-1)%ord;
	k1=(i-1-k2)/ord;
	l2=(j-1)%ord;
	l1=(j-1-l2)/ord;
	int temp;
	if(k2>l2){temp=k1;k1=l1;l1=temp;temp=k2;k2=l2;l2=temp;}
	if(k2==l2){if((k1==l1-1)||(k1==l1+1))ret=1;}
	else {if(((k1==l1)&&(k2==l2-1))||((k1==l1+1)&&(k2==l2-1)))ret=1;}
	}

	return ret;
}

double Pp::s(Pos i,Pos j,Pos k)
{
	double x1,y1,x2,y2,x3,y3;
	x1=i.getX();y1=i.getY();
	x2=j.getX();y2=j.getY();
	x3=k.getX();y3=k.getY();
	double temp=fabs((x2*y3-x3*y2)-(x1*y3-x3*y1)+(x1*y2-x2*y1))/2;
	return temp;
}

double Pp::ai(Pos i,Pos j,Pos k)
{
	double x2,y2,x3,y3;
	x2=j.getX();y2=j.getY();
	x3=k.getX();y3=k.getY();
	return x2*y3-x3*y2;
}

double Pp::bi(Pos i,Pos j,Pos k)
{
	double y2,y3;
	y2=j.getY();y3=k.getY();
	return y2-y3;
}

double Pp::ci(Pos i,Pos j,Pos k)
{
	double x2,x3;
	x2=j.getX();x3=k.getX();
	return x3-x2;
}

double Pp::point_i(Pos i)         //内点的共i点的三角形"和"
{
	double ret=0.0;
	int k1,k2,l1,l2,m1,m2;
	double b,c;
	k2=(i.getN()-1)%ord;
	k1=(i.getN()-1-k2)/ord;

	int u=0;
	while(u++<6){
		if(u==1){l1=k1;l2=k2-1;m1=k1+1;m2=k2-1;}
		else {l1=m1;l2=m2;}
		if(u==2){m1=k1+1;m2=k2;}
		if(u==3){m1=k1;m2=k2+1;}
		if(u==4){m1=k1-1;m2=k2+1;}
		if(u==5){m1=k1-1;m2=k2;}
		if(u==6){m1=k1;m2=k2-1;}
	    b=bi(i,A[l1][l2],A[m1][m2]);c=ci(i,A[l1][l2],A[m1][m2]);
	    ret+=(b*b+c*c)/fabs(s(i,A[l1][l2],A[m1][m2]));
	}

	return ret;
}

double Pp::edg_ij(Pos i,Pos j)
{
	double ret=0.0;
	int k1,k2,l1,l2,m11,m12,m21,m22;
	double b1,c1,b2,c2;
	k2=(i.getN()-1)%ord;
	k1=(i.getN()-1-k2)/ord;
	l2=(j.getN()-1)%ord;
	l1=(j.getN()-1-l2)/ord;

	if(k2==l2){
		if(k1==l1+1){
			m11=k1-1;m12=k2+1;m21=k1;m22=k2-1;
		}
		else{
			m11=k1+1;m12=k2-1;m21=k1;m22=k2+1;
		}
	}
	else if(k2>l2){
		if(k1==l1){
			m11=k1-1;m12=k2;m21=k1+1;m22=k2-1;
		}
		else {
			m11=k1;m12=k2-1;m21=k1+1;m22=k2;
		}
	}
	else {
		if(k1==l1){
		    m11=k1+1;m12=k2;m21=k1-1;m22=k2+1;
		}
		else {
			m11=k1;m12=k2+1;m21=k1-1;m22=k2;
		}
	}	
	b1=bi(i,A[m11][m12],j);c1=ci(i,A[m11][m12],j);b2=bi(j,i,A[m11][m12]);c2=ci(j,i,A[m11][m12]);
	ret+=(b1*b2+c1*c2)/fabs(s(i,j,A[m11][m12]));
    b1=bi(i,j,A[m21][m22]);c1=ci(i,j,A[m21][m22]);b2=bi(j,A[m21][m22],i);c2=ci(j,A[m21][m22],i);
	ret+=(b1*b2+c1*c2)/fabs(s(i,j,A[m11][m12]));

	return ret;
}


bool Pp::is(int i,int j)   //判断是否有效节点
{
	int c;
	c=i*ord+j;
	if(c>=1&&c<nodN)return true;
	return false;

}




#endif //PP_HPP