// ■円柱周りの流れ
//   右ボタンをクリックすると停止。左ボタンで再開。
//   右ボタンをクリックしながら画面キャプチャすることができる。
#include "myWin.h"
#include "stdio.h"
#include "stdlib.h"
#include "time.h"
#include "math.h"
#define frand()((double)rand()/(RAND_MAX+1))
#define PI 3.14159265358979
#define MX 101
#define MY 101
static double VX[MX][MY], VY[MX][MY], P[MX][MY], Q[MX][MY], TMP[MX][MY];
static double VXN[MX][MY], VYN[MX][MY], X[MX][MY],  Y[MX][MY];
static double DT,RE,DX, DY, EPS, ACCPC;
static int NX0, NY0, NX, NY, KK;
static double rDX, rDY, rDT, rRE, DX2,	DY2, rDX2, rDY2, FCT;
static HDC hBuff; static HBITMAP hBM;     //ビットマップ表示用
double average(double X, double Y){ return 0.5*(X+Y);}
double average(double D1, double D2, double D3, double D4){return 0.25*(D1+D2+D3+D4);}
void drawLine(int X1, int Y1, int X2, int Y2, int W, int Color){//■線分
	HPEN pen=CreatePen(PS_SOLID,W, Color)  ; SelectObject(hBuff, pen)  ;
	MoveToEx(hBuff,X1,Y1,NULL);LineTo(hBuff, X2, Y2); DeleteObject(pen);
}
void drawPolygon(POINT poly[],int N, int lineW, int lColor, int bColor){//■多角形
	HPEN pen=CreatePen(PS_SOLID, lineW, lColor)    ; SelectObject(hBuff, pen)  ;
	HBRUSH brush=CreateSolidBrush(bColor)      ; SelectObject(hBuff, brush);
	Polygon(hBuff,poly,N); DeleteObject(brush) ; DeleteObject(pen);
}

