C语言绘图流体力学小作业(圆柱绕流层流部分)

//--------------------------------------------------------------------------------------------------
//C语言绘图流体力学小作业(圆柱绕流层流部分)
//题目描述:
//01.参考书目:ISBN 978-7-112-11121-3.
//02.原书采用FORTRAN语言,此处参考原书采用C++语言.
//------------------------------------------------------------------------------------------------------------
//程序编制要求:
//1.计算流函数的值.
//2.绘制流函数分布示意图.
//--------------------------------------------------------------------------------------------------------------------
#include <graphics.h>//包含Easyx模拟TC的BGI绘图库头文件
#include <math.h>//包含数学运算头文件
#include <iostream.h>//
#include <stdio.h>//包含标准输入输出头文件
#define DD 40//宏定义标题边框距离
#define DW 10//宏定义宽度边框距离
#define DH 10//宏定义高度边框距离
#define PI 3.1415926//宏定义圆周率
//宏定义几何参数及计算数据------------------------------------------------------------
#define XRANGE 3.5//X方向边长
#define YRANGE 2.0//Y方向边长
#define RR 1.0//圆柱半径
#define NX 100//
#define E 0.00001//允许误差值
#define HH 0.025//网格高度
//#define ITER 0//迭代次数
#define NI (int(XRANGE/HH+1.1))//X方向总节点数
#define NJ (int(YRANGE/HH+1.1))//Y方向总节点数
#define N1 (int((XRANGE-RR)/HH+0.1))//下底边节点数
#define N2 (int((RR/HH+1.1))//圆柱占用的节点数
#define N3 (NJ-1)//总Y节点数-1
//辅助颜色------------------------------------------
#define c1 RGB(38,47,86)//宏定义辅助线颜色1
#define c2 RGB(38,47,86)//宏定义辅助线颜色2
#define c3 RGB(255,255,255)//宏定义文字及边框颜色
#define cc1 RGB(255,0,0)//宏定义曲线颜色1(红色)
#define cc2 RGB(0,255,0)//宏定义曲线颜色2(绿色)
#define cc3 RGB(0,0,255)//宏定义曲线颜色3(蓝色)
#define cc4 RGB(255,255,0)//宏定义曲线颜色4(黄色)
#define cc5 RGB(0,255,255)//宏定义曲线颜色5(青色)
#define cc6 RGB(0,255,0)//宏定义曲线颜色6(绿色)
#define cc7 RGB(255,0,255)//宏定义曲线颜色7(品红)
#define cc8 RGB(0,0,255)//宏定义曲线颜色8(蓝色)
#define cc9 RGB(255,0,255)//宏定义曲线颜色9(品红)
#define cc10 RGB(0,0,255)//宏定义曲线颜色10(蓝色)
#define d 3//节点标记大小
//窗口全局变量-----------------------------------------------------------------------------------------------------
int W = GetSystemMetrics(SM_CXSCREEN);//获取整个屏幕右下角X坐标
int H = GetSystemMetrics(SM_CYSCREEN);//屏幕Y坐标 
double dxa=(W-2.0*DW)/XRANGE;//实际每个单位对应绘图尺寸(以X方向为主)
//全局变量--------------------------------------------------------------------------------------------------------
int i,j,ITER;//
double A,B,AA,BB,AAA,dx,dy;//
double PSI[NI][NJ];//定义双精度数组(流函数的值)
//函数声明---------------------------------------------------------------------------------------------------------
void sub_frame();//绘制框架子函数
void sub_calculate();//计算数据子函数
void draw();//绘图子函数
//主函数-----------------------------------------------------------------------------------------------------------
void main()
{
	sub_calculate();
	sub_frame();//绘图初始化
	MOUSEMSG m;//定义鼠标消息 
	while(true)//循环
	{
		m=GetMouseMsg();//获取一条鼠标消息 
		switch(m.uMsg)//根据获得的消息选择分支
		{
		case WM_LBUTTONDOWN://鼠标左键单击时判断数据输入
			{
				sub_calculate();
				break;//
			}//结束case(结束鼠标左键单击事件消息处理)
		case WM_RBUTTONDOWN://鼠标移动的时候画个空心的圆
			return;//返回while
		}//结束switch(结束鼠标消息)
	}//结束while*/
}//结束主函数
//绘制框架子函数---------------------------------------------------------------------------------------------------
void sub_frame()
{
	//1.全屏效果---------------------------------------
	initgraph(W,H);//初始化绘图窗口
	HWND hWnd = GetHWnd();//获取窗口句柄
	LONG style = GetWindowLong(hWnd,-16);//获得窗口风格
	style = style & ~WS_CAPTION & ~WS_SIZEBOX; //窗口全屏显示且不可改变大小
	SetWindowLong(hWnd,-16,style);//设置窗口风格
	SetWindowPos(hWnd, NULL,0,0,W,H,SWP_NOZORDER);//改变窗口位置,尺寸和Z序	
	//2.绘制工作区边框----------------------------------
	setcolor(c3);rectangle(DW,DD,W-DW,H-DH);//外边框
	setfont(16,0,"黑体");outtextxy(5,10,"C语言绘图流体力学小作业(圆柱绕流层流部分)");//输出标题
	draw();
}//结束子程序
//计算数据子函数-----------------------------------------------------------------------------------------
void sub_calculate()
{
  //1.所有网格点赋初值---------------------------------------------
  for(i=1;i<=NI;i++)
  {
	  for(j=1;j<=NJ;j++)
	  {
		  PSI[i][j]=0.0;//对每个网格点,初值设为0
	  }
  }
  //2.左边界上的网格点赋初值---------------------------------------
  for(j=2;j<=N3;j++)
  {
	  PSI[1][j]=(j-1)*HH;//左边界上的点
  }
  //3.上边界上的网格点赋初值---------------------------------------
  for(i=1;i<=NI;i++)
  {
	  PSI[i][NJ]=N3*HH;//上边界上的点
  }
  //4.计算内点的值---------------------------------------------
  ITER=1;
  AA=1;
  while(AA>=E)
  {
  ITER=ITER+1;//迭代次数+1
  for(i=2;i<=NI;i++)
  {
	  for(j=2;j<=N3;j++)
	  {
		  BB=PSI[i][j];//存储当前值
		  if(j <= N2))//圆柱下侧 
		  {
			  if(i <= N1)//圆柱左侧
			  {
			     PSI[i][j]=0.25*(PSI[i-1][j]+PSI[i][j-1]+PSI[i+1][j]+PSI[i][j+1]);//左下侧内点
			  }
			  else//右侧与圆柱相邻
			  {
				  A=(j-1)*HH-sqrt(1-(NI-i)*(NI-i)*HH*HH);//系数A
				  B=(NI-i)*HH-sqrt(1-(j-1)*(j-1)*HH*HH);//系数B
				  //cout<<"A="<<A<<" "<<"B="<<B;
				  if((A<=0.0)&&(B<=0.0))//在圆柱体内
					  PSI[i][j]=0.0;
				  else
				  {
				    if(A>=HH) A=HH;
				    if(B>=HH) B=HH;
				    //非正则内点的值
                    PSI[i][j]=HH*((PSI[i-1][j]/HH+PSI[i+1][j]/B)/(B+HH)+(PSI[i][j-1]/A+PSI[i][j+1]/HH)/(A+HH))/(1/A+1/B);
				  }//结束else
			  }//结束else
		  }//结束if
		  else//否则节点在圆柱上侧
		  {
             if (i != NI)
               PSI[i][j]=0.25*(PSI[i-1][j]+PSI[i][j-1]+PSI[i+1][j]+PSI[i][j+1]);//非右边界点
			 else
			 {
				 PSI[i][j]=0.25*(2*PSI[i-1][j]+PSI[i][j-1]+PSI[i][j+1]);//右边界点
			 }
		  }
		  AAA=fabs(PSI[i][j]-BB);//存储前后两次计算差值
          if((i==2)&&(j==2)) AA=AAA;//存储最大差值
		  //if(AAA<=AA) break;
		  AA=AAA;//存储差值
	}//结束for
  }//结束for
  }//结束while
  /*//显示计算数值-------------------------------------------------------
  for(i=1;i<=NI;i++)
  {
	  for(j=1;j<=NJ;j++)
	  {
        cout<<PSI[i][j]<<" ";
	  }
	  cout<<endl;
  }
  cout<<"ITER="<<ITER<<endl;*/  
}//结束子函数
//绘图子函数--------------------------------------------------------------------------
void draw()
{
	for(i=1;i<=NI;i++)
	{
	  for(j=1;j<=NJ;j++)
	  {
		  dx=DW+dxa*((double)i/NI)*3.5;//当前节点在绘图区的位置x
		  dy=-dxa*((double)j/NJ)*2.0+H-DH;//当前节点在绘图区的位置y
		  //line(dx,dy,W/2,H/2);
		  if(PSI[i][j]==0) {}//
		  if((PSI[i][j]>0)&&(PSI[i][j]<=0.2)) 
		  {
			  setcolor(cc1);			  
              line(dx,dy-d,dx,dy+d);//画竖直小刻度线
              line(dx-d,dy,dx+d,dy);//画水平小刻度线

		  }//
		  if((PSI[i][j]>0.2)&&(PSI[i][j]<=0.4)) 
		  {
			  setcolor(cc2);			  
              rectangle(dx-d,dy-d,dx+d,dy+d);//绘制小矩形
		  }//
		  if((PSI[i][j]>0.4)&&(PSI[i][j]<=0.6)) 
		  {
			  setcolor(cc3);			  
              line(dx,dy-d,dx,dy+d);//画竖直小刻度线
              line(dx-d,dy,dx+d,dy);//画水平小刻度线
		  }//
		  if((PSI[i][j]>0.6)&&(PSI[i][j]<=0.8))
		  {
			  setcolor(cc4);			  
              rectangle(dx-d,dy-d,dx+d,dy+d);//绘制小矩形
		  }//
		  if((PSI[i][j]>0.8)&&(PSI[i][j]<=1.0))
		  {
			  setcolor(cc5);			  
              line(dx,dy-d,dx,dy+d);//画竖直小刻度线
              line(dx-d,dy,dx+d,dy);//画水平小刻度线
		  }//
		  if((PSI[i][j]>1.0)&&(PSI[i][j]<=1.2)) 
		  {
			  setcolor(cc6);			  
              rectangle(dx-d,dy-d,dx+d,dy+d);//绘制小矩形
		  }//
		  if((PSI[i][j]>1.2)&&(PSI[i][j]<=1.4))
		  {
			  setcolor(cc7);			  
              line(dx,dy-d,dx,dy+d);//画竖直小刻度线
              line(dx-d,dy,dx+d,dy);//画水平小刻度线
		  }//
          if((PSI[i][j]>1.4)&&(PSI[i][j]<=1.6))
		  {
			  setcolor(cc8);			  
              rectangle(dx-d,dy-d,dx+d,dy+d);//绘制小矩形
		  }//
		  if((PSI[i][j]>1.6)&&(PSI[i][j]<=1.8)) 
		  {
              setcolor(cc9);			  
              line(dx,dy-d,dx,dy+d);//画竖直小刻度线
              line(dx-d,dy,dx+d,dy);//画水平小刻度线
		  }//
		  if((PSI[i][j]>1.8)&&(PSI[i][j]<=2.0)) 
		  {
              setcolor(cc10);			  
              rectangle(dx-d,dy-d,dx+d,dy+d);//绘制小矩形
		  }//
	  }//结束for
	}//结束for
}//结束子函数
//-------------------------------------------------------------------------------------------------------------------