#ifndef Onedimension_h
#define Onedimension_h
#include "FinityElement.h"
#include "mymath.h"
#define filename1 "oneinitial.txt"
#define E 1.0e-8
double h;
double* tout;
class Onedimension:public FinityElement<double>{
private:
 	double height,low,high,steplen,left,right;//區間長度，端點，步長，左右邊值
	int N;//剖分點個數
	double *t,*sol,*exactsol,*noncoef;//自變量、解、精確解、非齊次系數
	double stiffness[100][100];//剛度矩陣
    double exactf(double x);//精確解函數
protected:
public:
	Onedimension();
	~Onedimension(){}
	void Equation();//方程
	void Function();//變分問題
	void Dispart();//剖分
	void Radix();//基函數
	void Finity();//有限元方程
	void Calculate();//計算
    void Exact();//精確解;
	void Error();//誤差
	friend double F(double x);
    //double f3(double x);
    friend double f4(double x);
    friend double f5(double x); 
};
Onedimension::Onedimension(){
    ifstream fin(filename1,ios::in|ios::nocreate);
	if(!fin){
		cout<<filename1<<" cannot be openned!"<<endl;exit(1);
	}
	cout<<"Input low:";fin>>low;	
	cout<<"Input high:";fin>>high;
	cout<<"Input leftedgevalue:";fin>>left;
	cout<<"Input rightedgevalue:";fin>>right;
	cout<<endl<<"輸入剖分點數 N（如5,10,20,40）:";cin>>N;
	fin.close();
}
void Onedimension::Dispart(){
	int i;
	height=high-low;
	h=steplen=height/N;
	t=new double[N+1];
	tout=new double[N+1];
	sol=new double[N+1];
	exactsol=new double[N+1];
	noncoef=new double[N+1];
	for(i=0;i<N+1;i++)
	t[i]=low+i*steplen;
	for(i=0;i<N+1;i++)
	tout[i]=t[i];
}
double F(double x){
	return 2*sin(x*pi/2);
	//return 0;
}
double f1(double x){
	return -1.0/h+h*pi*pi*(1-x)*x/4.0;
}
double f2(double x){
	return 1.0/h+h*pi*pi*x*x/4.0;
}
double f3(double x){
	return 1.0/h+h*pi*pi*(1-x)*(1-x)/4.0;
}
double f4(double x,int j){
	return F(tout[j-1]+h*x)*x;
}
double f5(double x,int j){
	return F(tout[j]+h*x)*(1-x);
}
void Onedimension::Radix(){//初始化剛度矩陣為0
	int i,j;
    for(i=0;i<N+1;i++){
		for(j=0;j<N+1;j++)
			stiffness[i][j]=0;
	}
	for(i=0;i<N+1;i++)
		sol[i]=0;
}
double Onedimension::exactf(double x){
	return (4.0*sin((pi*x)/2))/(pi*pi)+
		((left*exp(-pi/2.0)+2.0*right/pi)/(exp(pi/2)+exp(-pi/2)))*exp(pi*x/2.0)+
		((left*exp(pi/2.0)-2.0*right/pi)/(exp(pi/2)+exp(-pi/2)))*exp(-pi*x/2.0);
	//return 0;
}
void Onedimension::Exact(){
	int i;
	for(i=0;i<N+1;i++)
		exactsol[i]=exactf(t[i]);
}
void Onedimension::Calculate(){
	int j,i;
 	//以下是高斯勒讓德積分求剛度矩陣;
	//stiffness[1][1]=GL(f2,0,1,4)+GL(f3,0,1,4);因為有u[0],所以可以統一到下面的循環
	for(j=1;j<N;j++){
		stiffness[j-1][j]=GL(f1,0,1,4);
		stiffness[j][j]=GL(f2,0,1,4)+GL(f3,0,1,4);
		stiffness[j+1][j]=GL(f1,0,1,4);
	}//因為沒有u[n+1]所以要單獨考慮以下兩個元素
	stiffness[N-1][N]=GL(f1,0,1,4);
	stiffness[N][N]=GL(f2,0,1,4);
	//以下用帶兩個參數的龍貝格積分公式計算有限元方程的非齊次項系數
	for(j=1;j<N;j++)
	{noncoef[j]=h*(LBG2(f4,j,0,1)+LBG2(f5,j,0,1));}
	noncoef[N]=h*LBG2(f4,N,0,1);
	noncoef[1]=noncoef[1]-left*stiffness[0][1];//對左邊界作處理
    noncoef[N]=noncoef[N]+right*(1+(right-t[N])/h);//對右邊界作處理，這里有可能擴充p(x)
	delete []tout;
    cout<<endl<<"以下輸出剛度矩陣及其非齊次項系數："<<endl;
    for(i=1;i<N+1;i++){
		for(j=1;j<N+1;j++)
			cout<<stiffness[i][j]<<" ";cout<<setw(14)<<noncoef[i]<<endl;//
	}
	//用LU分解法解三對角系數方程組理論上精度不高,假如有a[0][0]=0則無法進行
    /*for(i=1;i<N+1;i++){
		for(j=1;j<N+1;j++)
			stiffness[i-1][j-1]=stiffness[i][j];
	}
	for(i=1;i<N+1;i++)
		noncoef[i-1]=noncoef[i]; 
    sol=LU(stiffness,noncoef,N);
    //GS(stiffness,noncoef,sol,N,15);//很難保證高斯勒讓德迭代法收斂
	for(i=N;i>0;i--)
	{sol[i]=sol[i-1];
	}
	sol[0]=left;*/
    double *a,*b,*c;
	a=new double[N+1];
	b=new double[N+1];
	c=new double[N+1];
	for(i=2;i<N+1;i++)a[i]=stiffness[i][i-1];
	for(i=1;i<N+1;i++)b[i]=stiffness[i][i];
	for(i=1;i<N;i++)c[i]=stiffness[i][i+1];
    sol[0]=left;
    Chase(sol, a, b, c,noncoef,N);//用追趕法解三對角系數方程組	
	delete []a;
	delete []b;
	delete []c;/**/
}
void Onedimension::Error(){
	double x=0;
	cout<<"以下輸出誤差分析："<<endl;
    cout<<setw(6)<<"自變量"<<setw(14)<<"有限元解"<<setw(14)<<"精確解"
		<<setw(16)<<"誤差"<<endl;
	for(int i=0;i<N+1;i++)
	{
		x=fabs(sol[i]-exactsol[i]);
		cout<<setw(6)<<t[i]<<setw(14)<<sol[i]<<setw(14)<<exactsol[i]<<
			setw(16)<<x<<endl;
	}

}
#endif