#include"iostream"
#include"fstream"
#include"math.h"
#define fi  100//點數,是到倍數加
#define fj  40//點數
#define di  1//點距一般為，修改做圖方程
#define dj  1//點距一般為
#define ve  10//速度
#define pw  10//壓力
#define rg  1.4//絕熱指數
#define Ma0    0.5 //進口馬赫數
#define v0   0//邊界速度
#define po   101325.1//出口壓力
#define fvu  30  //動力粘度
#define Rmg  287.1//氣體常數
#define pi   3.1415826//圓周率
#define data 0.5//切向動量適應常數
using std::cout;
using std::cin;
using std::endl;
using std::ofstream ;

double bar(double x1,double x2,double d);//求偏導數
double abr(double x1,double x2,double y1,double y2);//求abr
double jj(double xe,double xn,double ye,double yn);//轉換函數J
void kcc(int i,int j,double visx[fi+2][fj+2],double visy[fi+2][fj+2],double u,double v,double p,double rou,char k/*小寫字母*/,double (*kcc)[4]);//A,B矩陣
double temp(double p,double rou);//求溫度
void FmB(double (*kcca)[4],double *kccg,double *kccm);//矩陣乘以向量
void FmS(double *k,double *kbar,double cita);//向量乘以數值
void Fma(double *k,double *m,double *kdd,char kk);//向量相加
void ffu(double rou,double vu,double vv,double p,double *ubar);//求向量U
void fff(double rou,double vu,double vv,double p,double *ubar);//求向量F
void ffg(double rou,double vu,double vv,double p,double *ubar);//求向量G
void midd(double (*kcc1)[4],double (*kcc2)[4],double (*kcc)[4]);//求矩陣中值
void jubar(int i,int j,double visx[fi+2][fj+2],double visy[fi+2][fj+2],double rou[fi+2][fj+2],double vu[fi+2][fj+2],double vv[fi+2][fj+2],double p[fi+2][fj+2],double *ubar,double *fbar,double *gbar);//求轉換向量
void fadao(int i,int j,double visx[fi+2][fj+2],double visy[fi+2][fj+2],double *n);//求向導數
void uvp(double vu,double vv,double rou,double p,double *n,double *ubar);//固壁面邊界處理
void jin(double vu0,double vv0,double rou0,double p0,double vu,double vv,double rou,double p,double *n0,double *n,double *ubar);//進口
void chukou(double vu0,double vv0,double rou0,double p0,double vu,double vv,double rou,double p,double *n0,double *n,double *ubar);//出口

