mesh.h

C语言实现二维无粘欧拉方程的求解

#include "iostream.h"
#include "stdio.h"
#include "math.h"
#include "define.h"

void mesh(double coordinate_x[N+4][M+4],double coordinate_y[N+4][M+4],double Area[N+3][M+3],
		  double L_x[N+4][M+3],double L_y[N+4][M+3],double H_x[N+3][M+4],double H_y[N+3][M+4])
{
	double Y0[M+4],YN[M+4];
	double X0,XM;
	double m,n;
	int i,j;

	FILE *fp;

	if((fp=fopen("mesh.dat","w"))==NULL)
		cout<<"can't open file!"<<endl;

	X0=0.0;
	XM=6.28;
	m=double(M-1);
	n=double(N-1);

	for(i=2;i<N+2;i++)
	{
		for(j=2;j<M+2;j++)
		{
			coordinate_x[i][j]=(1-((j-2)/m))*X0+((j-2)/m)*XM;
		}
	}

	for(j=2;j<M+2;j++)
	{
		//YN[j]=0.1*coordinate_x[N][j]+0.1;
		YN[j]=0.178575*cos(coordinate_x[N][j])+0.321425;
		//YN[j]=0.178575*cos(coordinate_x[N][j])+1.0;
		//YN[j]=1.0;
		//YN[j]=1+2.2*(coordinate_x[N][j]-1.5)*(coordinate_x[N][j]-1.5);
		//YN[j]=1.0;
		Y0[j]=-YN[j];
		/*if(coordinate_x[N][j]<=0.0)
			Y0[j]=0.0;
		else 
		if(coordinate_x[N][j]>=0.0&&coordinate_x[N][j]<=0.5)
			//Y0[j]=tan(15.0*pi/180.0)*coordinate_x[N][j];
		else 
		if(coordinate_x[N][j]>0.5)
			Y0[j]=0.13397459;*/
	}

	for(i=2;i<N+2;i++)
	{
		for(j=2;j<M+2;j++)
		{
			coordinate_y[i][j]=(1-((i-2)/n))*Y0[j]+((i-2)/n)*YN[j];
		}
	}
///*
//////////////////////////////////////////////////////////////////////////////////////////////////
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
		{
			coordinate_x[i][j]=coordinate_x[2][2]+(j-2)*(coordinate_x[2][3]-coordinate_x[2][2]);
			coordinate_y[i][j]=coordinate_y[2][2]+(i-2)*(coordinate_y[3][2]-coordinate_y[2][2]);
		}

	for(i=0;i<2;i++)
		for(j=2;j<M+2;j++)
		{
			coordinate_x[i][j]=coordinate_x[2][j];
			coordinate_y[i][j]=coordinate_y[2][j]+(i-2)*(coordinate_y[3][j]-coordinate_y[2][j]);
		}

	for(i=0;i<2;i++)
		for(j=M+2;j<M+4;j++)
		{
			coordinate_x[i][j]=coordinate_x[2][M+1]+(j-M-1)*(coordinate_x[2][M+1]-coordinate_x[2][M]);
			coordinate_y[i][j]=coordinate_y[2][M+1]+(i-2)*(coordinate_y[3][M+1]-coordinate_y[2][M+1]);
		}
////////////////////////////////////////////////////////////////////////////////////////////////
	for(i=2;i<N+2;i++)
		for(j=0;j<2;j++)
		{
			coordinate_x[i][j]=coordinate_x[i][2]+(j-2)*(coordinate_x[i][3]-coordinate_x[i][2]);
			coordinate_y[i][j]=coordinate_y[i][2];
		}

	for(i=2;i<N+2;i++)
		for(j=M+2;j<M+4;j++)
		{
			coordinate_x[i][j]=coordinate_x[i][M+1]+(j-M-1)*(coordinate_x[i][M+1]-coordinate_x[i][M]);
			coordinate_y[i][j]=coordinate_y[i][M+1];
		}
//////////////////////////////////////////////////////////////////////////////////////////////////
	for(i=N+2;i<N+4;i++)
		for(j=0;j<2;j++)
		{
			coordinate_x[i][j]=coordinate_x[N+1][2]+(j-2)*(coordinate_x[N+1][3]-coordinate_x[N+1][2]);
			coordinate_y[i][j]=coordinate_y[N+1][2]+(i-N-1)*(coordinate_y[N+1][2]-coordinate_y[N][2]);
		}

	for(i=N+2;i<N+4;i++)
		for(j=2;j<M+2;j++)
		{
			coordinate_x[i][j]=coordinate_x[N+1][j];
			coordinate_y[i][j]=coordinate_y[N+1][j]+(i-N-1)*(coordinate_y[N+1][j]-coordinate_y[N][j]);
		}

	for(i=N+2;i<N+4;i++)
		for(j=M+2;j<M+4;j++)
		{
			coordinate_x[i][j]=coordinate_x[N+1][M+1]+(j-M-1)*(coordinate_x[N+1][M+1]-coordinate_x[N+1][M]);
			coordinate_y[i][j]=coordinate_y[N+1][M+1]+(i-N-1)*(coordinate_y[N+1][M+1]-coordinate_y[N][M+1]);
		}

/////////////////////////////////////////////////////////////////////////////////////////////////
//*/
		
	for(i=0;i<N+3;i++)
		for(j=0;j<M+3;j++)
		{
			Area[i][j]=0.5*((coordinate_x[i+1][j+1]-coordinate_x[i][j])
						   *(coordinate_y[i+1][j]-coordinate_y[i][j+1])
					       -(coordinate_y[i+1][j+1]-coordinate_y[i][j])
					       *(coordinate_x[i+1][j]-coordinate_x[i][j+1]));
		}
	for(i=0;i<N+3;i++)
		for(j=0;j<M+4;j++)
		{
			H_x[i][j]=coordinate_y[i+1][j]-coordinate_y[i][j];
			H_y[i][j]=-(coordinate_x[i+1][j]-coordinate_x[i][j]);
		}

	for(i=0;i<N+4;i++)
		for(j=0;j<M+3;j++)
		{
			L_x[i][j]=coordinate_y[i][j]-coordinate_y[i][j+1];
		    L_y[i][j]=-(coordinate_x[i][j]-coordinate_x[i][j+1]);
		}

/////////////////////////////////////////////////////////////////////////////////////////////////

	fprintf(fp,"ZONE F=POINT I=%d J=%d\n",M,N);
	
	for(i=2;i<3;i++)
	{
		for(j=2;j<M+2;j++)
		{
			fprintf(fp,"%f %f %f\n",coordinate_x[i][j],coordinate_y[i][j],0);
		}
	}
	fclose(fp);
	
}