fdtd_3d_pml_ez.cpp

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

#include "FDTD_3D_PML_Ez.h"

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ez field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D_PML::calc_Ez_1()
{
	long  i, j, k;
	double Fz_r, Gz_r;

	//1
	for (i = 1; i < nPML_x_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			for (k = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_1[i][j][k];

				Fz_1[i][j][k] = K_a[Ind[i][j][k]]*Fz_1[i][j][k] +
				                K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                    (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_1[i][j][k];

				Gz_1[i][j][k] = K_Gz_a_1[i]*Gz_1[i][j][k] + 
					            K_Gz_b_1[i]*(Fz_1[i][j][k] - Fz_r);

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*( K_Ez_c_1[k]*Gz_1[i][j][k] + K_Ez_d_1[k]*Gz_r);
			}
		}
	}
}

void CFDTD_3D_PML::calc_Ez_2()
{
	long  i, j, k, ia;
	double Fz_r;

	//2
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			for (k = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_2[ia][j][k];

				Fz_2[ia][j][k] = K_a[Ind[i][j][k]]*Fz_2[ia][j][k] +
				                 K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                     (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*( K_Ez_c_1[k]*Fz_2[ia][j][k] + K_Ez_d_1[k]*Fz_r);
			}
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ez_3()
{
	long  i, j, k, ii, ia;
	double Fz_r, Gz_r;

	//3
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			for (k = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_3[ia][j][k];

				Fz_3[ia][j][k] = K_a[Ind[i][j][k]]*Fz_3[ia][j][k] +
				                 K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                     (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_3[ia][j][k];

				Gz_3[ia][j][k] = K_Gz_a_2[ii]*Gz_3[ia][j][k] + 
					             K_Gz_b_2[ii]*(Fz_3[ia][j][k] - Fz_r);

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*( K_Ez_c_1[k]*Gz_3[ia][j][k] + K_Ez_d_1[k]*Gz_r);
			}
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ez_4()
{
	long  i, j, k, ja;
	double Fz_r, Gz_r;

	//4
	for (i = 1; i < nPML_x_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			for (k = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_4[i][ja][k];

				Fz_4[i][ja][k] = K_a[Ind[i][j][k]]*Fz_4[i][ja][k] +
				                 K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                     (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_4[i][ja][k];

				Gz_4[i][ja][k] = K_Gz_a_1[i]*Gz_4[i][ja][k] + 
					             K_Gz_b_1[i]*(Fz_4[i][ja][k] - Fz_r);

				Ez[i][j][k] += inv_TwoOREp0*( K_Ez_c_1[k]*Gz_4[i][ja][k] + K_Ez_d_1[k]*Gz_r);
			}
			ja++;
		}
	}
}

void CFDTD_3D_PML::calc_Ez_5()
{
	long  i, j, k, ia, ja;
	double Fz_r;
	
	//5
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			for (k = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_5[ia][ja][k];

				Fz_5[ia][ja][k] = K_a[Ind[i][j][k]]*Fz_5[ia][ja][k] +
				                  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

				Ez[i][j][k] += inv_TwoOREp0*( K_Ez_c_1[k]*Fz_5[ia][ja][k] + K_Ez_d_1[k]*Fz_r);
			}
			ja++;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ez_6()
{
	long  i, j, k, ii, ia, ja;
	double Fz_r, Gz_r;

	//6
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_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 = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_6[ia][ja][k];

				Fz_6[ia][ja][k] = K_a[Ind[i][j][k]]*Fz_6[ia][ja][k] +
				                  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_6[ia][ja][k];

				Gz_6[ia][ja][k] = K_Gz_a_2[ii]*Gz_6[ia][ja][k] + 
					              K_Gz_b_2[ii]*(Fz_6[ia][ja][k] - Fz_r);

				Ez[i][j][k] += inv_TwoOREp0*( K_Ez_c_1[k]*Gz_6[ia][ja][k] + K_Ez_d_1[k]*Gz_r);
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ez_7()
{
	long  i, j, k, jj, ja;
	double Fz_r, Gz_r;

	//7
	for (i = 1; 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 = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_7[i][ja][k];

				Fz_7[i][ja][k] = K_a[Ind[i][j][k]]*Fz_7[i][ja][k] +
				                 K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                     (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_7[i][ja][k];

				Gz_7[i][ja][k] = K_Gz_a_1[i]*Gz_7[i][ja][k] + 
					             K_Gz_b_1[i]*(Fz_7[i][ja][k] - Fz_r);

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*( K_Ez_c_1[k]*Gz_7[i][ja][k] + K_Ez_d_1[k]*Gz_r);
			}
			ja++;
			jj--;
		}
	}
}

void CFDTD_3D_PML::calc_Ez_8()
{
	long  i, j, k, jj, ia, ja;
	double Fz_r;

	//8
	ia = 0;
	for (i = nPML_x_1; i < nx_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 = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_8[ia][ja][k];

				Fz_8[ia][ja][k] = K_a[Ind[i][j][k]]*Fz_8[ia][ja][k] +
				                  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*( K_Ez_c_1[k]*Fz_8[ia][ja][k] + K_Ez_d_1[k]*Fz_r);
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ez_9()
{
	long  i, j, k, ii, jj, ia, ja;
	double Fz_r, Gz_r;

	//9
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_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 = 0; k < nPML_z_1; k++)
			{
				Fz_r = Fz_9[ia][ja][k];

				Fz_9[ia][ja][k] = K_a[Ind[i][j][k]]*Fz_9[ia][ja][k] +
				                  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_9[ia][ja][k];

				Gz_9[ia][ja][k] = K_Gz_a_2[ii]*Gz_9[ia][ja][k] + 
					              K_Gz_b_2[ii]*(Fz_9[ia][ja][k] - Fz_r);

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*( K_Ez_c_1[k]*Gz_9[ia][ja][k] + K_Ez_d_1[k]*Gz_r);
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ez_10()
{
	long  i, j, k, ka;
	double Fz_r, Gz_r;

	//10
	for (i = 1; i < nPML_x_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
			{
				Fz_r = Fz_10[i][j][ka];

				Fz_10[i][j][ka] = K_a[Ind[i][j][k]]*Fz_10[i][j][ka] +
				                  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_10[i][j][ka];

				Gz_10[i][j][ka] = K_Gz_a_1[i]*Gz_10[i][j][ka] + 
					              K_Gz_b_1[i]*(Fz_10[i][j][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*TwoOREp0*( Gz_10[i][j][ka] - Gz_r);
				ka++;
			}
		}
	}
}

void CFDTD_3D_PML::calc_Ez_11()
{
	long  i, j, k, ia, ka;
	double Fz_r;

	//11
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
			{
				Fz_r = Fz_11[ia][j][ka];

				Fz_11[ia][j][ka] = K_a[Ind[i][j][k]]*Fz_11[ia][j][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*TwoOREp0*( Fz_11[ia][j][ka] - Fz_r);
				ka++;
			}
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ez_12()
{
	long  i, j, k, ii, ia, ka;
	double Fz_r, Gz_r;

	//12
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
			{
				Fz_r = Fz_12[ia][j][ka];

				Fz_12[ia][j][ka] = K_a[Ind[i][j][k]]*Fz_12[ia][j][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_12[ia][j][ka];

				Gz_12[ia][j][ka] = K_Gz_a_2[ii]*Gz_12[ia][j][ka] + 
					               K_Gz_b_2[ii]*(Fz_12[ia][j][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*TwoOREp0*( Gz_12[ia][j][ka] - Gz_r);
				ka++;
			}
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ez_13()
{
	long  i, j, k, ja, ka;
	double Fz_r;

	//13
	for (i = 1; i < nPML_x_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
			{
				Fz_r = Fz_13[i][ja][ka];

				Fz_13[i][ja][ka] = K_a[Ind[i][j][k]]*Fz_13[i][ja][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Ez[i][j][k] = K_Gz_a_1[i]*Ez[i][j][k] + 
					          K_Gz_b_1[i]*(Fz_13[i][ja][ka] - Fz_r);
				ka++;
			}
			ja++;
		}
	}
}

void CFDTD_3D_PML::calc_Ez_14(long nx_a, long nx_b, long ny_a, long ny_b, long nz_a, long nz_b)
{
	long  i, j, k;
	
	//14 Outside of the PML region
	for (i = nx_a; i < nx_b; i++)
	{
		for (j = ny_a; j < ny_b; j++)
		{	
			for (k = nz_a; k < nz_b; k++)
			{
				Ez[i][j][k] = K_a[Ind[i][j][k]]*Ez[i][j][k] +
				              K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
							                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );
			}
		}
	}
}

void CFDTD_3D_PML::calc_Ez_15()
{
	long  i, j, k, ii, ia, ja, ka;
	double Fz_r;

	//15
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			ka = 0;