main.cpp

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

#include"position.hpp"
#include"pp.hpp"
#include<time.h>

//到这里设置相关函数
double f(double x,double y)            //右端函数f(x,y)
{
	return -6*(x+y);//(x)*(y);
}

double rfa1(double x,double y)         //边值条件1,上横边
{
	return (x)*(x)*(x)+1;//(x)*(y);
}
double rfa2(double x,double y)         //边值条件2，下横边
{
	return (x)*(x)*(x)-1;//(x)*(y);
}
double rfa3(double x,double y)         //边值条件3，左横边
{
	return (y)*(y)*(y)-1;//(x)*(y);
}
double rfa4(double x,double y)          //边值条件4，右横边
{
	return (y)*(y)*(y)+1;//(x)*(y);
}

//划分矩形区域，参数分别为（左上横坐标，左上纵坐标，横坐标区间长度，纵坐标区间长度，
//等分份数/默认为9，输出划分区域的各点坐标/默认为in.txt);
void ini_val(double x,double y,double a,double b,int n=9,char *fn="in.txt")
{
	ofstream ini;
	ini.open(fn);
	cout<<"◎"<<endl;
	cout<<"◎你现在解的是矩形域有限元,请输入等分次数"<<endl;
	cout<<"◎(2,3,4,5,...,9,...等分,分别有9,16,25,36,...,100,...个节点!)"<<endl<<"◎>";
	cin>>n;
	while(n<=0)cin>>n;
	int i,j;
	ini<<(n+1)*(n+1)<<"\n";
	double h1=a/n,h2=b/n;
	double x1,y1;
	for(i=0;i<=n;i++){
		y1=y-i*h2;
		for(j=0;j<=n;j++){
			x1=x+j*h1;
			ini<<"("<<setprecision(10)<<x1<<","<<setprecision(10)<<y1<<")";
		}
		ini<<endl;
	}

	ini.close();
	cout<<"◎数据被保存到文件*"<<fn<<"*中"<<endl;
	cout<<"◎\n";
}


//主函数
void main()
{
	
	char c;
	int hu=0;	
	long t1,t2;
	t1=clock();
	do{
	
	if(hu>0){ //界面有关
		cout<<endl<<endl;
		for(int i=0;i<39;i++)cout<<"◎";
		cout<<endl<<"◎";
		cout<<endl<<"◎               欢迎继续使用本程序!"<<endl;
		cout<<"◎"<<endl;
		for(i=0;i<39;i++)cout<<"◎";
		cout<<endl<<"◎\n";
	}


	//double (*fp)(double ,double);
	//fp=&f;
	Pp p(&f,"in.txt");
	if(hu++==0){           //第一次使用时给出
		p.about();              //程序信息
		p.help();               //帮助信息
	}
	ini_val(-1,1,2,2);          //矩形域的自动剖分函数
    //设置边值条件函数,0:齐次,1:一次非齐次,2:2次非齐次
	p.setE1(&rfa1,1);          //上边
	p.setE2(&rfa2,1);          //下边
	p.setE3(&rfa3,1);          //左边
	p.setE4(&rfa4,1);           //右边

	p.Set_SP(1);                //在这里设置是否要(有)输出精确解,如果设置了大于0，
	                            //你还要到pp.hpp文件设置原函数
	p.do_it();                  //做它
	//p.Uh();                   //逼近求值

	/*********界面有关***********/
	t2=clock();	
	cout<<"◎共用时:"<<(t2-t1)/1000.0<<"(秒)"<<endl<<"◎";
	cout<<endl;
	
	cout<<"◎";cout<<"\n◎";
	cout<<"'y'/'Y'继续(其它结束):";
	cin.get(c);
	cin.get(c);
	}while(c=='y'||c=='Y');
	
	for(int i=0;i<39;i++)cout<<"◎";
	cout<<endl<<"◎";
	cout<<endl<<"◎          谢谢使用本有限元求解程序，欢迎下次再来!"<<endl;
	cout<<"◎"<<endl;
	for(i=0;i<39;i++)cout<<"◎";
	cout<<endl;
	return;
}