📄 exam6.c
字号:
#include "engine.h"
#include "matrix.h"
#include "global_var.h"
void user_could_reset_some_setting_here()
{
int i;
char *string1="-=temprature=-";
char *string2="-=vel u=-";
char *string3="-=str fn=-";
char *string4="-=vel v=-";
char *string5="-=pressure=-";
char *string6="-=dens=-";
title[1]=strdup(string2);
title[2]=strdup(string4);
title[3]=strdup(string3);
title[4]=strdup(string1);
title[11]=strdup(string5);
title[12]=strdup(string6);
relax[1]=0.5;
relax[2]=0.5;
relax[11]=0.8;
for (i=1;i<=4;i++)
{
lsolve[i]=1;
lprint[i]=1;
}
lprint[11]=1;
lprint[12]=1;
lprintdebuginfo=0;
mode=1;
last=65;
xl=0.5;
yl=2.0;
l1=28;
m1=28;
pr=0.7;
amu=1.0;
amup=amu/pr;
tref=300;
rhoref=1;
rhot=rhoref*tref;
dt=0.0001;
for (i=1;i<=nfmax;i++) ntimes[i]=30;
}
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_initial_grid_value()
{
int i,j;
double tin=500.0,tw=300.0,vin=100.0;
double vout;
//vout=vin*xcv[l2]/x[l1]*tw/tin;
int inletgridindex;
double inletwidth=0;
inletgridindex=l1/8;
for (i=l2;i>=(l2-inletgridindex);i--)
{
inletwidth+=xcv[i];
}
vout=vin*inletwidth/x[l1]*tw/tin;
for (j=1;j<=m1;j++)
for (i=1;i<=l1;i++)
{
u(i,j)=0;
v(i,j)=vout;
v(1,j)=0;
v(i,2)=0;
t(i,j)=tw;
}
// v(l2,2)=vin;
// t(l2,1)=tin;
for (i=l2;i>=(l2-inletgridindex);i--)
{
v(i,2)=vin;
t(i,1)=tin;
}
}
void user_correct_density()
{
int i,j;
for (j=1;j<=m1;j++)
for (i=1;i<=l1;i++)
rho(i,j)=rhot/t(i,j);
}
void user_bound_condition_setup()
{
static double flowin;
double factor,fl;
int i;
//if (iter==0) flowin=rho(l2,1)*v(l2,2)*xcv[l2];
int inletgridindex;
double inletwidth=0;
inletgridindex=l1/8;
if (iter1==0) //每一个时间片开始后，重新计算入口流量。因为入口条件可能会变
{ flowin=0;
for (i=l2;i>=(l2-inletgridindex);i--)
{
flowin+=rho(i,1)*v(i,2)*xcv[i];
}
}
fl=0;
for (i=2;i<=l2;i++)
fl=fl+rho(i,m1)*v(i,m2)*xcv[i];
factor=flowin/fl;
for (i=2;i<=l2;i++)
{
v(i,m1)=v(i,m2)*factor;
t(i,m1)=t(i,m2);
}
}
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;
for (j=1;j<=m1;j++)
for (i=1;i<=l1;i++)
{
gam(i,j)=amu;
if (nf==4) gam(i,j)=amup;
if (nf!=1) gam(l1,j)=0;
gam(i,m1)=0;
}
}
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,velocv_ibeg,velocv_iend,velocv_jbeg,velocv_jend,maincv_ibeg,maincv_iend,maincv_jbeg,maincv_jend;
double matrixX[MAXGRIDSIZE*MAXGRIDSIZE],matrixY[MAXGRIDSIZE*MAXGRIDSIZE],matrixU[MAXGRIDSIZE*MAXGRIDSIZE];
double matrixV[MAXGRIDSIZE*MAXGRIDSIZE],matrixP[MAXGRIDSIZE*MAXGRIDSIZE],matrixT[MAXGRIDSIZE*MAXGRIDSIZE],matrixSFN[MAXGRIDSIZE*MAXGRIDSIZE];
static Engine *ep;
velocv_ibeg=2;
velocv_iend=l1-1;
velocv_jbeg=2;
velocv_jend=m1-1;
maincv_ibeg=1;
maincv_iend=l1-1;
maincv_jbeg=1;
maincv_jend=m1-1;
if (iter==0)
{
ep = engOpen(NULL);
give_matlab_a_interger(ep,"last",last);
give_matlab_a_interger(ep,"iter",iter);
engEvalString(ep,"M=moviein(last+1);iter=1");
}
matrixcount=0;
for (j=velocv_jbeg;j<=velocv_jend;j++)
for (i=velocv_ibeg;i<=velocv_iend;i++)
{
matrixX[matrixcount]=xu[i];
matrixY[matrixcount]=yv[j];
matrixU[matrixcount]=u(i,j);
matrixV[matrixcount]=v(i,j);
matrixcount++;
}
give_matlab_a_double_matrix(ep,"VELOX",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixX);
give_matlab_a_double_matrix(ep,"VELOY",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixY);
give_matlab_a_double_matrix(ep,"U",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixU);
give_matlab_a_double_matrix(ep,"V",velocv_iend-velocv_ibeg+1,velocv_jend-velocv_jbeg+1,matrixV);
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);
matrixP[matrixcount]=p(i,j);
matrixSFN[matrixcount]=pc(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);
give_matlab_a_double_matrix(ep,"P",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixP);
give_matlab_a_double_matrix(ep,"SFN",maincv_iend-maincv_ibeg+1,maincv_jend-maincv_jbeg+1,matrixSFN);
engEvalString(ep,"pcolor(X,Y,T); shading interp;caxis([300 500]);colorbar;axis equal");
engEvalString(ep,"M(:,iter)=getframe;iter=iter+1;");
if (iter==last)
{
engEvalString(ep,"figure;pcolor(X,Y,T);shading interp;caxis([300 500]);colorbar;axis equal;hold on;quiver(VELOX,VELOY,U,V);");
engEvalString(ep,"map=colormap;mpgwrite(M,map,'exam6.mpg');");
engEvalString(ep,"figure;contour(X,Y,SFN);colorbar;axis equal;");
}
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -