uplw.cpp

一个一维溃坝模拟动画

// UPLW.cpp: implementation of the UPLW class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "cube.h"
#include "UPLW.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
LWLF::LWLF(long double dtt,long double dxx)
{
	    
     	n=500;
	    //注意不要指向同一快内存
			gv=new long double*[n+3];
	    for(int k=0;k<n+3;k++)
			gv[k]=new long double[2+1];
        	cv=new long double*[n+3];
	    for(k=0;k<n+3;k++)
			cv[k]=new long double[2+1];
		x=new long double[n+2];
		tv1=new long double[n+2];
		tv2=new long double[n+2];
		time=0.0;
		dx=dxx;
		dt=dtt;
		c=dtt/dxx;
		g=9.8;

}
LWLF::~LWLF()
{
	for(int i=0;i<n+3;i++)
		delete []gv[i];
	delete []gv;
	for( i=0;i<n+3;i++)
		delete []cv[i];
	delete []cv;
		delete[]tv1;
		delete[]tv2;
		delete[]x;
}
void LWLF::form()
{
		for(int i=1;i<=n+1;i++)
		{
			x[i]=1.0*dx*(i-1);
		}
		for(int j=1;j<=n+1;j++)
		{
			gv[j][2]=0.0;
			cv[j][2]=0.0;
			if((x[j]>=0.0)&&(x[j]<500.0))
			{
				gv[j][1]=10.0;
			    cv[j][1]=0.0;
			}
			else if((x[j]>=500.0)&&(x[j]<=1000.0))
			{
		    	gv[j][1]=2.0;
		     	cv[j][1]=0.0;
	
			}
			else
			{
				cv[j][1]=0.0;
				gv[j][1]=0.0;
			}
		}


}
long double	LWLF ::fun(long double a)
{
	long double y;
	y=g*a*a/2.0;
	return y;
}
long double	LWLF ::fun1(long double a,long double b)
{
	long double y;
	y=a*a/b+g*b*b/2.0;
	return y;
}
long double LWLF ::sv1(long double h1,long double uh1,long double h2,long double uh2, double dt)
{
	long double y;
	y=0.5*(h1+h2)-0.5*c*(uh1-uh2);
	return y;
}
long double LWLF ::sv2(long double h1,long double uh1,long double h2,long double uh2, double dt)
{
	long double y;
	y=0.5*(uh1+uh2)-0.5*c*(fun1(uh1,h1)-fun1(uh2,h2));
	return y;
}
long double LWLF ::mymax(long double x1,long double x2)
{
	if(x1>x2)
		return x1;
	else 
		return x2;
}
long double LWLF ::mymin(long double x1,long double x2)
{
	if(x1>x2)
		return x2;
	else 
		return x1;
}
void LWLF ::cvpre(long double dt)
	{

    		long double temp1,temp2,temp4,temp5,temp6,temp7;
			long double d1,d2,d3,d4,limited;
				for(int i=2;i<=n-1;i++)
				{

					temp1=0.0;
					temp2=0.0;
					temp1=gv[i][1];
			        temp2=gv[i][2];
			        temp4=gv[i+1][1];
			        temp5=gv[i+1][2];
			        temp6=gv[i-1][1];
			        temp7=gv[i-1][2];
					d1=gv[i+1][1]-gv[i][1];
					d2=gv[i][1]-gv[i-1][1];
					d3=gv[i+1][2]-gv[i][2];
					d4=gv[i][2]-gv[i-1][2];

//					cout<<"gv[i][1]="<<gv[i][1]<<endl;
					limited=mymin(mymin(fabs(d1),fabs(d2)),0.5*mymax(fabs(d1),fabs(d2)));
//					cout<<"limited="<<limited<<endl;
					if((d1*d2)<0)
					{
						if(d1>0)
//								temp1=temp1+limited*d1;
                                temp1=temp1+limited;
      					else if(d1<0)
//								temp1=temp1+limited*d1;
                                temp1=temp1-limited;
//		        				cout<<"temp1="<<temp1<<endl;
						if(fabs(d1)>fabs(d2))
						{
							if(d1>0)
							temp4=temp4-limited;//保存校正后的值
							else if(d1<0)
							temp4=temp4+limited;//保存校正后的值
						}
						else if(fabs(d1)<fabs(d2))
						{
							if(d1>0)
								temp6=temp6-limited;//保存校正后的值
							else if(d1<0)
								temp6=temp6+limited;//保存校正后的值
						}
					}

					limited=mymin(mymin(fabs(d3),fabs(d4)),0.5*mymax(fabs(d3),fabs(d4)));
//					cout<<"d3="<<d3<<" "<<"d4="<<d4<<endl;
					if((d3*d4)<0)
					{
						if(d3>0)
//						temp2=temp2+limited*d3;
						temp2=temp2+limited;
						else if(d3<0)
//						temp2=temp2+limited-d3;
						temp2=temp2-limited;
//						cout<<"temp2="<<temp2<<endl;
///*
					    if(fabs(d3)>fabs(d4))
						{
							if(d3>=0)
								temp5=temp5-limited;//保存校正后的值
							else if(d3<0)
								temp5=temp5+limited;//保存校正后的值
						}
						else if(fabs(d3)<fabs(d4))
						{
							if(d3>=0)
								temp7=temp7-limited;//保存校正后的值
						    else if(d3<0)
								temp7=temp7+limited;//保存校正后的值
						}
					}
//*/
					gv[i][1]=temp1;//为何不保存则有误???是否三个节点均要保存？？结果一样！！！
			        gv[i][2]=temp2;
			        gv[i+1][1]=temp4;
			        gv[i+1][2]=temp5;
			        gv[i-1][1]=temp6;
        	        gv[i-1][2]=temp7;
//					//用校正后的值迭代
//					temp1=gv[i][1];
//			        temp2=gv[i][2];
//			        temp4=gv[i+1][1];
//			        temp5=gv[i+1][2];
//			        temp6=gv[i-1][1];
//			        temp7=gv[i-1][2];
				    tv1[i]=sv1(temp4,temp5,temp1,temp2,dt);
				    tv2[i]=sv2(temp4,temp5,temp1,temp2,dt);
//					cout<<"tv1["<<i<<"]="<<tv1[i]<<"  "<<"tv2["<<i<<"]="<<tv2[i]<<endl;
				}
				tv1[1]=tv1[2];
			    tv2[1]=tv2[2];
				tv1[n]=tv1[n-1];
			    tv2[n]=tv2[n-1];
}
void LWLF ::lf(long double dt)
{
				for(int i=1;i<=n-2;i++)
				{
					cv[i+1][1]=0.5*(tv1[i]+tv1[i+1])-0.5*c*(tv2[i+1]-tv2[i]);
                    cv[i+1][2]=0.5*(tv2[i]+tv2[i+1])-0.5*c*(fun1(tv2[i+1],tv1[i+1])-fun1(tv2[i],tv1[i]));
				}
				cv[1][1]=cv[2][1];
			    cv[1][2]=cv[2][2];
				cv[n][1]=cv[n-1][1];
			    cv[n][2]=cv[n-1][2];
				for(int k=1;k<=n;k++)
				{
					gv[k][1]=0.0;
					gv[k][2]=0.0;
				}
				for( k=1;k<=n;k++)
				{
					gv[k][1]=cv[k][1];
					gv[k][2]=cv[k][2];
				}
				for(k=1;k<=n;k++)
				{
//					cv[k][1]=0.0;
//					cv[k][2]=0.0;
					tv1[k]=0.0;
					tv2[k]=0.0;
				}

}
void LWLF ::lw(long double dt)
{
				for(int i=1;i<=n-2;i++)
				{
					cv[i+1][1]=gv[i+1][1]-c*(tv2[i+1]-tv2[i]);
                    cv[i+1][2]=gv[i+1][2]-c*(fun1(tv2[i+1],tv1[i+1])-fun1(tv2[i],tv1[i]));
				}
				cv[1][1]=cv[2][1];
			    cv[1][2]=cv[2][2];
				cv[n][1]=cv[n-1][1];
			    cv[n][2]=cv[n-1][2];
				for(int k=1;k<=n;k++)
				{
					gv[k][1]=0.0;
					gv[k][2]=0.0;
				}
				for( k=1;k<=n;k++)
				{
					gv[k][1]=cv[k][1];
					gv[k][2]=cv[k][2];
				}
				for(k=1;k<=n;k++)
				{
//					cv[k][1]=0.0;
//					cv[k][2]=0.0;
					tv1[k]=0.0;
					tv2[k]=0.0;
				}

}
void LWLF::caculate()
{
		
		for(int i=1;i<=2;i++)
		{
          cvpre(dt);
			lw(dt);
		    time=time+dt;
		}
	    cvpre(dt);
			lf(dt);
		    time=time+dt;
			gv[n][1]=gv[n-1][1];
			gv[n][2]=gv[n-1][2];
//			cout<<"time="<<time<<endl;
}
/*
void main()
{
	long double dt=0.05,dx=2.0;
	    LWLF test(dt,dx);
		test.form();
		while(test.time<60)
		test.caculate();
	ofstream out("lWLF31.plt");//创建输出流并创建打开文件
	for(int i=1;i<=test.n;i++)//写文件
	{
		out<<setw(10)<<setprecision(5)<<test.x[i]<<" ";
 	    out<<setw(10)<<setprecision(5)<<test.gv[i][1]<<endl;
// 	    out<<setw(10)<<setprecision(5)<<test.gv[i][2]<<endl;
	}
		out<<setw(10)<<setprecision(5)<<1000.0<<"    "<<0.0;
}
*/