void main()
{
double x[fi][fj],y[fi][fj],xx1[fi+1][fj+1],yy1[fi+1][fj+1];
int i,j,km,kn;

ofstream out;
x[0][0]=0;
y[0][0]=0;



for(i=1;i<fi;i++)/*給x賦值*/
x[i][0]=x[i-1][0]+di;

for(i=0;i<fi;i++)
for(j=1;j<fj;j++)
x[i][j]=x[i][0];



for(j=1;j<fj;j++)/*給y賦值*/
y[0][j]=y[0][j-1]+dj;


for(i=1;i<fi;i++)
for(j=0;j<fj;j++)
y[i][j]=y[0][j];

out.open("wanggexy.txt");//將x,y坐標點寫進文件
out<<"title=\"sample mesh\"\nvariables=\"x\",\"y\"\nzone i="<<fi<<" j="<<fj<<" f=point\n"<<endl;
for(j=0;j<fj;j++)
for(i=0;i<fi;i++)

out<<x[i][j]<<" \t"<<y[i][j]<<endl;
out.close();



//求內點
//四頂點
xx1[0][0]=x[0][0];
yy1[0][0]=y[0][0];
xx1[fi][fj]=x[fi-1][fj-1];
yy1[fi][fj]=y[fi-1][fj-1];

xx1[fi][0]=x[fi-1][0];
yy1[fi][0]=y[fi-1][0];
xx1[0][fj]=x[0][fj-1];
yy1[0][fj]=y[0][fj-1];

for(i=1;i<fi;i++)//邊界點
for(j=1;j<fj;j++)
{

xx1[i][0]=(x[i][0]+x[i-1][0])/2;
xx1[i][fj]=(x[i][fj-1]+x[i-1][fj-1])/2;
xx1[fi][j]=(x[fi-1][j-1]+x[fi-1][j])/2;
xx1[0][j]=(x[0][j-1]+x[0][j])/2;



yy1[i][0]=(y[i][0]+y[i-1][0])/2;
yy1[i][fj]=(y[i][fj-1]+y[i-1][fj-1])/2;
yy1[fi][j]=(y[fi-1][j-1]+y[fi-1][j])/2;
yy1[0][j]=(y[0][j-1]+y[0][j])/2;

}

for(i=1;i<fi;i++)//中間點
for(j=1;j<fj;j++)
{

xx1[i][j]=(x[i][j]+x[i][j-1]+x[i-1][j]+x[i-1][j-1])/4;
yy1[i][j]=(y[i][j]+y[i][j-1]+y[i-1][j]+y[i-1][j-1])/4;

}

out.open("wanggexy2.txt");//將x,y坐標點寫進文件
out<<"title=\"sample mesh\"\nvariables=\"x\",\"y\"\nzone i="<<fi+1<<" j="<<fj+1<<" f=point\n"<<endl;
for(j=0;j<=fj;j++)
for(i=0;i<=fi;i++)

out<<xx1[i][j]<<" \t"<<yy1[i][j]<<endl;
out.close();
//實際網格
double le,se,du,dv,lt,st,lh,lh0;//圖形邊長
cout<<"long   le=";
cin>>le;
cout<<endl<<"short se=";
cin>>se;
cout<<endl;
double u[fi][fj],v[fi][fj],uu1[fi+1][fj+1],vv1[fi+1][fj+1];
u[0][0]=0;
v[0][0]=0;
u[fi-1][0]=u[0][0]+le;
v[fi-1][0]=v[0][0];
u[0][fj-1]=u[0][0];
v[0][fj-1]=v[0][0]+se;
du=le/(fi-1.);
dv=se/(fj-1.);
for(j=1;j<fj;j++)/*給寬邊賦值*/
{   
	u[0][j]=0;
	v[0][j]=v[0][j-1]+dv;
	u[fi-1][j]=u[0][0]+le;
	v[fi-1][j]=v[0][j];
}

for(i=1;i<fi;i++)//給長邊賦值
{
u[i][fj-1]=u[i-1][fj-1]+du;
v[i][fj-1]=v[0][0]+se;
}

lt=le/3.;
lh0=asin(0.5/1.3);
lh=lh0*2/((fi-1)/3.);

for(i=1;i<=(fi-1)/3;i++)
{
	v[i][0]=v[0][0];
	u[i][0]=u[i-1][0]+du;
	v[(fi-1)/3+i][0]=(1.3*cos(lh0-lh*i)-1.2)*le/3.;
	u[(fi-1)/3+i][0]=lt+(0.5+1.3*sin(-lh0+lh*i))*le/3.;
	v[2*(fi-1)/3+i][0]=v[0][0];
	u[2*(fi-1)/3+i][0]=u[i-1][0]+lt*2+du;   
}

for(i=1;i<fi-1;i++)//內部賦初值
for(j=1;j<fj-1;j++)

{
	u[i][j]=u[i][0];
	v[i][j]=v[0][j];
}

out.open("wanggexy3.txt");//將x,y坐標點寫進文件
out<<"title=\"sample mesh\"\nvariables=\"u\",\"v\"\nzone i="<<fi<<" j="<<fj<<" f=point\n"<<endl;
for(j=0;j<fj;j++)
for(i=0;i<fi;i++)

out<<u[i][j]<<" \t"<<v[i][j]<<endl;
out.close();
 
double k,kk,n,cc,aa,bb,b1,rr,xx,xy,yx,yy;
k=1;
n=100000;

while(k>se*10e-6&&n>1)
{
	cc=0;
for(i=1;i<fi-1;i++)//內部賦求解
for(j=1;j<fj-1;j++)
{

xx=bar(u[i+1][j],u[i-1][j],2.);
xy=bar(u[i][j+1],u[i][j-1],2.);
yx=bar(v[i+1][j],v[i-1][j],2.);
yy=bar(v[i][j+1],v[i][j-1],2.);
aa=abr(xy,xy,yy,yy);
bb=-abr(xx,xy,yx,yy);
rr=abr(xx,xx,yx,yx);


kk=2.*(aa+rr);
b1=aa*(u[i+1][j]+u[i-1][j])+bb*(u[i+1][j+1]-u[i+1][j-1]-u[i-1][j+1]+u[i-1][j-1])/2.+rr*(u[i][j+1]+u[i][j-1]);//一般的
if(kk==0) kk=1; 
b1=b1/kk;
cc=cc+fabs(b1-u[i][j]);
u[i][j]=b1;


b1=aa*(v[i+1][j]+v[i-1][j])+bb*(v[i+1][j+1]-v[i+1][j-1]-v[i-1][j+1]+v[i-1][j-1])/2.+rr*(v[i][j+1]+v[i][j-1]);
if(kk==0) kk=1; 
b1=b1/kk;
cc=cc+fabs(b1-v[i][j]);
v[i][j]=b1;

}

if(cc<=k)k=cc;
cout<<k<<"\t"<<cc<<endl;
n=n-1;
}
out.open("wanggexy4.txt");//將x,y坐標點寫進文件
out<<"title=\"sample mesh\"\nvariables=\"u\",\"v\"\nzone i="<<fi<<" j="<<fj<<" f=point\n"<<endl;
for(j=0;j<fj;j++)
for(i=0;i<fi;i++)

out<<u[i][j]<<" \t"<<v[i][j]<<endl;
out.close();
//求虛擬擴充邊界
double hvix[fi][2],hviy[fi][2],svix[fj][2],sviy[fj][2];
for(i=0;i<fi;i++)
{
	hvix[i][0]=u[i][0];
	hvix[i][1]=u[i][fj-1];
	hviy[i][0]=2*v[i][0]-v[i][1];
	hviy[i][1]=2*v[i][fj-1]-v[i][fj-2];
	
}
for(j=0;j<fj;j++)
{
	svix[j][0]=2*u[0][j]-u[1][j];
	svix[j][1]=2*u[fi-1][j]-u[fi-2][j];
	sviy[j][0]=v[0][j];
	sviy[j][1]=v[fi-1][j];
}

//虛擬網格
double visx[fi+2][fj+2],visy[fi+2][fj+2];
visx[0][0]=svix[0][0];//四個頂點
visy[0][0]=hviy[0][0];
visx[fi+1][0]=svix[fj-1][1];
visy[fi+1][0]=hviy[fi-1][0];

visx[0][fj+1]=svix[fj-1][0];
visy[0][fj+1]=hviy[fi-1][1];
visx[fi+1][fj+1]=svix[fj-1][1];
visy[fi+1][fj+1]=hviy[fi-1][1];

for(i=0;i<fi;i++)
{
	visx[i+1][0]=hvix[i][0];
	visx[i+1][fj+1]=hvix[i][1];
	visy[i+1][0]=hviy[i][0];
	visy[i+1][fj+1]=hviy[i][1];
	
}
for(j=0;j<fj;j++)
{
	visx[0][j+1]=svix[j][0];
	visx[fi+1][j+1]=svix[j][1];
	visy[0][j+1]=sviy[j][0];
	visy[fi+1][j+1]=sviy[j][1];
}

for(j=0;j<fj;j++)
for(i=0;i<fi;i++)
{
	visx[i+1][j+1]=u[i][j];
	visy[i+1][j+1]=v[i][j];
}




out.open("vistual.txt");//將x,y坐標點寫進文件
out<<"title=\"sample mesh\"\nvariables=\"u\",\"v\"\nzone i="<<fi+2<<" j="<<fj+2<<" f=point\n"<<endl;
for(j=0;j<fj+2;j++)
for(i=0;i<fi+2;i++)

out<<visx[i][j]<<" \t"<<visy[i][j]<<endl;
out.close();

//求速度場
double vu[fi+2][fj+2], vut[fi+2][fj+2],vv[fi+2][fj+2],vvt[fi+2][fj+2],en[fi+2][fj+2],rou[fi+2][fj+2],/*密度*/rout[fi+2][fj+2],p[fi+2][fj+2],pt[fi+2][fj+2],tem1,w/*松弛因子*/;
double vu0[fi][fj],vv0[fi][fj],pp[fi][fj],rour[fi][fj];
double a10[3][4][4],b10[3][4][4];
double midda[2][4][4],middb[2][4][4];
double fk[4],gk[4],uk[4];
double fbar[9][4],gbar[9][4],ubar[9][4],ubar0[4];
double middf[5][4],middg[5][4],quf[2][4],qug[2][4];
double nkk[2],n0[2];//任意一點的法相導數數
int num1;//維數
double jj1,jj0/*表示i,j處的J*/,fai1,cit1,rou0,p0,Ma,vc;
double xe, xn, ye, yn;//導數
//邊界速度
double dt,nl;
cout<<"時間步長"<<endl<<"dt=";
cin>>dt;
cout<<"迭代步長數"<<endl<<"nl=";
cin>>nl;
cout<<"初始密度"<<endl<<"rou0=";
cin>>rou0;
cout<<"初始密度"<<endl<<"p0=";
cin>>p0;
cout<<"初始初始Ma數"<<endl<<"Ma=";
cin>>Ma;
cout<<"松弛因子w"<<endl<<"w=";
cin>>w;


tem1=temp(p0,rou0);
cout<<tem1<<endl;
vc=Ma*sqrt(p0*rg/rou0);
cout<<vc;
//賦值初場
for(j=0;j<fj+2;j++)
for(i=0;i<fi+2;i++)
{   
	vu[i][j]=vc;
	vv[i][j]=0;
	rou[i][j]=rou0;
	p[i][j]=p0;
}
//邊界賦值
for(j=0;j<fj+2;j++)
{
	vu[0][j]=vc;
	vv[0][j]=0;
	p[0][j]=p0;
	rou[0][j]=rou0;

}


while(nl>1)
{
	
	for(i=2;i<fi;i++)
	for(j=2;j<fj;j++)
	{
		//out.open("記錄.txt");//記錄寫進文件
//求轉換系數
xe=(visx[i+1][j]-visx[i-1][j])*0.5/di;
xn=(visx[i][j+1]-visx[i][j-1])*0.5/dj;
ye=(visy[i+1][j]-visy[i-1][j])*0.5/di;
yn=(visy[i][j+1]-visy[i][j-1])*0.5/dj;
jj1=jj(xe, xn, ye, yn);
jj0=jj1;
//求轉換向量
num1=0;
jubar(i-1,j+1,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=1;
jubar(i,j+1,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=2;
jubar(i+1,j+1,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=3;
jubar(i-1,j,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=4;
jubar(i,j,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=5;
jubar(i+1,j,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=6;
jubar(i-1,j-1,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=7;
jubar(i,j-1,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);
num1=8;
jubar(i+1,j-1,visx,visy,rou, vu, vv,p,ubar[num1],fbar[num1],gbar[num1]);


/*out<<endl<<endl<<"u  f   g";

for(k=0;k<9;k++)
{
	num1=k;
out<<endl<<ubar[num1][0]<<"  "<<ubar[num1][1]<<"  "<<ubar[num1][2]<<"  "<<ubar[num1][3];
out<<endl<<fbar[num1][0]<<"  "<<fbar[num1][1]<<"  "<<fbar[num1][2]<<"  "<<fbar[num1][3];
out<<endl<<gbar[num1][0]<<"  "<<gbar[num1][1]<<"  "<<gbar[num1][2]<<"  "<<gbar[num1][3]<<endl;
}
out<<endl<<endl;*/

//cout<<vu[i][j]<<"    "<<vv[i][j]<<"   "<<i<<"     "<<j<<endl;

//求A，B值
kcc(i,j,visx, visy, vu[i][j], vv[i][j], p[i][j], rou[i][j], 'a'/*小寫字母*/,a10[0]);
kcc(i-1,j,visx,visy, vu[i-1][j], vv[i-1][j], p[i-1][j], rou[i-1][j], 'a'/*小寫字母*/,a10[1]);
kcc(i+1,j,visx,visy, vu[i+1][j], vv[i+1][j], p[i+1][j], rou[i+1][j], 'a'/*小寫字母*/,a10[2]);

kcc(i,j,visx, visy, vu[i][j], vv[i][j], p[i][j], rou[i][j], 'b'/*小寫字母*/,b10[0]);
kcc(i,j-1,visx, visy, vu[i][j-1], vv[i][j-1], p[i][j-1], rou[i][j-1], 'b'/*小寫字母*/,b10[1]);
kcc(i,j+1,visx,visy,  vu[i][j+1], vv[i][j+1], p[i][j+1],rou[i][j+1],'b'/*小寫字母*/,b10[2]);

/*out<<"   a值"<<endl;
for(km=0;km<4;km++)
for(kn=0;kn<4;kn++)
out<<a10[0][km][kn]<<"       "<<a10[1][km][kn]<<"   "<<a10[2][km][kn]<<endl<<endl;
out<<endl<<endl;
out<<"    b值"<<endl;
for(km=0;km<4;km++)
for(kn=0;kn<4;kn++)
out<<b10[0][km][kn]<<"       "<<b10[1][km][kn]<<"   "<<b10[2][km][kn]<<endl;*/
//求A，B中值

midd(a10[0],a10[1],midda[0]);//A[i-1/2][j]
midd(a10[0],a10[2],midda[1]);//A[i+1/2][j]

/*out<<endl<<"midd值a"<<endl;
for(km=0;km<4;km++)
for(kn=0;kn<4;kn++)
out<<midda[0][km][kn]<<"       "<<midda[1][km][kn]<<endl;*/

midd(b10[0],b10[1],middb[0]);//B[i][j-1/2]
midd(b10[0],b10[2],middb[1]);//B[i][j+1/2]

/*out<<endl<<"midd值b";
for(km=0;km<4;km++)
for(kn=0;kn<4;kn++)
out<<middb[0][km][kn]<<"       "<<middb[1][km][kn]<<endl;*/

//求f,g中值
Fma(fbar[5],fbar[3],middf[0],'-');
FmS(middf[0],middf[0],1/2./di);

Fma(gbar[1],gbar[7],middg[0],'-');
FmS(middg[0],middg[0],1/2./dj);

Fma(fbar[5],fbar[4],middf[1],'-');
FmS(middf[1],middf[1],1/di);

Fma(gbar[1],gbar[2],qug[0],'+');
Fma(gbar[4],qug[0],qug[0],'+');
Fma(gbar[5],qug[0],qug[0],'+');
FmS(qug[0],qug[0],0.25);
Fma(gbar[4],gbar[5],qug[1],'+');
Fma(gbar[7],qug[1],qug[1],'+');
Fma(gbar[8],qug[1],qug[1],'+');
FmS(qug[1],qug[1],0.25);
Fma(qug[0],qug[1],middg[1],'-');
FmS(middg[1],middg[1],1/dj);


Fma(fbar[4],fbar[3],middf[2],'-');
FmS(middf[2],middf[2],1/di);

Fma(gbar[0],gbar[1],qug[0],'+');
Fma(gbar[3],qug[0],qug[0],'+');
Fma(gbar[4],qug[0],qug[0],'+');
FmS(qug[0],qug[0],0.25);
Fma(gbar[3],gbar[4],qug[1],'+');
Fma(gbar[6],qug[1],qug[1],'+');
Fma(gbar[7],qug[1],qug[1],'+');
FmS(qug[1],qug[1],0.25);
Fma(qug[0],qug[1],middg[2],'-');
FmS(middg[2],middg[2],1/dj);

Fma(gbar[2],gbar[4],middg[3],'-');
FmS(middg[3],middg[3],1/dj);

Fma(fbar[1],fbar[2],quf[0],'+');
Fma(fbar[4],quf[0],quf[0],'+');
Fma(fbar[5],quf[0],quf[0],'+');
FmS(quf[0],quf[0],0.25);
Fma(fbar[0],fbar[1],quf[1],'+');
Fma(fbar[3],quf[1],quf[1],'+');
Fma(fbar[4],quf[1],quf[1],'+');
FmS(quf[1],quf[1],0.25);
Fma(quf[0],quf[1],middf[3],'-');
FmS(middf[3],middf[3],1/di);


Fma(gbar[4],gbar[7],middg[4],'-');
FmS(middg[4],middg[4],1/dj);

Fma(fbar[4],fbar[5],quf[0],'+');
Fma(fbar[7],quf[0],quf[0],'+');
Fma(fbar[8],quf[0],quf[0],'+');
FmS(quf[0],quf[0],0.25);
Fma(fbar[3],fbar[4],quf[1],'+');
Fma(fbar[6],quf[1],quf[1],'+');
Fma(fbar[7],quf[1],quf[1],'+');
FmS(quf[1],quf[1],0.25);
Fma(quf[0],quf[1],middf[4],'-');
FmS(middf[4],middf[4],1/di);



/*out<<"midd f 值"<<endl;
for(kn=0;kn<4;kn++)
out<<middf[0][kn]<<"   "<<middf[1][kn]<<"   "<<middf[2][kn]<<"   "<<middf[3][kn]<<"    "<<middf[4][kn]<<endl;
out<<"midd g 值"<<endl;
for(kn=0;kn<4;kn++)
out<<middg[0][kn]<<"   "<<middg[1][kn]<<"   "<<middg[2][kn]<<"   "<<middg[3][kn]<<"    "<<middg[4][kn]<<endl;*/
//求和
Fma(middf[0],middg[0],middf[0],'+');
Fma(middf[1],middg[1],middf[1],'+');
Fma(middf[2],middg[2],middf[2],'+');
Fma(middf[3],middg[3],middf[3],'+');
Fma(middf[4],middg[4],middf[4],'+');

/*out<<"midd f1 和值"<<endl;
for(kn=0;kn<4;kn++)
out<<middf[0][kn]<<"   "<<middf[1][kn]<<"   "<<middf[2][kn]<<"   "<<middf[3][kn]<<"    "<<middf[4][kn]<<endl;*/


FmS(middf[0],middf[0],dt);

FmB(midda[1],middf[1],ubar0);
FmS(ubar0,middf[1],1);

FmB(midda[0],middf[2],ubar0);
FmS(ubar0,middf[2],1);