int setX(double X){	return (int)(X*40+400);}
int setY(double Y) { return (int)(Y*40+200);}
void initialize(){
	NX0=40; NY0=40; NX=NX0+1; NY=NY0+1;
	RE = 5000; DT = 0.005; DY = 0.075; KK=100; ACCPC=1; EPS=0.01;
	DX=2*PI/((double)(NX-2)); rDX=1/DX; rDY=1/DY;	rDT=1/DT; rRE=1/RE;
	DX2=rDX*rDX; DY2=rDY*rDY; rDX2=1/(4*DX); rDY2=1/(4*DY); FCT=2/(2*DX2+2*DY2);
	double THT, YY;//座標設定
	for(int j=0;j<NY;j++)for(int i=0;i<NX; i++){
		THT=DX*((double)i-0.5); YY=DY*((double)j-0.5); THT=DX*((double)i-0.5);
		X[i][j]=exp(YY)*cos(THT); Y[i][j]=exp(YY)*sin(THT);
	}
	for(int j=0;j<NY;j++)for(int i=0;i<NX; i++){//圧力と速度の初期設定
		THT=DX*((double)i-1);
		VX[i][j]=-cos(THT); VY[i][j]=sin(THT); P[i][j]=0.0;
	}
}
void tempVelocity(){//■速度設定
	double VYT, VXR, VYA,VYB, VYC, VYD,YY, E1,E2, VYS, VXA,VXB, VXC, VXD,VXS;
	for(int j=1;j<NY0;j++)for(int i=1;i<NX; i++){// VY[][]の計算
		VYT=VY[i][j];
		VXR=average(VX[i-1][j],VX[i][j],VX[i-1][j+1],VX[i][j+1]);
		VYN[i][j]=
			VYT*(VY[i+1][j]-VY[i-1][j])*rDX*0.5 +
			VXR*(VY[i][j+1]-VY[i][j-1])*rDY*0.5 +VY[i][j]*VXR -
			abs(VYT)*(VY[i+1][j]-2*VYT+VY[i-1][j])*rDX*0.5 -
			abs(VXR)*(VY[i][j+1]-2*VYT+VY[i][j-1])*rDY*0.5 ;
	}
	for(int j=2;j<NY-2;j++)for(int i=2;i<NX0; i++){
		VYT=VY[i][j]; 
		VXR= average(VX[i-1][j],VX[i][j],VX[i-1][j+1],VX[i][j+1]);
		VYA=-VY[i+2][j]+8.0*(VY[i+1][j]-VY[i-1][j])+VY[i-2][j];
		VYB= VY[i+2][j]+VY[i-2][j]-4*(VY[i+1][j]+VY[i-1][j]) +6*VYT;
		VYC=-VY[i][j+2]+8.0*(VY[i][j+1]-VY[i][j-1]) +VY[i][j-2];
		VYD= VY[i][j+2]+VY[i][j-2]-4*(VY[i][j+1]+VY[i][j-1]) +6*VYT;
		VYN[i][j]=(VYT*VYA*rDX+abs(VYT)*VYB*rDX
			+VXR*VYC*rDY+abs(VXR)*VYD*rDY)/12.0 +VY[i][j]*VXR;
	}
	for(int j=1;j<NY0;j++)for(int i=1;i<NX; i++){
		YY=DY*((double)j-0.5); E1=exp(-YY); E2=exp(-2.0*YY);
		VYS=(VY[i+1][j]-2*VY[i][j]+VY[i-1][j])*DX2 +
			(VY[i][j+1]-2*VY[i][j]+VY[i][j-1])*DY2-VY[i][j] +
			(VX[i][j]+VX[i][j+1]-VX[i-1][j]-VX[i-1][j+1])*rDX;
		TMP[i][j]=VY[i][j]+DT*(-E1*VYN[i][j]+rRE*E2*VYS);
	}
	for(int j=1;j<NY0;j++)for(int i=1;i<NX; i++)VY[i][j]=TMP[i][j];
	for(int j=1;j<NY;j++)for(int i=1;i<NX0; i++){// VX[][]の計算
		VYT=average(VY[i][j-1],VY[i][j],VY[i+1][j-1],VY[i+1][j]);
		VXR=VX[i][j]; 
		VXN[i][j]=
			VYT*(VX[i+1][j]-VX[i-1][j])*rDX*0.5+
			VXR*(VX[i][j+1]-VX[i][j-1])*rDY*0.5 -VYT*VYT-
			abs(VYT)*(VX[i+1][j]-2*VXR+VX[i-1][j])*rDX*0.5-
			abs(VXR)*(VX[i][j+1]-2*VXR+VX[i][j-1])*rDY*0.5 ;
	}
	for(int j=2;j<NY0;j++)for(int i=2;i<NX-2; i++){
		VYT= average(VY[i][j-1],VY[i][j],VY[i+1][j-1],VY[i+1][j]);
		VXR= VX[i][j]; 
		VXA=-VX[i+2][j]+8*(VX[i+1][j]-VX[i-1][j])+VX[i-2][j];
		VXB= VX[i+2][j]+VX[i-2][j]-4*(VX[i+1][j]+VX[i-1][j]) +6*VXR;
		VXC=-VX[i][j+2]+8*(VX[i][j+1]-VX[i][j-1])+VX[i][j-2];
		VXD= VX[i][j+2]+VX[i][j-2]-4*(VX[i][j+1]+VX[i][j-1]) +6*VXR;
		VXN[i][j]=(VYT*VXA*rDX+abs(VYT)*VXB*rDX+
			VXR*VXC*rDY+abs(VXR)*VXD*rDY)/12.0-VYT*VYT;
	}
	for(int j=1;j<NY;j++)for(int i=1;i<NX0; i++){
		YY=DY*((double)j-1); E1=exp(-YY); E2=exp(-2.0*YY);
		VXS=(VX[i+1][j]-2*VX[i][j]+VX[i-1][j])*DX2 +
			(VX[i][j+1]-2*VX[i][j]+VX[i][j-1])*DY2 -VX[i][j] -
			(VY[i+1][j]+VY[i+1][j-1]-VY[i][j]-VY[i][j-1])*rDX;
		TMP[i][j]=VX[i][j]+DT*(-E1*VXN[i][j]+rRE*E2*VXS);
	}
	for(int j=1;j<NY0;j++)for(int i=1;i<NX; i++)VX[i][j]=TMP[i][j];
}
double solvePoisson(){//■ポアソン方程式の解
	double VXR, DVYDT,DVXDY,ERR, RHS, dRHS;
	for(int j=1;j<NY0;j++)for(int i=1;i<NX0;i++){
		VXR=average(VX[i][j+1],VX[i][j]);
		DVYDT=(VY[i+1][j]-VY[i][j])*rDX;
		DVXDY=(VX[i][j+1]-VX[i][j])*rDY;
		Q[i][j]=exp(DY*((double)j-0.5))*rDT*(DVXDY+VXR+DVYDT);
	}
	FCT=1.0/(2*DX2+2*DY2);
	for(int k=0;k<KK;k++){
		ERR=0;int j=0;
		for(int i=0;i<NX;i++){
			P[i][j]=P[i][j+1]-2*rRE*exp(-(DY*j))*rDY*VX[i][j+2];
			P[i][NY0]=0;
		}
		for(int j=0;j<NY;j++){ P[0][j]=P[NX-2][j]; P[NX0][j]=P[1][j];}
		for(int j=1;j<NY0;j++)for(int i=1;i<NX0;i++){
			RHS=((P[i+1][j]+P[i-1][j])*DX2+
				 (P[i][j+1]+P[i][j-1])*DY2-Q[i][j])*FCT;
			dRHS = RHS-P[i][j]; ERR += dRHS*dRHS;
			P[i][j]=P[i][j]*(1-ACCPC)+RHS*ACCPC;
		}
		if(ERR<0.00001) break;
	}
	return ERR;
}
void newVelocity(){//■圧力による速度変更
	double YY, E1;
	for(int j=1;j<NY0;j++) for(int i=1; i<NX;i++){
		YY=DY*((double)j - 0.5); E1=exp(-YY);
		VY[i][j]+= DT*(-E1*(P[i][j]-P[i-1][j])*rDX);
	}
	for(int j=1;j<NY;j++) for(int i=1; i<NX0;i++){
		YY=DY*((double)j - 1); E1=exp(-YY);
		VX[i][j]+= DT*(-E1*(P[i][j]-P[i][j-1])*rDY);
	}
}
void boundaryCondition(){//■境界条件
	double THT;
	for(int i=0;i<NX;i++){
		VY[i][0]=-VY[i][1];	VX[i][0]=VX[i][2]; VX[i][1]=0;
	}
	for(int i=0;i<NX;i++){
		THT=DX*((double)i-1);
		VY[i][NY0]=sin(THT); VX[i][NY0]=-cos(THT); VX[i][NY]=-cos(THT);
	}
	for(int j=0;j<NY;j++){
		VX[0  ][j] = VX[NX-2][j]; VY[0  ][j] = VY[NX-2][j];
		VX[NX0][j] = VX[1   ][j]; VY[NX0][j] = VY[1   ][j];
		VX[NX ][j] = VX[2   ][j]; VY[NX ][j] = VY[2   ][j];
	}
}
void velocityVector(){//■速度ベクトルの設定
	double THT, VXA, VYA;
	for(int j=0; j<NY; j++)for(int i=0; i<NX; i++){
		THT=DX*((double)i-0.5);
		VXA=average(VX[i][j],VX[i+1][j]);
		VYA=average(VY[i][j],VY[i][j+1]);
		VXN[i][j]=VXA*cos(THT)-VYA*sin(THT);
		VYN[i][j]=VXA*sin(THT)+VYA*cos(THT);
	}
}
void drawGrid(){
	for(int j=0;j<NY0; j++)for(int i=0;i<NX0; i++){
		drawLine(setX(X[i  ][j  ]), setY(Y[i  ][j  ]),
			     setX(X[i+1][j  ]), setY(Y[i+1][j  ]),1,0xFF0000);
		drawLine(setX(X[i  ][j  ]), setY(Y[i  ][j  ]),
			     setX(X[i  ][j+1]), setY(Y[i  ][j+1]),1,0xFF);
	}
	for(int i=0;i<NX0; i++){
		drawLine(setX(X[i  ][NY0]), setY(Y[i  ][NY0]),
			     setX(X[i+1][NY0]), setY(Y[i+1][NY0]),1,0xFF0000);
	}
}
double maxValue(double P[MX][MY], int nx, int ny){
	double MV=-999;
	for(int i=1; i<nx-1; i++) for(int j=1;j<ny-1; j++) 
		if(P[i][j]>MV)	MV=P[i][j];
	return MV;
}
double minValue(double P[MX][MY], int nx, int ny){
	double MV=1E+64;
	for(int i=1; i<nx-1; i++) for(int j=1;j<ny-1; j++) 
		if(P[i][j]<MV)	MV=P[i][j];
	return MV;
}
void exec(){
	for(int k=0;k<KK;k++){
		tempVelocity();	solvePoisson();
		newVelocity() ;	boundaryCondition();
	}
	velocityVector();
}
int color(double V, double minV, double maxV){  //■カラー設定
	double DV=maxV-minV, T=0.5;	if(DV>0.00001) T=(V-minV)/DV;
	int R=(int)(T*512);
	if(R>=512) return RGB(255      , 0         , 0  );
	if(R>=384) return RGB(255      , 511- R * 2, 0  );
	if(R>=256) return RGB((R-256)*2, 255       , 0  );
	if(R<0)    return RGB( 0       , 255       , 0  );
	return            RGB( 0       ,   R        , 0 );	
}
void crearBitmap(){                                         //■ビットマップクリア
	SelectObject(hBuff,GetStockObject(NULL_PEN));
	PatBlt(hBuff,0, 0,1000,400,WHITENESS);
}
void drawP(){
	crearBitmap();
	POINT poly[4];double minV=minValue(P, NX, NY), maxV=maxValue(P, NX,NY);
	for(int j=1;j<NY0; j++)for(int i=1;i<NX0; i++){
		poly[0].x=setX(X[i  ][j  ]);poly[0].y=setY(Y[i  ][j  ]);
		poly[1].x=setX(X[i  ][j+1]);poly[1].y=setY(Y[i  ][j+1]);
		poly[2].x=setX(X[i+1][j+1]);poly[2].y=setY(Y[i+1][j+1]);
		poly[3].x=setX(X[i+1][j  ]);poly[3].y=setY(Y[i+1][j  ]);
		int C=color(P[i][j],minV,maxV); drawPolygon(poly,4, 1, C, C);
	}
	double CX, CY;
	for(int j=1;j<NY0; j++)for(int i=1;i<NX0; i++){
		CX=average(X[i][j],X[i][j+1],X[i+1][j],X[i+1][j+1]);
		CY=average(Y[i][j],Y[i][j+1],Y[i+1][j],Y[i+1][j+1]);
//		CX=X[i][j];CY=Y[i][j];
		int X1=setX(CX), Y1=setY(CY), X2=setX(CX+VXN[i][j]*1), Y2=setY(CY+VYN[i][j]*1);
		drawLine(X1,Y1, X2,Y2,1,0xCC0000);
	}
}
void procTimer(HWND hw, WPARAM wp,LPARAM lp){               //■タイマー処理
    exec(); drawP();
	InvalidateRect(hw,NULL,TRUE);
}
void procLButtonDown(HWND hw, WPARAM wp,LPARAM lp){         //■左ボタン処理
 	SetTimer(hw,1,20,NULL);  //画面キャプチャするときはここをコメント*/
}
void procRButtonDown(HWND hw, WPARAM wp,LPARAM lp){         //■左ボタン処理
    KillTimer(hw,1);
	exec(); drawP();
	InvalidateRect(hw,NULL,TRUE);
}
void initBitmap(HWND hw){                                   //■ビットマップ初期設定
	HDC hdc=GetDC(hw); hBM=CreateCompatibleBitmap(hdc,1000,400);
	hBuff=CreateCompatibleDC(hdc); SelectObject(hBuff,hBM);
	crearBitmap(); ReleaseDC(hw,hdc);
}
void procCreate(HWND hw, WPARAM wp,LPARAM lp){              //■Create時処理
	SetWindowText(hw,TEXT("円柱周りの粘性のある流れ（３次精度上流差分法）"));
	MoveWindow(hw,0,0,500,450,TRUE);
	initBitmap(hw);
	initialize();drawGrid();InvalidateRect(hw,NULL,TRUE);
	procLButtonDown(hw,wp,lp);
}
void procPaint(HWND hw, WPARAM wp,LPARAM lp){               //■Paint時処理
	PAINTSTRUCT ps;	HDC hdc=BeginPaint(hw, &ps);
	BitBlt(hdc, 0,0,1000,400,hBuff,0,0,SRCCOPY);
	EndPaint(hw,&ps);
}
LRESULT CALLBACK WndProc(HWND hw, UINT msg, WPARAM wp,LPARAM lp){//■Window Procedure
	switch(msg){
	case WM_DESTROY       : PostQuitMessage(0)         ; return 0;
	case WM_CREATE        : procCreate(hw,wp,lp)       ; return 0;
	case WM_PAINT         : procPaint(hw,wp,lp)        ; return 0;
	case WM_TIMER         : procTimer(hw,wp,lp)        ; return 0;
	case WM_LBUTTONDOWN   : procLButtonDown(hw,wp,lp)  ; break;
	case WM_RBUTTONDOWN   : procRButtonDown(hw,wp,lp)  ; break;
	}
	return DefWindowProc(hw,msg,wp,lp);
}