?? fdtd_3d.cpp
字號:
K_Bx_b[j]*( (Ey_recv_k[i][j] - Ey[i][j][k])*inv_dz -
(Ez[i][j+1][k] - Ez[i][j][k])*inv_dy);
Hx[i][j][k] = K_Hx_a[k]*Hx[i][j][k] +
K_Hx_b[k]*( K_Hx_c[i]*Bx[i][j][k] +
K_Hx_d[i]*Bx_r )/mu_r[Ind[i][j][k]];
}
}
}
j = nlyMIN1;
if ( nly > 0 )
{
#pragma omp for schedule(dynamic,nr_threads) nowait
for (i = 0; i < nlx; i++)
{
for (k = 0; k < nlzMIN1; k++)
{
Bx_r = Bx[i][j][k];
Bx[i][j][k] = K_Bx_a[j]*Bx[i][j][k] +
K_Bx_b[j]*( (Ey[i][j][k+1] - Ey[i][j][k])*inv_dz -
(Ez_recv_j[i][k] - Ez[i][j][k])*inv_dy);
Hx[i][j][k] = K_Hx_a[k]*Hx[i][j][k] +
K_Hx_b[k]*( K_Hx_c[i]*Bx[i][j][k] +
K_Hx_d[i]*Bx_r )/mu_r[Ind[i][j][k]];
}
}
}
#pragma omp for schedule(dynamic,nr_threads)
for (i = 0; i < nlx; i++)
{
for (j = 0; j < nlyMIN1; j++)
{
for (k = 0; k < nlzMIN1; k++)
{
Bx_r = Bx[i][j][k];
Bx[i][j][k] = K_Bx_a[j]*Bx[i][j][k] +
K_Bx_b[j]*( (Ey[i][j][k+1] - Ey[i][j][k])*inv_dz -
(Ez[i][j+1][k] - Ez[i][j][k])*inv_dy);
Hx[i][j][k] = K_Hx_a[k]*Hx[i][j][k] +
K_Hx_b[k]*( K_Hx_c[i]*Bx[i][j][k] +
K_Hx_d[i]*Bx_r )/mu_r[Ind[i][j][k]];
}
}
}
}
////////////////////////////////////////////////////////////////////
//Point source
////////////////////////////////////////////////////////////////////
if (jel_plane_wave == 0 && pt_source_Hx == 1 && n_Coord_ptSource>0 && iter <= switch_off_time)
{
PtSource_J(Hx, time);
}
////////////////////////////////////////////////////////////////////
//Total field scattered field formulation
////////////////////////////////////////////////////////////////////
if (jel_plane_wave == 1 && jel_TS == 1)
{
//the parallelization can be done this way because in case of
//magnetic field there aro no overlapping elements on the interface
#pragma omp parallel sections default(shared)
{
#pragma omp section
{
//face j0
if (jel_TS_planes[2] == 1)
TS_Hx_j0();
}
#pragma omp section
{
//face j1
if (jel_TS_planes[3] == 1)
TS_Hx_j1();
}
#pragma omp section
{
//face k0
if (jel_TS_planes[4] == 1)
TS_Hx_k0();
}
#pragma omp section
{
//face k1
if (jel_TS_planes[5] == 1)
TS_Hx_k1();
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
//Update the Hx field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Update_Hx_send()
{
long i, j, k;
//update Hx_send_j
if ( nly_Hx > 0 )
{
for (i = 0; i < nlx_Hx; i++)
{
for (k = 0; k < nlz_Hx; k++)
{
Hx_send_j[i][k] = Hx[i][nly_HxMIN1][k];
}
}
}
//update Hx_send_k
if ( nlz_Hx > 0 )
{
for (i = 0; i < nlx_Hx; i++)
{
for (j = 0; j < nly_Hx; j++)
{
Hx_send_k[i][j] = Hx[i][j][nlz_HxMIN1];
}
}
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Hy field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Calc_Hy(long nlx, long nly, long nlz)
{
long i, j, k;
double By_r;
long nlxMIN1 = nlx - 1;
long nlzMIN1 = nlz - 1;
#pragma omp parallel default(shared) private(i,j,k,By_r)
{
i = nlxMIN1;
k = nlzMIN1;
if ( nlx > 0 && nlz > 0 )
{
#pragma omp for schedule(dynamic,nr_threads) nowait
for (j = 0; j < nly; j++)
{
By_r = By[i][j][k];
By[i][j][k] = K_By_a[k]*By[i][j][k] +
K_By_b[k]*( (Ez_recv_i[j][k] - Ez[i][j][k])*inv_dx -
(Ex_recv_k[i][j] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] +
K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] +
K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
}
}
i = nlxMIN1;
if ( nlx > 0 )
{
#pragma omp for schedule(dynamic,nr_threads) nowait
for (j = 0; j < nly; j++)
{
for (k = 0; k < nlzMIN1; k++)
{
By_r = By[i][j][k];
By[i][j][k] = K_By_a[k]*By[i][j][k] +
K_By_b[k]*( (Ez_recv_i[j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] +
K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] +
K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
}
}
}
k = nlzMIN1;
if ( nlz > 0 )
{
#pragma omp for schedule(dynamic,nr_threads) nowait
for (i = 0; i < nlxMIN1; i++)
{
for (j = 0; j < nly; j++)
{
By_r = By[i][j][k];
By[i][j][k] = K_By_a[k]*By[i][j][k] +
K_By_b[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex_recv_k[i][j] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] +
K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] +
K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
}
}
}
#pragma omp for schedule(dynamic,nr_threads)
for (i = 0; i < nlxMIN1; i++)
{
for (j = 0; j < nly; j++)
{
for (k = 0; k < nlzMIN1; k++)
{
By_r = By[i][j][k];
By[i][j][k] = K_By_a[k]*By[i][j][k] +
K_By_b[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a[i]*Hy[i][j][k] +
K_Hy_b[i]*( K_Hy_c[j]*By[i][j][k] +
K_Hy_d[j]*By_r )/mu_r[Ind[i][j][k]];
}
}
}
}
////////////////////////////////////////////////////////////////////
//Point source
////////////////////////////////////////////////////////////////////
if (jel_plane_wave == 0 && pt_source_Hy == 1 && n_Coord_ptSource>0 && iter <= switch_off_time)
{
PtSource_J(Hy, time);
}
////////////////////////////////////////////////////////////////////
//Total field scattered field formulation
////////////////////////////////////////////////////////////////////
if (jel_plane_wave == 1 && jel_TS == 1)
{
//the parallelization can be done this way because in case of
//magnetic field there aro no overlapping elements on the interface
#pragma omp parallel sections default(shared)
{
#pragma omp section
{
//i0 face
if(jel_TS_planes[0] == 1)
TS_Hy_i0();
}
#pragma omp section
{
//i1 face
if(jel_TS_planes[1] == 1)
TS_Hy_i1();
}
#pragma omp section
{
//k0 face
if(jel_TS_planes[4] == 1)
TS_Hy_k0();
}
#pragma omp section
{
//k1 face
if(jel_TS_planes[5] == 1)
TS_Hy_k1();
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
//Update the Hy field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Update_Hy_send()
{
long i, j, k;
//update Hy_send_i
if ( nlx_Hy > 0 )
{
for (j = 0; j < nly_Hy; j++)
{
for (k = 0; k < nlz_Hy; k++)
{
Hy_send_i[j][k] = Hy[nlx_HyMIN1][j][k];
}
}
}
//update Hy_send_k
if ( nlz_Hy > 0 )
{
for (i = 0; i < nlx_Hy; i++)
{
for (j = 0; j < nly_Hy; j++)
{
Hy_send_k[i][j] = Hy[i][j][nlz_HyMIN1];
}
}
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Calculate the Hz field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Calc_Hz(long nlx, long nly, long nlz)
{
long i, j, k;
double Bz_r;
long nlxMIN1 = nlx - 1;
long nlyMIN1 = nly - 1;
#pragma omp parallel default(shared) private(i,j,k,Bz_r)
{
i = nlxMIN1;
j = nlyMIN1;
if( nlx > 0 && nly > 0 )
{
#pragma omp for schedule(dynamic,nr_threads) nowait
for (k = 0; k < nlz; k++)
{
Bz_r = Bz[i][j][k];
Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +
K_Bz_b[i]*( (Ex_recv_j[i][k] - Ex[i][j][k])*inv_dy -
(Ey_recv_i[j][k] - Ey[i][j][k])*inv_dx );
Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] +
K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] +
K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
}
}
j = nlyMIN1;
if ( nly > 0 )
{
#pragma omp for schedule(dynamic,nr_threads) nowait
for (i = 0; i < nlxMIN1; i++)
{
for (k = 0; k < nlz; k++)
{
Bz_r = Bz[i][j][k];
Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +
K_Bz_b[i]*( (Ex_recv_j[i][k] - Ex[i][j][k])*inv_dy -
( Ey[i+1][j][k] - Ey[i][j][k])*inv_dx );
Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] +
K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] +
K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
}
}
}
i = nlxMIN1;
if ( nlx > 0 )
{
#pragma omp for schedule(dynamic,nr_threads) nowait
for (j = 0; j < nlyMIN1; j++)
{
for (k = 0; k < nlz; k++)
{
Bz_r = Bz[i][j][k];
Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +
K_Bz_b[i]*( (Ex[i][j+1][k] - Ex[i][j][k])*inv_dy -
(Ey_recv_i[j][k] - Ey[i][j][k])*inv_dx );
Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] +
K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] +
K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
}
}
}
#pragma omp for schedule(dynamic,nr_threads)
for (i = 0; i < nlxMIN1; i++)
{
for (j = 0; j < nlyMIN1; j++)
{
for (k = 0; k < nlz; k++)
{
Bz_r = Bz[i][j][k];
Bz[i][j][k] = K_Bz_a[i]*Bz[i][j][k] +
K_Bz_b[i]*( (Ex[i][j+1][k] - Ex[i][j][k])*inv_dy -
(Ey[i+1][j][k] - Ey[i][j][k])*inv_dx );
Hz[i][j][k] = K_Hz_a[j]*Hz[i][j][k] +
K_Hz_b[j]*( K_Hz_c[k]*Bz[i][j][k] +
K_Hz_d[k]*Bz_r )/mu_r[Ind[i][j][k]];
}
}
}
}
////////////////////////////////////////////////////////////////////
//Point source
////////////////////////////////////////////////////////////////////
if (jel_plane_wave == 0 && pt_source_Hz == 1 && n_Coord_ptSource > 0 && iter <= switch_off_time)
{
PtSource_J(Hz, time);
}
////////////////////////////////////////////////////////////////////
//Total field scattered field formulation
////////////////////////////////////////////////////////////////////
if (jel_plane_wave == 1 && jel_TS == 1)
{
//the parallelization can be done this way because in case of
//magnetic field there aro no overlapping elements on the interface
#pragma omp parallel sections default(shared)
{
#pragma omp section
{
//face i0
if(jel_TS_planes[0] == 1)
TS_Hz_i0();
}
#pragma omp section
{
//face i1
if(jel_TS_planes[1] == 1)
TS_Hz_i1();
}
#pragma omp section
{
//face j0
if(jel_TS_planes[2] == 1)
TS_Hz_j0();
}
#pragma omp section
{
//face j1
if(jel_TS_planes[3] == 1)
TS_Hz_j1();
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
//Update the Hz field components which will be send to other processes
///////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Update_Hz_send()
{
long i, j, k;
long nlx_HzMIN1 = nlx_Hz - 1;
long nly_HzMIN1 = nly_Hz - 1;
//update Hz_send_i
if ( nlx_Hz > 0 )
{
for (j = 0; j < nly_Hz; j++)
{
for (k = 0; k < nlz_Hz; k++)
{
Hz_send_i[j][k] = Hz[nlx_HzMIN1][j][k];
}
}
}
//update Hz_send_j
if ( nly_Hz > 0 )
{
for (i = 0; i < nlx_Hz; i++)
{
for (k = 0; k < nlz_Hz; k++)
{
Hz_send_j[i][k] = Hz[i][nly_HzMIN1][k];
}
}
}
}
void CFDTD_3D::Get_SendRecv(double **&hz_recv_j, double **&hy_recv_k,
double **&ex_send_j, double **&ex_send_k,
double **&hx_recv_k, double **&hz_recv_i,
double **&ey_send_i, double **&ey_send_k,
double **&hy_recv_i, double **&hx_recv_j,
double **&ez_send_i, double **&ez_send_j,
double **&ey_recv_k, double **&ez_recv_j,
double **&hx_send_j, double **&hx_send_k,
double **&ez_recv_i, double **&ex_recv_k,
double **&hy_send_i, double **&hy_send_k,
double **&ex_recv_j, double **&ey_recv_i,
double **&hz_send_i, double **&hz_send_j)
{
//Ex
hz_recv_j = Hz_recv_j; hy_recv_k = Hy_recv_k;
ex_send_j = Ex_send_j; ex_send_k = Ex_send_k;
//Ey
hx_recv_k = Hx_recv_k; hz_recv_i = Hz_recv_i;
ey_send_i = Ey_send_i; ey_send_k = Ey_send_k;
//Ez
hy_recv_i = Hy_recv_i; hx_recv_j = Hx_recv_j;
ez_send_i = Ez_send_i; ez_send_j = Ez_s?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -