fdtd_3d_pml_ex.cpp

利用c++语言写的三维FDTD

#include "fdtd_3d_pml_ex.h"

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ex field - decomposed in 27 blocks
///////////////////////////////////////////////////////////////////////////////////////
void calc_Ex_1(double ***Ex, double ***Hy, double ***Hz, double ***Fx_1, double ***Gx_1,
			   long  ***Ind, double *K_a, double *K_b, double *K_Gx_a_1, double *K_Gx_b_1,
			   double *K_Ex_a_1, double *K_Ex_b_1, double *K_Ex_c_1, double *K_Ex_d_1,
			   double inv_dy, double inv_dz, long  nPML_x_1, long  nPML_y_1, long  nPML_z_1)
{
	double Fx_r, Gx_r;
	long  i, j, k;
	
	//1
	for (i = 0; i < nPML_x_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_1[i][j][k];
				
				Fx_1[i][j][k] = K_a[Ind[i][j][k]]*Fx_1[i][j][k] +
					            K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                    (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_1[i][j][k];

				Gx_1[i][j][k] = K_Gx_a_1[j]*Gx_1[i][j][k] + 
					            K_Gx_b_1[j]*(Fx_1[i][j][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*( K_Ex_c_1[i]*Gx_1[i][j][k] + K_Ex_d_1[i]*Gx_r);
			}
		}
	}
}

void calc_Ex_2(double ***Ex, double ***Hy, double ***Hz, double ***Fx_2, double ***Gx_2,
			   long  ***Ind, double *K_a, double *K_b, double *K_Gx_a_1, double *K_Gx_b_1,
			   double *K_Ex_a_1, double *K_Ex_b_1, double inv_dy, double inv_dz, 
			   long  nPML_x_1, long  nx_MIN_1_MIN_nPML_x_2, long  nPML_y_1, long  nPML_z_1, 
			   double TwoOREp0)
{
	double Fx_r, Gx_r;
	long i, j, k, ia;
	
	//2
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_2[ia][j][k];
				
				Fx_2[ia][j][k] = K_a[Ind[i][j][k]]*Fx_2[ia][j][k] +
					             K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                     (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_2[ia][j][k];

				Gx_2[ia][j][k] = K_Gx_a_1[j]*Gx_2[ia][j][k] + 
					             K_Gx_b_1[j]*(Fx_2[ia][j][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*TwoOREp0*( Gx_2[ia][j][k] - Gx_r);
			}
		}
		ia++;
	}
}

void calc_Ex_3(double ***Ex, double ***Hy, double ***Hz, double ***Fx_3, double ***Gx_3,
			   long  ***Ind, double *K_a, double *K_b, double *K_Gx_a_1, double *K_Gx_b_1,
			   double *K_Ex_a_1, double *K_Ex_b_1, double *K_Ex_c_2, double *K_Ex_d_2,
			   double inv_dy, double inv_dz, long  nx_MIN_1_MIN_nPML_x_2, long  nx_MIN_1, 
			   long  nPML_y_1, long  nPML_z_1, long  nPML_x_2_MIN_1)
{
	double Fx_r, Gx_r;
	long i, j, k, ii, ia;
	
	//3
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_3[ia][j][k];
				
				Fx_3[ia][j][k] = K_a[Ind[i][j][k]]*Fx_3[ia][j][k] +
					             K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                     (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_3[ia][j][k];

				Gx_3[ia][j][k] = K_Gx_a_1[j]*Gx_3[ia][j][k] + 
					             K_Gx_b_1[j]*(Fx_3[ia][j][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*( K_Ex_c_2[ii]*Gx_3[ia][j][k] + K_Ex_d_2[ii]*Gx_r);
			}
		}
		ia++;
		ii--;
	}
}

void calc_Ex_4(double ***Ex, double ***Hy, double ***Hz, double ***Fx_4, long  ***Ind, 
			   double *K_a, double *K_b, double *K_Ex_a_1, double *K_Ex_b_1, 
			   double *K_Ex_c_1, double *K_Ex_d_1, double inv_dy, double inv_dz, 
			   long  nPML_x_1, long  nPML_y_1, long  ny_MIN_nPML_y_2, long  nPML_z_1)
{
	double Fx_r;
	long i, j, k, ja;

	//4
	for (i = 0; i < nPML_x_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_4[i][ja][k];
				
				Fx_4[i][ja][k] = K_a[Ind[i][j][k]]*Fx_4[i][ja][k] +
					             K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                     (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*( K_Ex_c_1[i]*Fx_4[i][ja][k] + K_Ex_d_1[i]*Fx_r);
			}
			ja++;
		}
	}
}


void calc_Ex_5(double ***Ex, double ***Hy, double ***Hz, double ***Fx_5, long  ***Ind, 
			   double *K_a, double *K_b, double *K_Ex_a_1, double *K_Ex_b_1, 
			   double inv_dy, double inv_dz, long  nPML_x_1, long  nx_MIN_1_MIN_nPML_x_2,
			   long  nPML_y_1, long  ny_MIN_nPML_y_2, long  nPML_z_1, double TwoOREp0)
{
	double Fx_r;
	long i, j, k, ia, ja;

	//5
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_5[ia][ja][k];
				
				Fx_5[ia][ja][k] = K_a[Ind[i][j][k]]*Fx_5[ia][ja][k] +
					              K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                      (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*TwoOREp0*( Fx_5[ia][ja][k] - Fx_r);
			}
			ja++;
		}
		ia++;
	}
}


void calc_Ex_6(double ***Ex, double ***Hy, double ***Hz, double ***Fx_6, long  ***Ind, 
			   double *K_a, double *K_b, double *K_Ex_a_1, double *K_Ex_b_1, 
			   double *K_Ex_c_2, double *K_Ex_d_2, double inv_dy, double inv_dz, 
			   long  nx_MIN_1_MIN_nPML_x_2, long  nx_MIN_1, long  nPML_y_1, 
			   long  ny_MIN_nPML_y_2, long  nPML_z_1, long  nPML_x_2_MIN_1)
{
	double Fx_r;
	long i, j, k, ii, ia, ja;

	//6
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_6[ia][ja][k];
				
				Fx_6[ia][ja][k] = K_a[Ind[i][j][k]]*Fx_6[ia][ja][k] +
					              K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                      (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*(K_Ex_c_2[ii]*Fx_6[ia][ja][k] + K_Ex_d_2[ii]*Fx_r);
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void calc_Ex_7(double ***Ex, double ***Hy, double ***Hz, double ***Fx_7, double ***Gx_7,
			   long  ***Ind, double *K_a, double *K_b, double *K_Gx_a_2, double *K_Gx_b_2,
			   double *K_Ex_a_1, double *K_Ex_b_1, double *K_Ex_c_1, double *K_Ex_d_1, 
			   double inv_dy, double inv_dz, long  nPML_x_1, long  ny_MIN_nPML_y_2, 
			   long  ny_MIN_1, long  nPML_z_1, long  nPML_y_2_MIN_1)
{
	double Fx_r, Gx_r;
	long i, j, k, jj, ja;

	//7
	for (i = 0; i < nPML_x_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_7[i][ja][k];
				
				Fx_7[i][ja][k] = K_a[Ind[i][j][k]]*Fx_7[i][ja][k] +
					             K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                     (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_7[i][ja][k];

				Gx_7[i][ja][k] = K_Gx_a_2[jj]*Gx_7[i][ja][k] + 
					             K_Gx_b_2[jj]*(Fx_7[i][ja][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*( K_Ex_c_1[i]*Gx_7[i][ja][k] + K_Ex_d_1[i]*Gx_r);
			}
			ja++;
			jj--;
		}
	}
}


void calc_Ex_8(double ***Ex, double ***Hy, double ***Hz, double ***Fx_8, double ***Gx_8,
			   long  ***Ind, double *K_a, double *K_b, double *K_Gx_a_2, double *K_Gx_b_2,
			   double *K_Ex_a_1, double *K_Ex_b_1, double inv_dy, double inv_dz, 
			   long  nPML_x_1, long  nx_MIN_1_MIN_nPML_x_2, long  ny_MIN_nPML_y_2, 
			   long  ny_MIN_1, long  nPML_z_1, long  nPML_y_2_MIN_1, double TwoOREp0)
{
	double Fx_r, Gx_r;
	long i, j, k, jj, ia, ja;

	//8
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_8[ia][ja][k];
				
				Fx_8[ia][ja][k] = K_a[Ind[i][j][k]]*Fx_8[ia][ja][k] +
					              K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                      (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_8[ia][ja][k];

				Gx_8[ia][ja][k] = K_Gx_a_2[jj]*Gx_8[ia][ja][k] + 
					              K_Gx_b_2[jj]*(Fx_8[ia][ja][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*TwoOREp0*( Gx_8[ia][ja][k] - Gx_r);
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void calc_Ex_9(double ***Ex, double ***Hy, double ***Hz, double ***Fx_9, double ***Gx_9,
			   long  ***Ind, double *K_a, double *K_b, double *K_Gx_a_2, double *K_Gx_b_2,
			   double *K_Ex_a_1, double *K_Ex_b_1, double *K_Ex_c_2, double *K_Ex_d_2, 
			   double inv_dy, double inv_dz, long  nx_MIN_1_MIN_nPML_x_2, long  nx_MIN_1,
			   long  ny_MIN_nPML_y_2, long  ny_MIN_1, long  nPML_z_1, long  nPML_x_2_MIN_1,
			   long  nPML_y_2_MIN_1)
{
	double Fx_r, Gx_r;
	long i, j, k, ii, jj, ia, ja;
	
	//9
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			for (k = 1; k < nPML_z_1; k++)
			{
				Fx_r = Fx_9[ia][ja][k];
				
				Fx_9[ia][ja][k] = K_a[Ind[i][j][k]]*Fx_9[ia][ja][k] +
					              K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                      (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_9[ia][ja][k];

				Gx_9[ia][ja][k] = K_Gx_a_2[jj]*Gx_9[ia][ja][k] + 
					              K_Gx_b_2[jj]*(Fx_9[ia][ja][k] - Fx_r);

				Ex[i][j][k] = K_Ex_a_1[k]*Ex[i][j][k] + 
					          K_Ex_b_1[k]*(K_Ex_c_2[ii]*Gx_9[ia][ja][k] + K_Ex_d_2[ii]*Gx_r);
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}

void calc_Ex_10(double ***Ex, double ***Hy, double ***Hz, double ***Fx_10, 
				double ***Gx_10, long  ***Ind, double *K_a, double *K_b, double *K_Gx_a_1,
				double *K_Gx_b_1, double *K_Ex_c_1, double *K_Ex_d_1, double inv_dy, 
				double inv_dz, long  nPML_x_1, long  nPML_y_1, long  nPML_z_1, 
				long  nz_MIN_nPML_z_2, double inv_TwoOREp0)
{
	double Fx_r, Gx_r;
	long i, j, k, ka;