#include <iostream.h>
#include <string.h>
#include <ostream.h>
#include <fstream.h>
#define N 6                                  //x軸方向
#define M 6                                  //y軸方向
main()
{
	float GAM(float,float,float);              //導熱系數隨坐標的變化值 
	float Sp(float,float,float);               //不同坐標下源項的值
	float Sc(float,float,float);               //不同坐標下源項的值
	float density(float,float,float);
	float c(float,float,float);
	float max(float,float);
	int i,j,k,num;
	float X_CV[N],X_point[N+1],Y_CV[M],Y_point[M+1],aw[N][M],as[N][M],ae[N][M],an[N][M],ap[N][M],b[N][M],ap0[N][M];
	float dx,dy,lx=0.12,ly=0.12,u,e,esp,t,dt=5,maxt=60;
	float midx[N+1],midy[M+1],T[N+1][M+1],T0[N+1][M+1];
	dx=lx/N;                              //取相同的微元
	dy=ly/M;
	X_CV[0]=0;                            //控制容積法劃分網格
	for(i=1;i<N;i++)
	{
		X_CV[i]=dx*i;
		X_point[i]=0.5*(X_CV[i-1]+X_CV[i]);         //結點劃分
	}
	X_point[N]=N*dx;
        Y_CV[0]=0;
	for(i=1;i<M;i++)
	{
		Y_CV[i]=dx*i;
		Y_point[i]=0.5*(Y_CV[i-1]+Y_CV[i]);
	}
	Y_point[M]=M*dx;
	midx[0]=midx[N-1]=dx/2;
	for(i=1;i<(N-1);i++)
		midx[i]=dx;
	midy[0]=midy[M-1]=dy/2;
	for(i=1;i<(M-1);i++)
		midy[i]=dy;
    for(j=1;j<M;j++)                     //賦迭代初值
	{
		for(i=1;i<N;i++)
		{
			T[i][j]=280;
			T0[i][j]=280;
		}
	}		
	for(i=1;i<N;i++)                      //邊界條件
	{
		T[i][0]=250;
		T[i][M]=250;
	}
	for(i=0;i<=M;i++)
	{
		T[0][i]=370;
		T[N][i]=250;
	}
	ofstream ft("text2.txt");                                                             //結果輸入進文件
		ft<<"---------------------------計算結果為---------------------------"<<endl;
	for(num=0;num<=maxt/dt;num++)
	{
	t=dt*num;
    
	for(j=1;j<(M-1);j++)                   //求四個系數
	{
		for(i=0;i<N;i++)
		{
			aw[i][j]=dy*GAM(X_point[i],Y_point[j],t)/midx[i];
			ae[i+1][j]=dy*GAM(X_point[i+1],Y_point[j],t)/midx[i+1];
		}
	}
	for(i=1;i<(N-1);i++)
	{
		for(j=0;j<M;j++)
		{
			as[i][j]=dx*GAM(X_point[i],Y_point[j],t)/midy[j];
			an[i][j+1]=dx*GAM(X_point[i],Y_point[j+1],t)/midy[j+1];
		}
	}
	for(j=1;j<M;j++)                      //確定ap和b
	{
		for(i=1;i<N;i++)
		{
			ap0[i][j]=density(X_point[i],Y_point[j],t)*c(X_point[i],Y_point[j],t)*dx*dy/dt;
			ap[i][j]=aw[i-1][j]+ae[i+1][j]+as[i][j-1]+an[i][j+1]-dx*dy*Sp(X_point[i],Y_point[j],t)+ap0[i][j];
			b[i][j]=Sc(X_point[i],Y_point[j],t)*dx*dy+ap0[i][j]*T0[i][j];
		}
	}
	k=0;
	esp=10;
	while(k<100&&esp>0.08)               //高斯迭代
	{
		esp=0.08;
		for(j=1;j<M;j++)
		{
			for(i=1;i<N;i++)
			{
				u=T[i][j];
				T[i][j]=(aw[i-1][j]*T[i-1][j]+ae[i+1][j]*T[i+1][j]+an[i][j+1]*T[i][j+1]+as[i][j-1]*T[i][j-1]+b[i][j])/ap[i][j];
				T0[i][j]=T[i][j];
				e=T[i][j]-u;
				esp=max(e,esp);
			}
		}
			k=k+1;
	}
	if(k==100)                           //屏幕輸出
		cout<<"迭代法不收斂!"<<endl;
	else
	{
		ft<<"-------------------------------------------------------------------"<<endl;
		for(j=0;j<=M;j++)
		{
			for(i=0;i<=N;i++)
			{
				cout.width(10);
				cout<<T[i][j];
				ft.width(10);
				ft<<T[i][j];
			}
			cout<<endl;
			ft<<endl;
		}
	}
	}
	ft.close();
}

float GAM(float x,float y,float t)
{
	float u;
	u=171;
	return u;
}
float Sp(float x,float y,float t)
{
	float u;
	u=0*x+0*y;
	return u;
}
float Sc(float x,float y,float t)
{
	float u;
	u=0*x+0*y;
	return u;
}
float max(float x,float y)
{
	if(x>y)
		return x;
	else 
		return y;
}
float density(float x,float y,float t)
{
	float u;
	u=1746;
	return u;
}
float c(float x,float y,float t)
{
	float u;
	u=971;
	return u;
}