exam1.c

用C++ 写的SIMPLE算法程序,还可以

#include "engine.h"
#include "matrix.h"

#include "global_var.h" 

void user_could_reset_some_setting_here()
{
 char *string="-=temprature=-";
 last=1;       //迭代次数改为10
 lsolve[4]=1;   //求解温度场
 lprint[4]=1;   //打印温度场
 title[4]=strdup(string);  //给出f[i][j][3]的说明字符串
 xl=1.0;        //求解域 x方向宽度＝1.0 
 yl=2.0;        //求解域 y方向宽度＝2.0 
 l1=7;    
 m1=7;          
 }

void user_generate_velocity_grid()  //速度u v的控制网格
{  int i,j;
   double dx,dy;  
   xu[2]=0;       //  第二个主控制容积在x方向的w界面位置，即速度u(i,j)所在位置的x坐标＝0   
                  //  难道说，第一个控制容积留给了虚拟的边壁？ -=???=-
   dx=xl/(double)(l1-2);  // x方向上的网格宽度
   for (i=3;i<=l1;i++)  // l1-1就是最后一个控制体（虚拟边界控制体）的i坐标了。
	                   //计数的习惯与C的数组的序列号保持一致，从0开始。
   {
	   xu[i]=xu[i-1]+dx; //网格在x方向是均匀，所以第三列～第七列控制容积的w面简单以dx递增
   }
   // NOTE : xu[1] 没有被赋值
   yv[2]=0;
   dy=yl/(double)(m1-2);
   for (j=3;j<=m1;j++) 
   {
	   yv[j]=yv[j-1]+dy;
   }
   // NOTE : yv[1] 没有被赋值
}

void user_correct_density()
{}

void user_bound_condition_setup()
{}


void user_output()
{   
	//fprintf(filehandle8,"\niter=%3d  SMAX=%7.3f  SSUM=%9.4f  v(4,7)=%9.4f   T(4,7)=%9.4f   F(4,7,3)=%9.4f\n",iter,smax,ssum,v(4,7),t(4,7),f[4][7][3]);
	fprintf(filehandle8,"\niter==%3d time=%7.3f\n",iter,time);	printinfo();
    user_plot();
}


void user_gamsor()
{   int i,j;
	if (iter>0) return;
	for (j=1;j<=m1;j++)
  	  for (i=1;i<=l1;i++)
	    gam(i,j)=1;  
}

void user_initial_grid_value()
{   

	int i,j;
    for (j=1;j<=m1;j++)
  	  for (i=1;i<=l1;i++)
	  {  
		  t(i,j)=0;     //#define t(i,j)   (f[i][j][4])
          if ((i==1)||(i==l1)||(j==1)||(j==m1) )
			  t(i,j)=(x[i]+y[j]+x[i]*y[j]);
	 
	  }
}

void give_matlab_a_double_matrix(Engine* ep, char * matrix_name, int matrix_x_dimension, int matrix_y_dimension, double * databuffer) 
{   mxArray *T;
	T=mxCreateDoubleMatrix(matrix_x_dimension,matrix_y_dimension,mxREAL);
	mxSetName(T, matrix_name);
	memcpy((void*)mxGetPr(T),(void*)databuffer,matrix_x_dimension*matrix_y_dimension*sizeof(double));
	engPutArray(ep, T);
    mxDestroyArray(T);
}

void give_matlab_a_interger(Engine* ep, char * matrix_name, int data)
{ char matlabsentence[255];
  sprintf(matlabsentence,"%s=%d;",matrix_name,data);
  engEvalString(ep,matlabsentence);
}
	
void user_plot() //MATLAB画图
{  
   
       int i,j,matrixcount,maincv_ibeg,maincv_iend,maincv_jbeg,maincv_jend;
       double matrixX[MAXGRIDSIZE*MAXGRIDSIZE],matrixY[MAXGRIDSIZE*MAXGRIDSIZE],matrixT[MAXGRIDSIZE*MAXGRIDSIZE];
	    static Engine *ep;
      
	  maincv_ibeg=1;
      maincv_iend=l1-1;
      maincv_jbeg=1;
      maincv_jend=m1-1;   

   if (iter==0)
	  {
		  ep = engOpen(NULL);
      
	  }
	 
	    
		matrixcount=0;
		for (j=maincv_jbeg;j<=maincv_jend;j++)
		   for (i=maincv_ibeg;i<=maincv_iend;i++)
		     {
			   matrixX[matrixcount]=x[i];
	           matrixY[matrixcount]=y[j];
			   matrixT[matrixcount]=t(i,j);
			   matrixcount++;
             }
		   give_matlab_a_double_matrix(ep,"X",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixX);
           give_matlab_a_double_matrix(ep,"Y",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixY);
		   give_matlab_a_double_matrix(ep,"T",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixT);
	
		 engEvalString(ep,"pcolor(X,Y,T);colorbar;shading interp;axis equal");
		 
       
}