📄 fdtd_3d_pml_ez.cpp
字号:
#include "fdtd_3d_pml_ez.h"
///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ez field
///////////////////////////////////////////////////////////////////////////////////////
void calc_Ez_1(double ***Ez, double ***Hx, double ***Hy, double ***Fz_1, double ***Gz_1,
long ***Ind, double *K_a, double *K_b, double *K_Gz_a_1, double *K_Gz_b_1,
double *K_Ez_a_1, double *K_Ez_b_1, double *K_Ez_c_1, double *K_Ez_d_1,
double inv_dx, double inv_dy, long nPML_x_1, long nPML_y_1, long nPML_z_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 calc_Ez_2(double ***Ez, double ***Hx, double ***Hy, double ***Fz_2, long ***Ind,
double *K_a, double *K_b, double *K_Ez_a_1, double *K_Ez_b_1,
double *K_Ez_c_1, double *K_Ez_d_1, double inv_dx, double inv_dy,
long nPML_x_1, long nx_MIN_nPML_x_2, long nPML_y_1, long nPML_z_1)
{
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 calc_Ez_3(double ***Ez, double ***Hx, double ***Hy, double ***Fz_3, double ***Gz_3,
long ***Ind, double *K_a, double *K_b, double *K_Gz_a_2, double *K_Gz_b_2,
double *K_Ez_a_1, double *K_Ez_b_1, double *K_Ez_c_1, double *K_Ez_d_1,
double inv_dx, double inv_dy, long nx_MIN_nPML_x_2, long nx_MIN_1,
long nPML_y_1, long nPML_z_1, long nPML_x_2_MIN_1)
{
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 calc_Ez_4(double ***Ez, double ***Hx, double ***Hy, double ***Fz_4, double ***Gz_4,
long ***Ind, double *K_a, double *K_b, double *K_Gz_a_1, double *K_Gz_b_1,
double *K_Ez_c_1, double *K_Ez_d_1, double inv_dx, double inv_dy,
long nPML_x_1, long nPML_y_1, long ny_MIN_nPML_y_2, long nPML_z_1,
double inv_TwoOREp0)
{
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 calc_Ez_5(double ***Ez, double ***Hx, double ***Hy, double ***Fz_5, long ***Ind,
double *K_a, double *K_b, double *K_Ez_c_1, double *K_Ez_d_1,
double inv_dx, double inv_dy, long nPML_x_1, long nx_MIN_nPML_x_2,
long nPML_y_1, long ny_MIN_nPML_y_2, long nPML_z_1, double inv_TwoOREp0)
{
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 calc_Ez_6(double ***Ez, double ***Hx, double ***Hy, double ***Fz_6, double ***Gz_6,
long ***Ind, double *K_a, double *K_b, double *K_Gz_a_2, double *K_Gz_b_2,
double *K_Ez_c_1, double *K_Ez_d_1, double inv_dx, double inv_dy,
long nx_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 inv_TwoOREp0)
{
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 calc_Ez_7(double ***Ez, double ***Hx, double ***Hy, double ***Fz_7, double ***Gz_7,
long ***Ind, double *K_a, double *K_b, double *K_Gz_a_1, double *K_Gz_b_1,
double *K_Ez_a_2, double *K_Ez_b_2, double *K_Ez_c_1, double *K_Ez_d_1,
double inv_dx, double inv_dy, long nPML_x_1, long ny_MIN_nPML_y_2,
long ny_MIN_1, long nPML_z_1, long nPML_y_2_MIN_1)
{
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 calc_Ez_8(double ***Ez, double ***Hx, double ***Hy, double ***Fz_8, long ***Ind,
double *K_a, double *K_b, double *K_Ez_a_2, double *K_Ez_b_2,
double *K_Ez_c_1, double *K_Ez_d_1, double inv_dx, double inv_dy,
long nPML_x_1, long nx_MIN_nPML_x_2, long ny_MIN_nPML_y_2, long ny_MIN_1,
long nPML_z_1, long nPML_y_2_MIN_1)
{
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 calc_Ez_9(double ***Ez, double ***Hx, double ***Hy, double ***Fz_9, double ***Gz_9,
long ***Ind, double *K_a, double *K_b, double *K_Gz_a_2, double *K_Gz_b_2,
double *K_Ez_a_2, double *K_Ez_b_2, double *K_Ez_c_1, double *K_Ez_d_1,
double inv_dx, double inv_dy, long nx_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)
{
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 calc_Ez_10(double ***Ez, double ***Hx, double ***Hy, double ***Fz_10,
double ***Gz_10, long ***Ind, double *K_a, double *K_b, double *K_Gz_a_1,
double *K_Gz_b_1, double *K_Ez_a_1, double *K_Ez_b_1, double inv_dx,
double inv_dy, long nPML_x_1, long nPML_y_1, long nPML_z_1,
long nz_MIN_1_MIN_nPML_z_2, double TwoOREp0)
{
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 calc_Ez_11(double ***Ez, double ***Hx, double ***Hy, double ***Fz_11, long ***Ind,
double *K_a, double *K_b, double *K_Ez_a_1, double *K_Ez_b_1,
double inv_dx, double inv_dy, long nPML_x_1, long nx_MIN_nPML_x_2,
long nPML_y_1, long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2, double TwoOREp0)
{
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 calc_Ez_12(double ***Ez, double ***Hx, double ***Hy, double ***Fz_12,
double ***Gz_12, long ***Ind, double *K_a, double *K_b, double *K_Gz_a_2,
double *K_Gz_b_2, double *K_Ez_a_1, double *K_Ez_b_1, double inv_dx,
double inv_dy, long nx_MIN_nPML_x_2, long nx_MIN_1,
long nPML_y_1, long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2,
long nPML_x_2_MIN_1, double TwoOREp0)
{
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 calc_Ez_13(double ***Ez, double ***Hx, double ***Hy, double ***Fz_13, long ***Ind,
double *K_a, double *K_b, double *K_Gz_a_1, double *K_Gz_b_1,
double inv_dx, double inv_dy, long nPML_x_1, long nPML_y_1,
long ny_MIN_nPML_y_2, long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2)
{
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 calc_Ez_14(double ***Ez, double ***Hx, double ***Hy, long ***Ind, double *K_a,
double *K_b, double inv_dx, double inv_dy, 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 calc_Ez_15(double ***Ez, double ***Hx, double ***Hy, double ***Fz_15, long ***Ind,
double *K_a, double *K_b, double *K_Gz_a_2, double *K_Gz_b_2,
double inv_dx, double inv_dy, long nx_MIN_nPML_x_2, long nx_MIN_1,
long nPML_y_1, long ny_MIN_nPML_y_2, long nPML_z_1,
long nz_MIN_1_MIN_nPML_z_2, long nPML_x_2_MIN_1)
{
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;
for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -