fdtd_3d_complex.cpp

采用UMPL边界的FDTD编程

#include "FDTD_3D_COMPLEX.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_3D_COMPLEX::CFDTD_3D_COMPLEX(void)
{
	Ex_re = NULL; Ey_re = NULL; Ez_re = NULL;
	Ex_im = NULL; Ey_im = NULL; Ez_im = NULL;

	Param_Source_re = NULL; Param_Source_im = NULL;
	
	#ifndef run_enviroment
		Path_Load_Workspace_Data = NULL;
		Path_Save_Workspace_Data = NULL;
	#endif
}

CFDTD_3D_COMPLEX::~CFDTD_3D_COMPLEX(void)
{
	if(Param_Source_re)
	    Param_Source_re = Free_Matrix_2D<double>(Param_Source_re);
	
	if(Param_Source_im)
	    Param_Source_im = Free_Matrix_2D<double>(Param_Source_im);

	FDTD_3D_re.Free_Mem(); //Free the memory allocated in the FDTD_3D_re object
	FDTD_3D_im.Free_Mem(); //Free the memory allocated in the FDTD_3D_im  object

	#ifndef run_enviroment
		if (Path_Load_Workspace_Data)
		{
			free(Path_Load_Workspace_Data);
			Path_Load_Workspace_Data = NULL;
		}
		
		if (Path_Save_Workspace_Data)
		{
			free(Path_Save_Workspace_Data);
			Path_Save_Workspace_Data = NULL;
		}
	#endif
}


///////////////////////////////////////////////////////////////////////////////////////
//Initialize the real and imaginary parts of the FDTD algorithm
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D_COMPLEX::Init(long ***&Index, long n_x, long n_y, long n_z, double **&Materials, long n_Mat, 
		                   long n_PML_x_1, long n_PML_x_2, long n_PML_y_1, long n_PML_y_2, long n_PML_z_1,
		                   long n_PML_z_2, double d_t, double d_x, double d_y, double d_z,
		                   double PML_eps_r_x1, double PML_mu_r_x1, double PML_eps_r_x2, 
		                   double PML_mu_r_x2, double PML_eps_r_y1, double PML_mu_r_y1, 
		                   double PML_eps_r_y2, double PML_mu_r_y2, double PML_eps_r_z1, 
		                   double PML_mu_r_z1, double PML_eps_r_z2, double PML_mu_r_z2)
{
	pi = 180.0*atan(1.0)/45.0;

	//dimensions of the computational space
    nx = n_x;
	ny = n_y; 
	nz = n_z;

	//the time steep
	dt = d_t;
	dtDIV2 = dt/2.0;

	//the space steep
	dx = d_x;
	dy = d_y;
	dz = d_z;

	nPML_x_1 = n_PML_x_1;
	nPML_x_2 = n_PML_x_2;
	nPML_y_1 = n_PML_y_1;
	nPML_y_2 = n_PML_y_2;
	nPML_z_1 = n_PML_z_1;
	nPML_z_2 = n_PML_z_2;

	if(FDTD_3D_re.Init(Index, n_x, n_y, n_z, Materials, n_Mat, n_PML_x_1, n_PML_x_2, 
	                   n_PML_y_1, n_PML_y_2, n_PML_z_1, n_PML_z_2, dt, dx, dy, dz, 
					   PML_eps_r_x1, PML_mu_r_x1, PML_eps_r_x2, PML_mu_r_x2, 
					   PML_eps_r_y1, PML_mu_r_y1, PML_eps_r_y2, PML_mu_r_y2, 
					   PML_eps_r_z1, PML_mu_r_z1, PML_eps_r_z2, PML_mu_r_z2))
	{
		cout << "Memory allocation problem in FDTD_3D_Complex.Init - FDTD_re" << endl;
		return 1;
	}

	if(FDTD_3D_im.Init(Index, n_x, n_y, n_z, Materials, n_Mat, n_PML_x_1, n_PML_x_2, 
	                   n_PML_y_1, n_PML_y_2, n_PML_z_1, n_PML_z_2, dt, dx, dy, dz, 
					   PML_eps_r_x1, PML_mu_r_x1, PML_eps_r_x2, PML_mu_r_x2, 
					   PML_eps_r_y1, PML_mu_r_y1, PML_eps_r_y2, PML_mu_r_y2, 
					   PML_eps_r_z1, PML_mu_r_z1, PML_eps_r_z2, PML_mu_r_z2))
	{
		cout << "Memory allocation problem in FDTD_3D_Complex.Init - FDTD_im" << endl;
		return 1;
	}

	return 0;
}


