📄 cflf.cpp
字号:
//partdam
//四点平均和五点平均不同，前者流速大一些
//此程序大坝下游测边界不对????????????????????????????????
//时间越小抹的越平
#include"CFLF.h"
CFLF::CFLF(long double dtt,long double nn_time,int nngrid,int nsizex,int nsizey,long double ndxy)
{
dt=dtt;
n_time=nn_time;
ngrid=nngrid;
sizex=nsizex;
sizey=nsizey;
dxy=ndxy;
// count=0;
gridx=new long double*[sizex+10];
for(int k=0;k<sizex+10;k++)
gridx[k]=new long double[sizey+10];
gridy=new long double*[sizex+10];
for(k=0;k<sizex+10;k++)
gridy[k]=new long double[sizey+10];
gv1=new long double*[sizex+10];
for(k=0;k<sizex+10;k++)
gv1[k]=new long double[sizey+10];
gv2=new long double*[sizex+10];
for(k=0;k<sizex+10;k++)
gv2[k]=new long double[sizey+10];
gv3=new long double*[sizex+10];
for(k=0;k<sizex+10;k++)
gv3[k]=new long double[sizey+10];
cv1=new long double*[sizex+10];
for(k=0;k<sizex+10;k++)
cv1[k]=new long double[sizey+10];
cv2=new long double*[sizex+10];
for(k=0;k<sizex+10;k++)
cv2[k]=new long double[sizey+10];
cv3=new long double*[sizex+10];
for(k=0;k<sizex+10;k++)
cv3[k]=new long double[sizey+10];
}
CFLF::~CFLF()
{
delete []gridx;//节点X坐标
delete []gridy;
for(int i=0;i<sizex+10;i++)
delete []gv1[i];
delete []gv1;
for(i=0;i<sizex+10;i++)
delete []gv2[i];
delete []gv2;
for(i=0;i<sizex+10;i++)
delete []gv3[i];
delete []gv3;
for(i=0;i<sizex+10;i++)
delete []cv1[i];
delete []cv1;
for(i=0;i<sizex+10;i++)
delete []cv2[i];
delete []cv2;
for(i=0;i<sizex+10;i++)
delete []cv3[i];
delete []cv3;
}
void CFLF::form()
{
int dummy1;
long double gridx1,gridy1;
//先装入结点文件
ifstream inn("CFLF.D");//创建输入流并打开文件
if(!inn)
{
cerr<<"cannt open salrbook file.\n";
exit(1);
}
cout.precision(20);
inn.precision(20);
inn>>dummy1;//读入节点数目
//要先定义行列数
for(int i=0;i<=sizex;i++)//分行列存储
for(int j=0;j<=sizey;j++)
{
inn>>gridx1>>gridy1;
gridx[i][j]=gridx1;
gridy[i][j]=gridy1;
}
//添加界外辅助节点
for(i=sizex+1;i<=sizex+2;i++)
for(int j=0;j<=sizey;j++)
{
gridx[i][j]=10*i;
gridy[i][j]=10*j;
}
for(int j=sizey+1;j<=sizey+2;j++)
for(int i=0;i<=sizex;i++)
{
gridx[i][j]=10*i;
gridy[i][j]=10*j;
}
for(i=0;i<=sizex+2;i++)
for(int j=0;j<=sizey+2;j++)
{
if((gridx[i][j]>=0.0)&&(gridx[i][j]<400.01))
{
gv1[i][j]=10.0;
gv2[i][j]=0.0;
gv3[i][j]=0.0;
}
else
{
gv1[i][j]=10.0;
gv2[i][j]=0.0;
gv3[i][j]=0.0;
}
}
//if ((rand() % 16) == 0)
// {
// wavemap[nmap][(rand() % (MAPY-2))+1][(rand() % (MAPX-2))+1] = 50;
// wavemap[cmap][(rand() % (MAPY-2))+1][(rand() % (MAPX-2))+1] = 50;
// }
gv1[50][50]=-10.0;
gv2[50][50]=0.0;
gv3[50][50]=0.0;
inn.close();
}
long double CFLF::fun1(long double a,long double b)
{
long double y;
y=a*a/b+9.8*b*b/2.0;
return y;
}
//等价于F的表达式的三个维向量（可能有误）
long double CFLF::sv1(long double h1,long double vh1,long double h2,long double vh2,long double dt,long double dxy)
{
long double y;
y=0.5*(h1+h2)+0.25*dt/dxy*(vh1-vh2);
return y;
}
long double CFLF::sv2(long double h1,long double uh1,long double vh1,long double h2,long double uh2,long double vh2,long double dt,long double dxy)
{
long double y;
y=0.5*(uh1+uh2)+0.25*dt/dxy*(uh1*vh1/h1-uh2*vh2/h2);
return y;
}
long double CFLF::sv3(long double h1,long double uh1,long double vh1,long double h2,long double uh2,long double vh2,long double dt,long double dxy)
{
long double y;
y=0.5*(vh1+vh2)+0.25*dt/dxy*(fun1(vh1,h1)-fun1(vh2,h2));
return y;
}
//等价于G的表达式的三个维向量（可能有误）
long double CFLF::sv4(long double h1,long double uh1,long double h2,long double uh2,long double dt,long double dxy)
{
long double y;
y=0.5*(h1+h2)+0.25*dt/dxy*(uh1-uh2);
return y;
}
long double CFLF::sv5(long double h1,long double uh1,long double vh1,long double h2,long double uh2,long double vh2,long double dt,long double dxy)
{
long double y;
y=0.5*(uh1+uh2)+0.25*dt/dxy*(fun1(uh1,h1)-fun1(uh2,h2));
return y;
}
long double CFLF::sv6(long double h1,long double uh1,long double vh1,long double h2,long double uh2,long double vh2,long double dt,long double dxy)
{
long double y;
y=0.5*(vh1+vh2)+0.25*dt/dxy*(uh1*vh1/h1-uh2*vh2/h2);
return y;
}
void CFLF::cvpre()
{
long double tv1,tv2,tv3,tv4,tv5,tv6,tv11,tv21,tv31,tv41,tv51,tv61;
for(int j=0;j<=sizey+2;j++)
for(int i=0;i<=sizex+2;i++)
{
cv1[i][j]=0.0;
cv2[i][j]=0.0;
cv3[i][j]=0.0;
}
long double temp1,temp2,temp3,temp4,temp5,temp6;
long double temp7,temp8,temp9,temp10,temp11,temp12;
for(j=0;j<=sizey;j++)
for(int i=0;i<=sizex;i++)
{
temp1=gv1[i][j];//i,j
temp2=gv2[i][j];
temp3=gv3[i][j];
temp4=gv1[i+1][j];//i+1,j
temp5=gv2[i+1][j];
temp6=gv3[i+1][j];
temp7=gv1[i][j+1];//i,j+1
temp8=gv2[i][j+1];
temp9=gv3[i][j+1];
temp10=gv1[i+1][j+1];//i+1,j+1
temp11=gv2[i+1][j+1];
temp12=gv3[i+1][j+1];
tv1=sv1(temp10,temp12,temp4,temp6,dt,dxy);//对么 ？？？？？？
tv2=sv2(temp10,temp11,temp12,temp4,temp5,temp6,dt,dxy);
tv3=sv3(temp10,temp11,temp12,temp4,temp5,temp6,dt,dxy);
tv11=sv1(temp7,temp9,temp1,temp3,dt,dxy);//对么 ？？？？？？
tv21=sv2(temp7,temp8,temp9,temp1,temp2,temp3,dt,dxy);
tv31=sv3(temp7,temp8,temp9,temp1,temp2,temp3,dt,dxy);
tv4=sv4(temp10,temp11,temp7,temp8,dt,dxy);
tv5=sv5(temp10,temp11,temp12,temp7,temp8,temp9,dt,dxy);
tv6=sv6(temp10,temp11,temp12,temp7,temp8,temp9,dt,dxy);
tv41=sv4(temp4,temp5,temp1,temp2,dt,dxy);
tv51=sv5(temp4,temp5,temp6,temp1,temp2,temp3,dt,dxy);
tv61=sv6(temp4,temp5,temp6,temp1,temp2,temp3,dt,dxy);
cv1[i][j]=0.25*(temp1+temp4+temp7+temp10)+0.5*dt/dxy*(tv1-tv11)+0.5*dt/dxy*(tv4-tv41);
cv2[i][j]=0.25*(temp2+temp5+temp8+temp11)+0.5*dt/dxy*(tv2-tv21)+0.5*dt/dxy*(tv5-tv51);
cv3[i][j]=0.25*(temp3+temp6+temp9+temp12)+0.5*dt/dxy*(tv3-tv31)+0.5*dt/dxy*(tv6-tv61);
}
}
void CFLF::lfgv()
{
long double tv1,tv2,tv3,tv4,tv5,tv6,tv11,tv21,tv31,tv41,tv51,tv61;
long double temp1,temp2,temp3,temp4,temp5,temp6;
long double temp7,temp8,temp9,temp10,temp11,temp12;
for(int j=1;j<sizey;j++)
for(int i=1;i<=sizex;i++)
{
temp1=cv1[i-1][j-1];//i-1/2,j-1/2
temp2=cv2[i-1][j-1];
temp3=cv3[i-1][j-1];
temp4=cv1[i][j-1];//i+1/2,j-1/2
temp5=cv2[i][j-1];
temp6=cv3[i][j-1];
temp7=cv1[i-1][j];//i-1/2,j+1/2
temp8=cv2[i-1][j];
temp9=cv3[i-1][j];
temp10=cv1[i][j];//i+1/2,j+1/2
temp11=cv2[i][j];
temp12=cv3[i][j];
tv1=sv1(temp10,temp12,temp4,temp6,dt,dxy);//对么 ？？？？？？
tv2=sv2(temp10,temp11,temp12,temp4,temp5,temp6,dt,dxy);
tv3=sv3(temp10,temp11,temp12,temp4,temp5,temp6,dt,dxy);
tv11=sv1(temp7,temp9,temp1,temp3,dt,dxy);//对么 ？？？？？？
tv21=sv2(temp7,temp8,temp9,temp1,temp2,temp3,dt,dxy);
tv31=sv3(temp7,temp8,temp9,temp1,temp2,temp3,dt,dxy);
tv4=sv4(temp10,temp11,temp7,temp8,dt,dxy);
tv5=sv5(temp10,temp11,temp12,temp7,temp8,temp9,dt,dxy);
tv6=sv6(temp10,temp11,temp12,temp7,temp8,temp9,dt,dxy);
tv41=sv4(temp4,temp5,temp1,temp2,dt,dxy);
tv51=sv5(temp4,temp5,temp6,temp1,temp2,temp3,dt,dxy);
tv61=sv6(temp4,temp5,temp6,temp1,temp2,temp3,dt,dxy);
gv1[i][j]=0.25*(temp1+temp4+temp7+temp10)+0.5*dt/dxy*(tv1-tv11)+0.5*dt/dxy*(tv4-tv41);
gv2[i][j]=0.25*(temp2+temp5+temp8+temp11)+0.5*dt/dxy*(tv2-tv21)+0.5*dt/dxy*(tv5-tv51);
gv3[i][j]=0.25*(temp3+temp6+temp9+temp12)+0.5*dt/dxy*(tv3-tv31)+0.5*dt/dxy*(tv6-tv61);
}
}
void CFLF::cfgv()
{
long double temp1,temp2,temp3,temp4,temp5,temp6;
long double temp7,temp8,temp9,temp10,temp11,temp12;
for(int j=1;j<sizey;j++)
for(int i=1;i<=sizex;i++)
{
//应如何取周围节点？？？可能不对
temp1=cv1[i-1][j-1];//i-1/2,j-1/2
temp2=cv2[i-1][j-1];
temp3=cv3[i-1][j-1];
temp4=cv1[i][j-1];//i+1/2,j-1/2
temp5=cv2[i][j-1];
temp6=cv3[i][j-1];
temp7=cv1[i-1][j];//i-1/2,j+1/2
temp8=cv2[i-1][j];
temp9=cv3[i-1][j];
temp10=cv1[i][j];//i+1/2,j+1/2
temp11=cv2[i][j];
temp12=cv3[i][j];
gv1[i][j]=gv1[i][j]+0.5*dt/dxy*(temp11+temp5-temp2-temp8+temp12+temp9-temp3-temp6);//疑问？？？
gv2[i][j]=gv2[i][j]+0.5*dt/dxy*(fun1(temp11,temp10)+fun1(temp5,temp4)-fun1(temp2,temp1)-fun1(temp8,temp7)+temp11*temp12/temp10+temp8*temp9/temp7-temp2*temp3/temp1-temp5*temp6/temp4);
gv3[i][j]=gv3[i][j]+0.5*dt/dxy*(temp11*temp12/temp10+temp5*temp6/temp4-temp2*temp3/temp1-temp8*temp9/temp7+fun1(temp12,temp10)+fun1(temp9,temp7)-fun1(temp3,temp1)-fun1(temp6,temp4));
}
}
void CFLF::compute(long double dt,long double& m_time)
{
for(int k=1;k<=1;k++)
{
cvpre();
cfgv();
}
for(int m=1;m<=1;m++)
{
cvpre();
lfgv();
}
m_time=m_time+2*dt;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -