pp.hpp

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

#ifndef PP_HPP
#define PP_HPP

#include<fstream.h>
#include"position.hpp"
#include<math.h>
#include<iomanip.h>

#define uu(x,y) ((x)*(x)*(x)+(y)*(y)*(y))            //验证的原方程
//
//
//
/* 存在问题:节点个数最好为一个可整开方的数，
 * 如4，9，16，25，36，49，64，81，100，。。
 *
 */

class Pp{
public:
	Pp(double (*func)(double ,double),char *filename=NULL);
	~Pp();
	void about();  //
	void help();   //
	
	//double Uh();
	void Set_SP(int p=0);
	void do_it();
	void Inputfile();     //从文件输入数据
	void OutFile();        //输出到文件
	double T6(double x1,double y1,double x2,double y2,double x3,double y3); //三角形的高斯6点积分计算式
	double a(Pos i,Pos j);             //与i，j节点有关的三角形求“和”
	double f_j(Pos j);                 //求(f,j)
	void formMat();                    //形成刚度矩阵
	void Matsolve();                   //用LU分解求解矩阵

	void setE1(double (*func)(double ,double)=NULL,int flag1=0);  // 指向边值值函数
	void setE2(double (*func)(double ,double)=NULL,int flag2=0);
	void setE3(double (*func)(double ,double)=NULL, int flag3=0);
	void setE4(double (*func)(double ,double)=NULL, int flag4=0);

private:	
	int next(int i,int j);             //判断两点关系，同一点:0,相异邻接:1,其余:-1
	double s(Pos i,Pos j,Pos k);      //面积
	double ai(Pos i,Pos j,Pos k);     
	double bi(Pos i,Pos j,Pos k);
	double ci(Pos i,Pos j,Pos k);
	double point_i(Pos i);             //计算点i的a(i,i)
	double edg_ij(Pos i,Pos j);        //计算点i_j的a(i,j)
	//double EP_ij(Pos i,Pos j,Pos k);   //与边节点有关的节点求和
	//double EP_ji(Pos i,Pos k,Pos j);   //与边节点有关的节点求和
	bool is(int i,int j);            //判断是否有效节点
    /**************************************/
	double (*f)(double ,double);  //指向函数指针
	char *fn;             //输入文件名
	char *ofn;            //输出文件名
	int nodN;             //节点数
	int ord;              //存放节点表的二维数组的阶
	Pos **A;               //节点表
	Pos *INN;
	double **Mat;          //刚度矩阵
	double *b;             //右边项
	double *X;             //最后的解
	int xN;                //解的个数(内点个数)

	int flag[4];           //标记边值条件类型
	int end[4];            //标记各边值的结束位置
	double *eV;            //存放边值条件
	int eN;                //边节点个数

	int show_precision;    //是否要输出精确解


	double (*Ef)(double ,double);   
	double (*Ef1)(double ,double);  //指向边值函数指针
	double (*Ef2)(double ,double);
	double (*Ef3)(double ,double);
	double (*Ef4)(double ,double);

};

Pp::Pp(double (*func)(double ,double),char *filename)
{
	show_precision=0;
	f=func;
	if(filename==NULL){
		cout<<"◎输入数据文件名:";
		fn=new char[50];
		cin>>fn;
	}
	else fn=filename;	
	for(int i=0;i<4;i++)flag[i]=0;
	return;

}

Pp::~Pp(){
	for(int i=0;i<ord;i++)
		delete[] A[i];
	    delete[] A;
		delete[] b;
		delete[] X;
		delete[] eV;
		delete[] INN;
}

//形成方程，求解任意点的值


void Pp::Set_SP(int p)
{
	if(p>=0)show_precision=p;
}
void Pp::about()
{
	cout<<"◎"<<endl;
	cout<<"◎欢迎使用、演示、...本二维问题有限元求解程序"<<endl;
	cout<<"◎"<<endl;
	cout<<"◎*******************版权所有*******************"<<endl;
	cout<<"◎第三组组员:(.....)"<<endl;
	cout<<"◎                                 2004-11"<<endl;
	int i=0;
	cout<<"◎"<<endl;
	while(i++<38)cout<<"◎";
	cout<<"◎"<<endl;
	cout<<"◎"<<endl;
}

void Pp::help()
{

	cout<<"◎本程序使用三角元方法求解二维问题"<<endl;
	cout<<"◎如果你要求解问题的区域为矩形,那么恭喜你,程序将回自动帮你划分单元"<<endl;
	cout<<"◎如果不是,那也不要紧,你将有机会自己(划分+输入)各节点坐标(文件输入)"<<endl;
	cout<<"◎第一行输入节点个数(n-1)^2,最好如2^2=4,3^2=9,4^2=16,5^2=25,...,10^2=100,...个"<<endl;
	cout<<"◎划分需要一定技巧,你需要将区域边界分为四段(按边值条件而定)"<<endl;
	cout<<"◎然后将边界等分n段,接下来就容易了:)"<<endl;
	cout<<"◎当然,你可以在main.cpp文件中设置f(x,y)和各边值条件函数"<<endl;
	cout<<"◎0:齐次,1:一次菲齐次,2:2次菲齐次"<<endl;
	cout<<"◎在使用本软件过程中如遇到问题,请邮件至:[e-mail],你会得到满意答复"<<endl;
	cout<<"◎"<<endl;
	cout<<"◎－－－－－－－－－祝你愉快－－－－－－－－－－"<<endl;
	cout<<"◎"<<endl;
	int i=0;
	while(i++<38)cout<<"◎";
	cout<<"◎"<<endl;
	
}


void Pp::do_it()
{
	Inputfile();                  //输入
	formMat();                   //形成刚度矩阵
	Matsolve();                  //解－>矩阵
	OutFile();
}

void Pp::Inputfile()
{
	ifstream input;
	input.open(fn);
	if(input==NULL){
		cout<<"Open file failed!"<<endl;
		return;
	}
	input>>nodN;         //输入节点总数

	float temp;
	temp=(float)sqrt(nodN);
	ord=(int)temp;
	if(temp-ord>0)ord+=1;   //适当放大
	//xN=ord-2;
	int i,j,count;  
/************************边值条件相关*************************************/
           //边值条件
	count=eN=4*(ord-1);         //边值条件节点总数
	xN=nodN-eN;                 //内点个数
	eV=new double[count];
	for(i=0;i<count;i++)eV[i]=0.0;
	end[0]=ord;end[1]=2*ord;end[2]=2*ord+ord-2;end[3]=4*(ord-1);  //设置各边值边的结束位置

	
/**********************************************************************/
	A=new Pos*[ord];          //申请空间
	for(i=0;i<ord;i++)
		A[i]=new Pos[ord];

	
	for(i=0,count=0;i<ord;i++)                //输入
		for(j=0;j<ord&&count<nodN;j++){
			input>>A[i][j];
			count++;
			A[i][j].setN(count);             //设置编号
		}

	input.close();
	cout<<"◎从文件*"<<fn<<"*输入数据成功!:)"<<endl;
	cout<<"◎正在求解中..."<<endl;

	
	return;
}

void Pp::OutFile()
{
	cout<<"◎\n◎输入*输出文件*名:";
	ofn=new char[50];
	cin>>ofn;
	ofstream outfile;
	outfile.open(ofn);
	int i,j;
	outfile<<"刚度矩阵:右端项"<<endl;
	for(i=0;i<xN;i++){
		for(j=0;j<xN;j++){
			outfile<<Mat[i][j]<<" ";
		}
		outfile<<":"<<b[i];
		outfile<<endl;
	}
	outfile<<endl;
	double mm,nn;

	outfile<<"方程组的解为:";
	if(show_precision>0)outfile<<"             精确解为：              误差为:"<<endl;
	outfile<<endl;
	for(i=0;i<xN;i++){            //输出解
			outfile<<'x'<<i+1<<'=';
			outfile.fill(' ');
			outfile.setf(ios::left);
			outfile.width(25);
			outfile.precision(10);
			outfile<<(mm=X[i]);
			outfile.fill(' ');
			if(show_precision>0){
				outfile.setf(ios::left);
			outfile.width(25);
			outfile.precision(10);
			outfile<<(nn=uu(INN[i].getX(),INN[i].getY()));
			outfile.fill(' ');
			outfile.setf(ios::left);
			outfile.width(25);
			outfile.precision(10);
			outfile<<fabs(mm-nn);
			}
			outfile<<endl;
	}
	outfile.unsetf(ios::left);
	outfile.close();
	cout<<"◎\n◎求解结果已经被保存到文件*"<<ofn<<"*中,敬请查看!"<<endl<<"◎"<<endl;
	return;
}

//6点计算公式
double Pp::T6(double x1,double y1,double x2,double y2,double x3,double y3)
{
	double ret=0.0;
	double aa,bb,cc;
	aa=x2*y3-x3*y2;bb=y2-y3;cc=x3-x2;
	double s=((x2*y3-x3*y2)-(x1*y3-x3*y1)+(x1*y2-x2*y1))/2;
	//s=fabs(s)/2;
	//ccout<<"s"<<s<<endl;
	double temp1,temp2;
	temp1=0.659027622374092*x1+0.231933368553031*x2+0.109039009072877*x3;
	temp2=0.659027622374092*y1+0.231933368553031*y2+0.109039009072877*y3;
    ret+=(*f)(temp1,temp2)*(aa+bb*temp1+cc*temp2);                           //第1项
	temp1=0.659027622374092*x1+0.231933368553031*x3+0.109039009072877*x2;
	temp2=0.659027622374092*y1+0.231933368553031*y3+0.109039009072877*y2;
	ret+=(*f)(temp1,temp2)*(aa+bb*temp1+cc*temp2);                        //第2项
	temp1=0.659027622374092*x2+0.231933368553031*x1+0.109039009072877*x3;
	temp2=0.659027622374092*y2+0.231933368553031*y1+0.109039009072877*y3;
	ret+=(*f)(temp1,temp2)*(aa+bb*temp1+cc*temp2);                          //第3项
	temp1=0.659027622374092*x2+0.231933368553031*x3+0.109039009072877*x1;
	temp2=0.659027622374092*y2+0.231933368553031*y3+0.109039009072877*y1;
	ret+=(*f)(temp1,temp2)*(aa+bb*temp1+cc*temp2);                          //第4项
	temp1=0.659027622374092*x3+0.231933368553031*x1+0.109039009072877*x2;
	temp2=0.659027622374092*y3+0.231933368553031*y1+0.109039009072877*y2;
	ret+=(*f)(temp1,temp2)*(aa+bb*temp1+cc*temp2);                          //第5项
	temp1=0.659027622374092*x3+0.231933368553031*x2+0.109039009072877*x1;
	temp2=0.659027622374092*y3+0.231933368553031*y2+0.109039009072877*y1;
	ret+=(*f)(temp1,temp2)*(aa+bb*temp1+cc*temp2);                          //第6项
	ret=ret/6*fabs(s);                   //乘权
	return ret;	
}

//求a(i,j)
double Pp::a(Pos i,Pos j)      
{
	double ret;
	if(next(i.getN(),j.getN())==-1)ret=0.0;
	else if(next(i.getN(),j.getN())==0)ret=point_i(i);
	else ret=edg_ij(i,j);

	return ret/4.0;
}

//求(f,j)