///////////////////////////////////////////////////////////////////////////////////////
//To measure the ellapsed time on IBM AIX
///////////////////////////////////////////////////////////////////////////////////////
#ifndef run_enviroment //IBM AIX
int CFDTD_3D_COMPLEX::Init_Mesure_Ellapsed_Time(double TimeStart, double LimitTime, int MyRank,
												int LoadWorkspaceData, char *path_load, char *path_save)
{
	time_start = TimeStart;
	limit_time = LimitTime;

	load_workspace_data = LoadWorkspaceData;

	Path_Load_Workspace_Data =(char *) calloc(512,sizeof(char));
	Path_Save_Workspace_Data =(char *) calloc(512,sizeof(char));
	if (!Path_Load_Workspace_Data || !Path_Save_Workspace_Data)
	{
		cout << "memory allocation problem in FDTD_3D_Complex.Init_Mesure_Ellapsed_Time - Path" << endl;
		return 1;
	}
	
	strcpy(Path_Load_Workspace_Data,path_load);
	strcpy(Path_Save_Workspace_Data,path_save);
    	
	myrank = MyRank;

	return 0;
}
#endif

///////////////////////////////////////////////////////////////////////////////////////
//Bloch type boundary conditions for Ex in j direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::bc_Bloch_Ex_jDir(double ***&Ex_re, double ***Ex_im, 
										double cos_kyORb, double sin_kyORb)
{
	long i, k;

	//j = 0; 
	for (i = 0; i < nx-1; i++)
	{	
		for (k = 0; k < nz; k++)
		{
			Ex_re[i][0][k] =  Ex_re[i][ny-2][k]*cos_kyORb + Ex_im[i][ny-2][k]*sin_kyORb;
			Ex_im[i][0][k] = -Ex_re[i][ny-2][k]*sin_kyORb + Ex_im[i][ny-2][k]*cos_kyORb;
		}
	}

	//j = ny - 1
	for (i = 0; i < nx-1; i++)
	{	
		for (k = 0; k < nz; k++)
		{
			Ex_re[i][ny-1][k] = Ex_re[i][1][k]*cos_kyORb - Ex_im[i][1][k]*sin_kyORb;
			Ex_im[i][ny-1][k] = Ex_re[i][1][k]*sin_kyORb + Ex_im[i][1][k]*cos_kyORb;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Bloch type boundary conditions for Ex in k direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::bc_Bloch_Ex_kDir(double ***&Ex_re, double ***Ex_im, 
										double cos_kzORc, double sin_kzORc)
{
	long i, j;

	//k = 0; 
	for (i = 0; i < nx-1; i++)
	{	
		for (j = 0; j < ny; j++)
		{
			Ex_re[i][j][0] =  Ex_re[i][j][nz-2]*cos_kzORc + Ex_im[i][j][nz-2]*sin_kzORc;
			Ex_im[i][j][0] = -Ex_re[i][j][nz-2]*sin_kzORc + Ex_im[i][j][nz-2]*cos_kzORc;
		}
	}

	//k = nz - 1;
	for (i = 0; i < nx-1; i++)
	{	
		for (j = 0; j < ny; j++)
		{
			Ex_re[i][j][nz-1] = Ex_re[i][j][1]*cos_kzORc - Ex_im[i][j][1]*sin_kzORc;
			Ex_im[i][j][nz-1] = Ex_re[i][j][1]*sin_kzORc + Ex_im[i][j][1]*cos_kzORc;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Bloch type boundary conditions for Ey in i direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::bc_Bloch_Ey_iDir(double ***&Ey_re, double ***&Ey_im,
										double cos_kxORa, double sin_kxORa)
{
	long j, k;

	//i = 0;
	for (j = 0; j < ny-1; j++)
	{	
		for (k = 0; k < nz; k++)
		{
			Ey_re[0][j][k] =  Ey_re[nx-2][j][k]*cos_kxORa + Ey_im[nx-2][j][k]*sin_kxORa;
			Ey_im[0][j][k] = -Ey_re[nx-2][j][k]*sin_kxORa + Ey_im[nx-2][j][k]*cos_kxORa;
		}
	}

	//i = nx-1;
	for (j = 0; j < ny-1; j++)
	{	
		for (k = 0; k < nz; k++)
		{
			Ey_re[nx-1][j][k] = Ey_re[1][j][k]*cos_kxORa - Ey_im[1][j][k]*sin_kxORa;
			Ey_im[nx-1][j][k] = Ey_re[1][j][k]*sin_kxORa + Ey_im[1][j][k]*cos_kxORa;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Bloch type boundary conditions for Ey in k direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::bc_Bloch_Ey_kDir(double ***&Ey_re, double ***&Ey_im,
										double cos_kzORc, double sin_kzORc)
{
	long i, j;

	//k = 0
	for (i = 1; i < nx; i++)
	{	
		for (j = 0; j < ny-1; j++)
		{	
			Ey_re[i][j][0] =  Ey_re[i][j][nz-2]*cos_kzORc + Ey_im[i][j][nz-2]*sin_kzORc;
			Ey_im[i][j][0] = -Ey_re[i][j][nz-2]*sin_kzORc + Ey_im[i][j][nz-2]*cos_kzORc;
		}
	}

	// k = nz-1;
	for (i = 0; i < nx; i++)
	{	
		for (j = 0; j < ny-1; j++)
		{
			Ey_re[i][j][nz-1] = Ey_re[i][j][1]*cos_kzORc - Ey_im[i][j][1]*sin_kzORc;
			Ey_im[i][j][nz-1] = Ey_re[i][j][1]*sin_kzORc + Ey_im[i][j][1]*cos_kzORc;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Bloch type boundary conditions for Ez in i direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::bc_Bloch_Ez_iDir(double ***&Ez_re, double ***&Ez_im,
										double cos_kxORa, double sin_kxORa)
{
	long j, k;

	//i = 0;
	for (j = 0; j < ny; j++) 
	{	
		for (k = 0; k < nz-1; k++)
		{
			Ez_re[0][j][k] =  Ez_re[nx-2][j][k]*cos_kxORa + Ez_im[nx-2][j][k]*sin_kxORa;
			Ez_im[0][j][k] = -Ez_re[nx-2][j][k]*sin_kxORa + Ez_im[nx-2][j][k]*cos_kxORa;
		}
	}

	//i = nx-1;
	for (j = 0; j < ny; j++)
	{	
		for (k = 0; k < nz-1; k++)
		{
			Ez_re[nx-1][j][k] = Ez_re[1][j][k]*cos_kxORa - Ez_im[1][j][k]*sin_kxORa;
			Ez_im[nx-1][j][k] = Ez_re[1][j][k]*sin_kxORa + Ez_im[1][j][k]*cos_kxORa;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Bloch type boundary conditions for Ez in j direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::bc_Bloch_Ez_jDir(double ***&Ez_re, double ***&Ez_im,
										double cos_kyORb, double sin_kyORb)
{
	long i, k;

	//j = 0
	for (i = 0; i < nx; i++)
	{	
		for (k = 0; k < nz-1; k++)
		{
			Ez_re[i][0][k] =  Ez_re[i][ny-2][k]*cos_kyORb + Ez_im[i][ny-2][k]*sin_kyORb;
			Ez_im[i][0][k] = -Ez_re[i][ny-2][k]*sin_kyORb + Ez_im[i][ny-2][k]*cos_kyORb;
		}
	}

	//j = ny-1
	for (i = 0; i < nx; i++)
	{	
		for (k = 0; k < nz-1; k++)
		{
			Ez_re[i][ny-1][k] = Ez_re[i][1][k]*cos_kyORb - Ez_im[i][1][k]*sin_kyORb;
			Ez_im[i][ny-1][k] = Ez_re[i][1][k]*sin_kyORb + Ez_im[i][1][k]*cos_kyORb;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init in case of symmetric or anti-symmetric boundary conditions
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_COMPLEX::Init_Symmetries(int jel_Sym_ia,  int jel_Sym_ib,  int jel_Sym_ja,
	                                   int jel_Sym_jb,  int jel_Sym_ka,  int jel_Sym_kb,
									   int jel_ASym_ia, int jel_ASym_ib, int jel_ASym_ja,
	                                   int jel_ASym_jb, int jel_ASym_ka, int jel_ASym_kb)
{
	FDTD_3D_re.Init_bc_Symmetric(jel_Sym_ia, jel_Sym_ib, jel_Sym_ja,
	                             jel_Sym_jb, jel_Sym_ka, jel_Sym_kb);

	FDTD_3D_re.Init_bc_AntiSymmetric(jel_ASym_ia, jel_ASym_ib, jel_ASym_ja,
	                                 jel_ASym_jb, jel_ASym_ka, jel_ASym_kb);

	FDTD_3D_im.Init_bc_Symmetric(jel_Sym_ia, jel_Sym_ib, jel_Sym_ja,
	                             jel_Sym_jb, jel_Sym_ka, jel_Sym_kb);

	FDTD_3D_im.Init_bc_AntiSymmetric(jel_ASym_ia, jel_ASym_ib, jel_ASym_ja,
	                                 jel_ASym_jb, jel_ASym_ka, jel_ASym_kb);
}


///////////////////////////////////////////////////////////////////////////////////////
//Init the point sources