📄 twodimension.h

📁 偏微分方程数值解有限元法面向对象变分问题剖分问题边值处理误差分析椭圆方程
💻 H
📖 第 1 页 / 共 2 页
字号:
12 下一页
#ifndef Twodimension_h
#define Twodimension_h
#include "FinityElement.h"
#include "mymath.h"
#include "Point.h"
#define filename2 "twoinitial.txt"
double H;
Point T[50][50];
class Twodimension:public FinityElement<Point>{
private:
 	Point leftlow,rightlow,lefthigh,righthigh;
	double steplen,wideth,height;
	int Nrow,Ncol;//剖分点个数
	Point t[50][50];//自变量、
	double sol[50][50],exactsol[50][50],noncoef[50][50];//解、精确解、非齐次系数
	double stiffness[16][16][16][16];//刚度矩阵
    double exactf(Point p);//精确解函数
	void calstiffness();//计算刚度矩阵
	void calnoncoef();//计算非齐次项系数
	void calequation();//计算有限元方程组
	double upedge1(Point p);//第一上边值条件
	double downedge1(Point p);//第一下边值条件
	friend double leftedge2(Point p);//第二左边值条件
	double rightedge2(Point p);//第二右边值条件
	void egdeprocess();//边值处理
protected:
public:
	Twodimension();
	~Twodimension(){}
	void Equation();//方程
	void Function();//变分问题
	void Dispart();//剖分
	void Radix();//基函数
	void Finity();//有限元方程
	void Calculate();//计算
    void Exact();//精确解;
	void Error();//误差

};
Twodimension::Twodimension(){
    ifstream fin(filename2,ios::in|ios::nocreate);
	if(!fin){
		cout<<filename2<<" cannot be openned!"<<endl;exit(1);
	}
	cout<<"Input leftlow:"<<endl;fin>>leftlow;	
	cout<<"Input rightlow:"<<endl;fin>>rightlow;
	cout<<"Input rightehigh:"<<endl;fin>>righthigh;
	cout<<"Input lefthigh:"<<endl;fin>>lefthigh;
	cout<<endl<<"输入x方向的剖分点数Nrow（如4，5，8,10）:";cin>>Nrow;
	fin.close();
}
void Twodimension::Dispart(){
	int i,j;
	wideth=rightlow.x-leftlow.x;	
	H=steplen=wideth/Nrow;
	height=righthigh.y-rightlow.y;
	Ncol=int(height/steplen);
    for(i=0;i<Nrow+1;i++){
		for(j=0;j<Ncol+1;j++){
		{T[i][j].x=t[i][j].x=leftlow.x+i*steplen;
		 T[i][j].y=t[i][j].y=leftlow.y+j*steplen;
		}
		}
	}
}
void Twodimension::Radix(){//刚度矩阵初始化为0
	int i,j,k,l;
    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++)
	stiffness[i][j][k][l]=0;
}
double commonx1(Point p,int r,int l){//下面函数代表一系列小公共操作
	return (1+(p.x-T[r][l].x)/H)/H;
}
double negcommonx1(Point p,int r,int l){
	return -commonx1(p,r,l);
}
double commonx2(Point p,int r,int l){
	return (1-(p.x-T[r][l].x)/H)/H;
}
double negcommonx2(Point p,int r,int l){
	return -commonx2(p,r,l);
}
double commony1(Point p,int r,int l){
	return (1+(p.y-T[r][l].y)/H)/H;
}
double negcommony1(Point p,int r,int l){
	return -commony1(p,r,l);
}
double commony2(Point p,int r,int l){
	return (1-(p.y-T[r][l].y)/H)/H;
}
double negcommony2(Point p,int r,int l){
	return -commony2(p,r,l);
}
double F11(Point p,int r,int l){//下面函数为计算刚度矩阵调用
	return
		negcommony2(p, r, l)*negcommony2(p, r, l)+negcommonx2(p, r, l)*negcommonx2(p, r, l);		
}
double F12(Point p,int r,int l){
	return 
		commony2(p, r, l)*commony2(p, r, l)+negcommonx1(p, r, l)*negcommonx1(p, r, l);
}
double F13(Point p,int r,int l){
	return 
		commony1(p, r, l)*commony1(p, r, l)+commonx1(p, r, l)*commonx1(p, r, l);
}
double F14(Point p,int r,int l){
	return 
		negcommony1(p, r, l)*negcommony1(p, r, l)+commonx2(p, r, l)*commonx2(p, r, l);
}
double F21(Point p,int r,int l){
	return negcommony2(p, r, l)*commony2(p, r, l)+negcommonx2(p, r, l)*negcommonx1(p, r+1, l);
}
double F22(Point p,int r,int l){
	return negcommony1(p, r, l)*commony1(p, r, l)+commonx2(p, r, l)*commonx1(p, r+1, l);
}
double F3(Point p,int r,int l){
	return negcommony2(p, r, l)*commony1(p, r, l+1)+negcommonx2(p, r, l)*commonx1(p, r+1, l);
}
double F41(Point p,int r,int l){
	return negcommony2(p, r, l)*negcommony1(p, r, l+1)+negcommonx2(p, r, l)*commonx2(p, r, l);
}
double F42(Point p,int r,int l){
	return commony2(p, r, l)*commony1(p, r, l+1)+negcommonx1(p, r, l)*commonx1(p, r, l);
}
double F5(Point p,int r,int l){
	return commony2(p, r, l)*negcommony1(p, r, l+1)+negcommonx1(p, r, l)*commonx2(p, r-1, l);
}
double F61(Point p,int r,int l){
	return commony2(p, r, l)*negcommony2(p, r, l)+negcommonx1(p, r, l)*negcommonx2(p, r-1, l);
}
double F62(Point p,int r,int l){
	return commony1(p, r, l)*negcommony1(p, r, l)+commonx1(p, r, l)*commonx2(p, r-1, l);
}
double F7(Point p,int r,int l){
	return commony1(p, r, l)*negcommony2(p, r, l-1)+commonx1(p, r, l)*negcommonx2(p, r-1, l);
}
double F81(Point p,int r,int l){
	return commony1(p, r, l)*commony2(p, r, l-1)+commonx1(p, r, l)*negcommonx1(p, r, l);
}
double F82(Point p,int r,int l){
	return negcommony1(p, r, l)*negcommony2(p, r, l-1)+commonx2(p, r, l)*negcommonx2(p, r, l);
}
double F9(Point p,int r,int l){
	return negcommony1(p, r, l)*commony2(p, r, l-1)+commonx2(p, r, l)*negcommonx1(p, r+1, l);
}
void Twodimension::calstiffness(){//计算刚度矩阵
	int i,j;
	//处理所有中央点
	for(i=1;i<Nrow;i++)
	{for(j=1;j<Ncol;j++)
	{
	stiffness[i][j][i][j]=Guass4(F11,i,j,T[i][j],T[i+1][j],T[i+1][j+1],T[i][j+1])+
	Guass4(F12,i,j,T[i-1][j],T[i][j],T[i][j+1],T[i-1][j+1])+
	Guass4(F13,i,j,T[i-1][j-1],T[i][j-1],T[i][j],T[i-1][j])+
	Guass4(F14,i,j,T[i][j-1],T[i+1][j-1],T[i+1][j],T[i][j]);
	stiffness[i][j][i+1][j]=Guass4(F21,i,j,T[i][j],T[i+1][j],T[i+1][j+1],T[i][j+1])+
	Guass4(F22,i,j,T[i][j-1],T[i+1][j-1],T[i+1][j],T[i][j]);
	stiffness[i][j][i+1][j+1]=Guass4(F3,i,j,T[i][j],T[i+1][j],T[i+1][j+1],T[i][j+1]);
	stiffness[i][j][i][j+1]=Guass4(F41,i,j,T[i][j],T[i+1][j],T[i+1][j+1],T[i][j+1])+
    Guass4(F42,i,j,T[i-1][j],T[i][j],T[i][j+1],T[i-1][j+1]);
	stiffness[i][j][i-1][j+1]=Guass4(F5,i,j,T[i-1][j],T[i][j],T[i][j+1],T[i-1][j+1]);
 	stiffness[i][j][i-1][j]=Guass4(F61,i,j,T[i-1][j],T[i][j],T[i][j+1],T[i-1][j+1])+
    Guass4(F62,i,j,T[i-1][j-1],T[i][j-1],T[i][j],T[i-1][j]);
	stiffness[i][j][i-1][j-1]=Guass4(F7,i,j,T[i-1][j-1],T[i][j-1],T[i][j],T[i-1][j]);
  	stiffness[i][j][i][j-1]=Guass4(F81,i,j,T[i-1][j-1],T[i][j-1],T[i][j],T[i-1][j])+
    Guass4(F82,i,j,T[i][j-1],T[i+1][j-1],T[i+1][j],T[i][j]);
	stiffness[i][j][i+1][j-1]=Guass4(F9,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];
    stiffness[0][0][0][0]=Guass4(F11,0,0,T[0][0],T[0+1][0],T[0+1][0+1],T[0][0+1]);
    stiffness[0][0][1][1]=Guass4(F3,0,0,T[0][0],T[0+1][0],T[0+1][0+1],T[0][0+1]);
    stiffness[0][0][0][1]=Guass4(F21,0,0,T[0][0],T[1][0],T[1][1],T[0][1]);
    stiffness[0][0][1][0]=Guass4(F41,0,0,T[0][0],T[1][0],T[1][1],T[0][1]);
	//左下点结束
    stiffness[Nrow][0][Nrow][0]=
		Guass4(F12,Nrow,0,T[Nrow-1][0],T[Nrow][0],T[Nrow][0+1],T[Nrow-1][0+1]);
	stiffness[Nrow][0][Nrow-1][1]=
		Guass4(F5,Nrow,0,T[Nrow-1][0],T[Nrow][0],T[Nrow][0+1],T[Nrow-1][0+1]);
    stiffness[Nrow][0][Nrow][1]=
		Guass4(F42,Nrow,0,T[Nrow-1][0],T[Nrow][0],T[Nrow][0+1],T[Nrow-1][0+1]);
    stiffness[Nrow][0][Nrow-1][0]=
		Guass4(F61,Nrow,0,T[Nrow-1][0],T[Nrow][0],T[Nrow][0+1],T[Nrow-1][0+1]);
    //右下点结束
    stiffness[Nrow][Ncol][Nrow][Ncol]=
		Guass4(F13,Nrow,Ncol,T[Nrow-1][Ncol-1],T[Nrow][Ncol-1],T[Nrow][Ncol],T[Nrow-1][Ncol]);
	stiffness[Nrow][Ncol][Nrow-1][Ncol-1]=
		Guass4(F7,Nrow,Ncol,T[Nrow-1][Ncol-1],T[Nrow][Ncol-1],T[Nrow][Ncol],T[Nrow-1][Ncol]);
     stiffness[Nrow][Ncol][Nrow-1][Ncol]=
		Guass4(F62,Nrow,Ncol,T[Nrow-1][Ncol-1],T[Nrow][Ncol-1],T[Nrow][Ncol],T[Nrow-1][Ncol]);
    stiffness[Nrow][Ncol][Nrow][Ncol-1]=
		Guass4(F81,Nrow,Ncol,T[Nrow-1][Ncol-1],T[Nrow][Ncol-1],T[Nrow][Ncol],T[Nrow-1][Ncol]);
    //右上点结束
    stiffness[0][Ncol][0][Ncol]=
	Guass4(F14,0,Ncol,T[0][Ncol-1],T[0+1][Ncol-1],T[0+1][Ncol],T[0][Ncol]);
	stiffness[0][Ncol][1][Ncol-1]=
		Guass4(F9,0,Ncol,T[0][Ncol-1],T[0+1][Ncol-1],T[0+1][Ncol],T[0][Ncol]);
     stiffness[0][Ncol][1][Ncol]=
		Guass4(F22,0,Ncol,T[0][Ncol-1],T[0+1][Ncol-1],T[0+1][Ncol],T[0][Ncol]);
    stiffness[0][Ncol][0][Ncol-1]=
		Guass4(F82,0,Ncol,T[0][Ncol-1],T[0+1][Ncol-1],T[0+1][Ncol],T[0][Ncol]);
    //左上点结束

	//处理下边界点
	for(i=1;i<Nrow;i++){
    stiffness[i][0][i][0]=Guass4(F11,i,0,T[i][0],T[i+1][0],T[i+1][0+1],T[i][0+1])+
	Guass4(F12,i,0,T[i-1][0],T[i][0],T[i][0+1],T[i-1][0+1]);
	stiffness[i][0][i+1][0]=Guass4(F21,i,0,T[i][0],T[i+1][0],T[i+1][0+1],T[i][0+1]);
	stiffness[i][0][i+1][0+1]=Guass4(F3,i,0,T[i][0],T[i+1][0],T[i+1][0+1],T[i][0+1]);
 	stiffness[i][0][i][0+1]=Guass4(F41,i,0,T[i][0],T[i+1][0],T[i+1][0+1],T[i][0+1])+
    Guass4(F42,i,0,T[i-1][0],T[i][0],T[i][0+1],T[i-1][0+1]);
	stiffness[i][0][i-1][0+1]=Guass4(F5,i,0,T[i-1][0],T[i][0],T[i][0+1],T[i-1][0+1]);
	stiffness[i][0][i-1][0]=Guass4(F61,i,0,T[i-1][0],T[i][0],T[i][0+1],T[i-1][0+1]);
	}
	//处理右边界点
	for(j=1;j<Ncol;j++){
    stiffness[Nrow][j][Nrow][j]=
	Guass4(F12,Nrow,j,T[Nrow-1][j],T[Nrow][j],T[Nrow][j+1],T[Nrow-1][j+1])+
	Guass4(F13,Nrow,j,T[Nrow-1][j-1],T[Nrow][j-1],T[Nrow][j],T[Nrow-1][j]);
	stiffness[Nrow][j][Nrow][j+1]=
    Guass4(F42,Nrow,j,T[Nrow-1][j],T[Nrow][j],T[Nrow][j+1],T[Nrow-1][j+1]);
	stiffness[Nrow][j][Nrow-1][j+1]=
		Guass4(F5,Nrow,j,T[Nrow-1][j],T[Nrow][j],T[Nrow][j+1],T[Nrow-1][j+1]);
	stiffness[Nrow][j][Nrow-1][j]=
		Guass4(F61,Nrow,j,T[Nrow-1][j],T[Nrow][j],T[Nrow][j+1],T[Nrow-1][j+1])+
    Guass4(F62,Nrow,j,T[Nrow-1][j-1],T[Nrow][j-1],T[Nrow][j],T[Nrow-1][j]);
	stiffness[Nrow][j][Nrow-1][j-1]=
		Guass4(F7,Nrow,j,T[Nrow-1][j-1],T[Nrow][j-1],T[Nrow][j],T[Nrow-1][j]);
 	stiffness[Nrow][j][Nrow][j-1]=
		Guass4(F81,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--){
	stiffness[i][Ncol][i][Ncol]=
	Guass4(F13,i,Ncol,T[i-1][Ncol-1],T[i][Ncol-1],T[i][Ncol],T[i-1][Ncol])+
	Guass4(F14,i,Ncol,T[i][Ncol-1],T[i+1][Ncol-1],T[i+1][Ncol],T[i][Ncol]);
	stiffness[i][Ncol][i+1][Ncol]=
	Guass4(F22,i,Ncol,T[i][Ncol-1],T[i+1][Ncol-1],T[i+1][Ncol],T[i][Ncol]);
	stiffness[i][Ncol][i-1][Ncol]=
    Guass4(F62,i,Ncol,T[i-1][Ncol-1],T[i][Ncol-1],T[i][Ncol],T[i-1][Ncol]);
	stiffness[i][Ncol][i-1][Ncol-1]=
		Guass4(F7,i,Ncol,T[i-1][Ncol-1],T[i][Ncol-1],T[i][Ncol],T[i-1][Ncol]);
 	stiffness[i][Ncol][i][Ncol-1]=
		Guass4(F81,i,Ncol,T[i-1][Ncol-1],T[i][Ncol-1],T[i][Ncol],T[i-1][Ncol])+
12 下一页
💿 文件大小 146 K
👤 上传用户 junbo2009
📂 所属分类数学计算
🏷️ 相关标签

#分 #偏微分方程 #数值 #对象
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -