From beaa522475f2874794fd4166c0c08cb8751d1ec8 Mon Sep 17 00:00:00 2001 From: maksimKravchenko1 Date: Sun, 14 May 2023 18:32:31 +0300 Subject: [PATCH] . --- CMakeLists.txt | 2 + README.md | Bin 2 -> 1334 bytes Wz.txt | 220 ++++ dist/index.js | 20 +- dist/test-addon.js | 140 +-- index.d.ts | 69 +- k_Fz_1_new3.txt | 2 + refactored_pml_tf_sf.m | 902 ++++++++++++++++ src/fdtd/1d-pml/fdtd-pml-1d.h | 1 + src/fdtd/2d-pml copy/fdtd-pml-2d.cpp | 187 ++++ src/fdtd/2d-pml copy/fdtd-pml-2d.h | 126 +++ src/fdtd/2d-pml/fdtd-pml-2d.cpp | 374 +++---- src/fdtd/2d-pml/fdtd-pml-2d.h | 991 ++++++++++++++++-- src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.cpp | 909 ++++++++++++++++ src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.h | 320 ++++++ src/index.cpp | 3 + src/transform-to-js/fdtd-2d copy/fdtd-2d.cpp | 323 ++++++ src/transform-to-js/fdtd-2d copy/fdtd-2d.h | 20 + .../fdtd-2d-tf-sf/fdtd-2d-tf-sf.cpp | 352 +++++++ .../fdtd-2d-tf-sf/fdtd-2d-tf-sf.h | 20 + src/transform-to-js/fdtd-2d/fdtd-2d.cpp | 142 ++- test-addon.ts | 24 +- upml.cpp | 0 upml_tf_sf.m | 983 +++++++++++++++++ 24 files changed, 5714 insertions(+), 416 deletions(-) create mode 100644 Wz.txt create mode 100644 k_Fz_1_new3.txt create mode 100644 refactored_pml_tf_sf.m create mode 100644 src/fdtd/2d-pml copy/fdtd-pml-2d.cpp create mode 100644 src/fdtd/2d-pml copy/fdtd-pml-2d.h create mode 100644 src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.cpp create mode 100644 src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.h create mode 100644 src/transform-to-js/fdtd-2d copy/fdtd-2d.cpp create mode 100644 src/transform-to-js/fdtd-2d copy/fdtd-2d.h create mode 100644 src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.cpp create mode 100644 src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.h create mode 100644 upml.cpp create mode 100644 upml_tf_sf.m diff --git a/CMakeLists.txt b/CMakeLists.txt index a89f97b..fdade0f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,8 +16,10 @@ include_directories(${CMAKE_JS_INC}) file(GLOB SOURCE_FILES "src/*.h" "src/*.cpp" "src/fdtd/2d-pml/*.h" "src/fdtd/2d-pml/*.cpp" "src/fdtd/1d-pml/*.h" "src/fdtd/1d-pml/*.cpp" + "src/fdtd/2d-upml-tf-sf/*.h" "src/fdtd/2d-upml-tf-sf/*.cpp" "src/transform-to-js/fdtd-1d/*.h" "src/transform-to-js/fdtd-1d/*.cpp" "src/transform-to-js/fdtd-2d/*.h" "src/transform-to-js/fdtd-2d/*.cpp" + "src/transform-to-js/fdtd-2d-tf-sf/*.h" "src/transform-to-js/fdtd-2d-tf-sf/*.cpp" ) # This line will tell CMake that we're building a shared library diff --git a/README.md b/README.md index 46b134b197f35e75e0784bedbf94a8dd124693b1..98af9f2339e671f77b4dc73636734b4b76b13665 100644 GIT binary patch literal 1334 zcmZ{k-7W+{5QXd7Q*>Oj;Q?Goga{HB?q#gOFh6T{SUfz=SEbWSRys2~UETF}&Z+(S ze6q?GT9qYhY^nE)-nHH8YqCM_ZS9|Bt+hV%V1Mx6mAqc{wQW1HmPiM-uYISJ%6gCa z*CW|V*;=in+@tbP!M&czBRB5XsEYM`mv&+2^1zpfgQ}XeC(FBje^*M3d)bm=+-cP| z*X^11XY$BB&)qxOjc2O7+E%tvzT^%|AhB-cRT`^WOi# z2WIg?h4hO?vM#W6uFZuq!8dt=6Le1v*y)4X!5Uq|J~)qViOE0#$xija?7GHk$At=( z${aJwgwVxiRR*6_j_w9|k)8X4V_o#LS+z6U(mHgHV3(d2jsiW99rj{ITc3wKVc1b` zOn>jJau&UUWaw+u319htBqktNSn_4GLn{##ly@~HbjI8S4(Oi=j9Cjj|7RzB=_4?C KRjgXf% x) && (i < (x + x_cells))... + && (j > y) && (j < (y + y_cells))) + Index(i,j) = 1; + end + if ((i+0.5 > x) && (i+0.5 < (x + x_cells))... + && (j > y) && (j < (y + y_cells))) + IndexX(i,j) = 1; + end + if ((i > x) && (i < (x + x_cells))... + && (j+0.5 > y) && (j+0.5 < (y + y_cells))) + IndexY(i,j) = 1; + end + end + end +end + + + + + +for ii = 0:(upml_width-1) + sigma_x = sigma_max*((upml_width - ii)/ upml_width).^m; + ka_x = 1 + (ka_max - 1)*((upml_width - ii)/ upml_width).^m; + + k_Fz_a(end-ii) = (2*epsilon_0*ka_x - sigma_x*dt)./... + (2*epsilon_0*ka_x + sigma_x*dt); + + k_Fz_b(end-ii) = 2*epsilon_0*dt./(2*epsilon_0*ka_x+sigma_x*dt)/dx; + k_Hz_a(end-ii) = (2*epsilon_0*ka_x - sigma_x*dt)./(2*epsilon_0*ka_x + sigma_x*dt); + k_Hz_b(end-ii) = 2*epsilon_0*dt./(2*epsilon_0*ka_x+sigma_x*dt)/(mu_0*Material(1,2)*dx); + + sigma_x = sigma_max*((upml_width - ii - 0.5)/ upml_width).^m; + ka_x = 1 + (ka_max - 1)*((upml_width - ii - 0.5)/ upml_width).^m; + + k_Hy_a(end-ii) = (2*epsilon_0*ka_x - sigma_x*dt)./(2*epsilon_0*ka_x + sigma_x*dt); + k_Hy_b(end-ii) = 2*epsilon_0*dt./(2*epsilon_0*ka_x + sigma_x*dt)/(mu_0*Material(1,2)*dx); + k_Ey_a(end-ii) = (2*epsilon_0*ka_x - sigma_x*dt)./(2*epsilon_0*ka_x + sigma_x*dt); + k_Ey_b(end-ii) = 2*dt./(2*epsilon_0*ka_x + sigma_x*dt)/(Material(1,1)*dx); +end + +%% Another transformation coefficients +for ii = mmm:number_of_materials + K_a(ii) = (2*epsilon_0*Material(ii,1) - ... + Material(ii,3)*dt)./(2*epsilon_0*Material(ii,1) + ... + Material(ii,3)*dt); + K_b(ii) = 2*epsilon_0*Material(ii,1)./ ... + (2*epsilon_0*Material(ii,1) + Material(ii,3)*dt); + + K_a1(ii) = 1./(2*epsilon_0*Material(ii,1) + ... + Material(ii,3)*dt); + K_b1(ii) = 2*epsilon_0*Material(ii,1) - ... + Material(ii,3)*dt; +end + + +for i = mmm:nx-2 + for j = mmm:ny-2 + K_a_new(i,j) = K_a_new(i,j) + 1; + K_b_new(i,j) = K_b_new(i,j) + 1; + end +end + +% Replace [1,2] -> [1,-1] +for i = mmm:nx-2 + for j = mmm:ny-2 + if(K_a_new(i,j) == 2) + K_a_new(i,j) = K_a(2); + end + if(K_b_new(i,j) == 2) + K_b_new(i,j) = K_b(2); + end + end +end + + + +for i = mmm:nx-2 + for j = mmm:ny-2 + M0(i,j) = Material(Index(i+1,j+1)+1, 1); + M1(i,j) = Material(Index(i+1,j+1)+1, 2); + end +end + +for i = mmm:nx + for j = mmm:ny-1 + M2(i,j) = Material(Index(i,j)+1, 2); + end +end + +for i = mmm:nx-1 + for j = mmm:ny + M3(i,j) = Material(Index(i,j)+1, 2); + end +end + + +%% Create fullscreen figure and set a double-buffer +figure('units','normalized','outerposition',[0 0 1 1]); +set(gcf,'doublebuffer','on'); + +tic; + +% /////////////////// LOOP ///////////////////////////////////////////// +%% Main FDTD UPML routine with TF/SF excitation interface calculates TE and +%% TM modes simultaneously +for T = 1:number_of_iterations + %% Calculate incident 1D plain waves for TF/SF implementation + % TE mode (Hz) + for i = mmm:nx_b + Ey_1D(i) = k_Ey_a(i)*Ey_1D(i) ... + - k_Ey_b(i)*(Hz_1D(i+1)-Hz_1D(i)); + end + + t0 = 5 + tau = 15 + E0 = 2 + % Hz_1D(1) = E0*sin(2*pi*frequency*(T-1)*dt); + Hz_1D(1) = -E0 * ((T - t0) / tau) * exp(-1.0 * (((T - t0) / tau) ^ 2)); + + for i = mmm+1:nx_b + Hz_1D(i) = k_Hz_a(i)*Hz_1D(i) ... + - k_Hz_b(i)*(Ey_1D(i)-Ey_1D(i-1)); + end + + % TM mode (Ez) + for i = mmm:nx_b + Hy_1D(i) = k_Hy_a(i)*Hy_1D(i) + ... + k_Hy_b(i)*(Ez_1D(i+1)-Ez_1D(i)); + end + + % Ez_1D(1) = E0*sin(2*pi*frequency*(T-1)*dt); + Ez_1D(1) = -E0 * ((T - t0) / tau) * exp(-1.0 * (((T - t0) / tau) ^ 2)); + + for i = mmm:nx_b-1 + Fz_1D_r(i) = Fz_1D(i+1); + end + + for i = mmm+1:nx_b + Fz_1D(i) = k_Fz_a(i)*Fz_1D(i) + ... + k_Fz_b(i)*( Hy_1D(i) - Hy_1D(i-1) ); + end + + for i = mmm+1:nx_b + Ez_1D(i) = k_Ez_a*Ez_1D(i) + k_Ez_b*( Fz_1D(i) - Fz_1D_r(i-1)); + end + + %% Calculate Ez (TM mode) and Hz (TE mode) total fields + % TE: Wz -> Hz + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Fz_r1(i-1,j-1) = Wz(i,j); + end + end + + + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + + + Wz(i,j) = k_Fz_1_new(i-1,j-1) * Wz(i,j) + ... + k_Fz_2_new(i-1,j-1) * ((Ex(i,j)-Ex(i,j-1))/dy - ... + (Ey(i,j)-Ey(i-1,j))/dx ); + + % if(T == 10 && i == 10 && j == 10 ) + % k_Fz_1_new(i-1,j-1) * Wz(i,j) + % ((Ex(i,j)-Ex(i,j-1))/dy - (Ey(i,j)-Ey(i-1,j))/dx ) + % Wz(i,j) + % end + + + Hz(i,j) = k_Hz_1_new(i-1,j-1)*Hz(i,j) + ... + k_Hz_2_new(i-1,j-1) * ( Wz(i,j)-Fz_r1(i-1,j-1) )/M1(i-1,j-1); + end + end + + + + % TM: Fz -> Tz -> Ez + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Fz_r(i-1,j-1) = Fz(i,j); + Tz_r(i-1,j-1) = Tz(i,j); + end + end + + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Fz(i,j) = k_Fz_1_new(i-1,j-1) * Fz(i,j) + k_Fz_2_new(i-1,j-1)*( (Hy(i,j) - ... + Hy(i-1,j))/dx - (Hx(i,j) - ... + Hx(i,j-1))/dy ); + end + end + + + for i = mmm:nx-2 + for j = mmm:ny-2 + Tz(i+1,j+1) = K_a_new(i,j)* ... + Tz(i+1,j+1) + ... + K_b_new(i,j)*( Fz(i+1,j+1) - Fz_r(i,j)); + end + end + + + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Ez(i,j) = k_Ez_1_new(i-1,j-1)*Ez(i,j) + ... + k_Ez_2_new(i-1,j-1)*( Tz(i,j) - Tz_r(i-1,j-1) )/ M0(i-1,j-1); + end + end + + %% Calculate scattered field Ez and Hz in TF/SF + % TE mode + i = nx_a+1; + for j = ny_a+1:ny_b+1 + Hz(i,j) = Hz(i,j) + ... + dt / (mu_0*Material(Index(i,j)+1, 2)*dx) * Ey_1D(1); + + end + + i = nx_b+1; + for j = ny_a+1:ny_b+1 + Hz(i,j) = Hz(i,j) - ... + dt/(mu_0*Material(Index(i,j)+1, 2)*dx)*Ey_1D(nx_b-nx_a+2); + % j + end + + % TM mode + i = nx_a+1; + for j = ny_a+1:ny_b+1 + Ez(i,j) = Ez(i,j) - ... + dt./(epsilon_0*Material(Index(i,j)+1,1)*dx)*Hy_1D(1); + end + + i = nx_b+1; + for j = ny_a+1:ny_b+1 + Ez(i,j) = Ez(i,j) + ... + dt/(epsilon_0*Material(Index(i,j)+1,1)*dx)*Hy_1D(nx_b-nx_a+2); + end + + %% Calculate Hx and Ex total fields + % TE mode + for i = mmm:nx + for j = mmm:ny-1 + Gx_r1(i,j) = Mx(i,j); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Mx(i,j) = k_Gx_1_new(i,j) * Mx(i,j) + ... + k_Gx_2_new(i,j)*( Hz(i,j+1)-Hz(i,j) ); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Ex(i,j) = K_a1(IndexX(i,j)+1)*( K_b1(IndexX(i,j)+1)*... + Ex(i,j) + k_Ex_1_new(i,j)*Mx(i,j)-k_Ex_2_new(i,j)*Gx_r1(i,j) ); + % if(T == 10 && i == 2 && j == 2 ) + % k_Fz_1_new(i-1,j-1) * Wz(i,j) + % ((Ex(i,j)-Ex(i,j-1))/dy - (Ey(i,j)-Ey(i-1,j))/dx ) + % ( K_b1(IndexX(i,j)+1)*... + % Ex(i,j) + k_Ex_1_new(i,j)*Mx(i,j)-k_Ex_2_new(i,j)*Gx_r1(i,j) ) + % end + end + end + + % /// tmp comm start + % TM mode + for i = mmm:nx + for j = mmm:ny-1 + Gx_r(i,j) = Gx(i,j); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Gx(i,j) = k_Gx_1_new(i,j) * Gx(i,j) - ... + k_Gx_2_new(i,j) * ( Ez(i,j+1)-Ez(i,j) ); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Hx(i,j) = Hx(i,j) + (k_Hx_1_new(i,j) * Gx(i,j) - k_Hx_2_new(i,j)*Gx_r(i,j)) / M2(i,j); + end + end + + % %% Calculate Hy and Ey total fields + % TE mode + for i = mmm:nx-1 + for j = mmm:ny + Gy_r1(i,j) = My(i,j); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + My(i,j) = k_Gy_1_new(i,j)*My(i,j) - ... + k_Gy_2_new(i,j)*( Hz(i+1,j)-Hz(i,j)); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + Ey(i,j) = K_a1(IndexY(i,j)+1)*( K_b1(IndexY(i,j)+1)*... + Ey(i,j) + k_Ey_1_new(i,j)*My(i,j)-k_Ey_2_new(i,j)*Gy_r1(i,j) ); + end + end + + % TM mode + for i = mmm:nx-1 + for j = mmm:ny + Gy_r(i,j) = Gy(i,j); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + Gy(i,j) = k_Gy_1_new(i,j)*Gy(i,j) + ... + k_Gy_2_new(i,j)*(Ez(i+1,j)-Ez(i,j)); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + Hy(i,j) = Hy(i,j) + ... + (k_Hy_1_new(i,j)*Gy(i,j) - k_Hy_2_new(i,j)*Gy_r(i,j))/M3(i,j); + end + end + + % %% Calculate scattered field Hx and Ex in TF/SF + % % TE mode + j = ny_a; + for i = nx_a+1:nx_b+1 + Ex(i,j) = Ex(i,j) - ... + 2*dt/dy*K_a1(Index(i,j)+1)*Hz_1D(i-nx_a+1); + end + + j = ny_b+1; + for i = nx_a+1:nx_b+1 + Ex(i,j) = Ex(i,j) + ... + 2*dt/dy*K_a1(Index(i,j)+1)*Hz_1D(i-nx_a+1); + + end + + % TM mode + j = ny_a; + for i = nx_a+1:nx_b+1 + Hx(i,j) = Hx(i,j) + ... + dt/(mu_0*dy*Material(Index(i,j)+1,2))*Ez_1D(i-nx_a+1); + end + + j = ny_b+1; + for i = nx_a+1:nx_b+1 + Hx(i,j) = Hx(i,j) - ... + dt/(mu_0*dy*Material(Index(i,j)+1,2))*Ez_1D(i-nx_a+1); + end + + %% Calculate scattered field Hy and Ey in TF/SF + % TE mode + i = nx_a; + for j = ny_a+1:ny_b+1 + Ey(i,j) = Ey(i,j) + ... + 2*dt/dx*K_a1(Index(i,j)+1,1)*Hz_1D(mmm+1); + end + + i = nx_b+1; + for j = ny_a+1:ny_b+1 + Ey(i,j) = Ey(i,j) - ... + 2*dt/dx*K_a1(Index(i,j)+1,1)*Hz_1D(nx_b-nx_a+2); + + end + + % TM mode + i = nx_a; + for j = ny_a+1:ny_b+1 + Hy(i,j) = Hy(i,j) - ... + dt/(mu_0*dx*Material(Index(i,j)+1,2))*Ez_1D(2); + end + + i = nx_b+1; + for j = ny_a+1:ny_b+1 + Hy(i,j) = Hy(i,j) + ... + dt/(mu_0*dx*Material(Index(i,j)+1,2))*Ez_1D(nx_b-nx_a+2); + + if(T == 10) + j + dt/(mu_0*dx*Material(Index(i,j)+1,2))*Ez_1D(nx_b-nx_a+2) + end + end + % /// tmp comm end + + %% Plot Ez and Hz fields dynamics + if (mod(T,8) == 0) + subplot(1,2,1); + pcolor(Y(upml_width:ny-upml_width),X(upml_width:nx-upml_width),Ez(upml_width:nx-upml_width,upml_width:ny-upml_width)); + shading interp; + caxis([-E0 E0]); + axis image; + colorbar; + title('E_{z}(x,y)', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + xlabel('x, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + ylabel('y, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + subplot(1,2,2); + pcolor(Y(upml_width:ny-upml_width),X(upml_width:nx-upml_width),Hz(upml_width:nx-upml_width,upml_width:ny-upml_width)); + shading interp; + caxis([-E0 E0]); + axis image; + colorbar; + title('H_{z}(x,y)', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + xlabel('x, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + ylabel('y, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + drawnow; + end +end +disp(['Time elapsed - ',num2str(toc/60),' minutes']); diff --git a/src/fdtd/1d-pml/fdtd-pml-1d.h b/src/fdtd/1d-pml/fdtd-pml-1d.h index 8e5b581..c0c0174 100644 --- a/src/fdtd/1d-pml/fdtd-pml-1d.h +++ b/src/fdtd/1d-pml/fdtd-pml-1d.h @@ -7,6 +7,7 @@ #include // std::fill #include #include +#include class FdtdPml1D { // Speed of light in nm/fs. diff --git a/src/fdtd/2d-pml copy/fdtd-pml-2d.cpp b/src/fdtd/2d-pml copy/fdtd-pml-2d.cpp new file mode 100644 index 0000000..8baf84a --- /dev/null +++ b/src/fdtd/2d-pml copy/fdtd-pml-2d.cpp @@ -0,0 +1,187 @@ +#include "fdtd-pml-2d.h" + +void FdtdPml2D::SetParams(std::vector> &eps, + std::vector> &mu, + std::vector> &sigma, + size_t new_src_position_row, + size_t new_src_position_col) { + time_step = 0; + + src_row = new_src_position_row + pml_width; + src_col = new_src_position_col + pml_width; + + // Init field arrays. + std::fill_n(&dz[0][0], rows * cols, 0.0); + std::fill_n(&ez[0][0], rows * cols, 0.0); + std::fill_n(&hx[0][0], rows * cols, 0.0); + std::fill_n(&hy[0][0], rows * cols, 0.0); + std::fill_n(&gaz[0][0], rows * cols, 1.0); + + // Init pml arrays. + std::fill_n(&gi2[0], rows, 1.0); + std::fill_n(&gi3[0], rows, 1.0); + std::fill_n(&fi1[0], rows, 0.0); + std::fill_n(&fi2[0], rows, 1.0); + std::fill_n(&fi3[0], rows, 1.0); + + std::fill_n(&gj2[0], rows, 1.0); + std::fill_n(&gj3[0], rows, 1.0); + std::fill_n(&fj1[0], rows, 0.0); + std::fill_n(&fj2[0], rows, 1.0); + std::fill_n(&fj3[0], rows, 1.0); + + // Fill pml arrays. + for (size_t i = 0; i < pml_width; ++i) { + int xnum = pml_width - i; + double xd = pml_width; + double xxn = xnum / xd; + double xn = 0.33 * std::pow(xxn, 3); + + gi2[i] = 1.0 / (1.0 + xn); + gi2[rows - 1 - i] = 1.0 / (1.0 + xn); + gi3[i] = (1.0 - xn) / (1.0 + xn); + gi3[rows - i - 1] = (1.0 - xn) / (1.0 + xn); + + gj2[i] = 1.0 / (1.0 + xn); + gj2[rows - 1 - i] = 1.0 / (1.0 + xn); + gj3[i] = (1.0 - xn) / (1.0 + xn); + gj3[rows - i - 1] = (1.0 - xn) / (1.0 + xn); + + xxn = (xnum - 0.5) / xd; + xn = 0.33 * std::pow(xxn, 3); + + fi1[i] = xn; + fi1[rows - 2 - i] = xn; + fi2[i] = 1.0 / (1.0 + xn); + fi2[rows - 2 - i] = 1.0 / (1.0 + xn); + fi3[i] = (1.0 - xn) / (1.0 - xn); + fi3[rows - 2 - i] = (1.0 - xn) / (1.0 + xn); + + fj1[i] = xn; + fj1[rows - 2 - i] = xn; + fj2[i] = 1.0 / (1.0 + xn); + fj2[rows - 2 - i] = 1.0 / (1.0 + xn); + fj3[i] = (1.0 - xn) / (1.0 - xn); + fj3[rows - 2 - i] = (1.0 - xn) / (1.0 + xn); + } + + // Fill medium array. + for (size_t i = 0; i < rows; ++i) { + for (size_t j = 0; j < cols; ++j) { + if (i > pml_width && i < (rows - pml_width) && j > pml_width && j < (cols - pml_width)) { + gaz[i][j] = 1.0 / (eps[i - pml_width][j - pml_width] + (sigma[i - pml_width][j - pml_width] * dt) / epsz); + } else { + gaz[i][j] = 1.0 / (epsilon1 + (sigma1 * dt) / epsz); + } + + // Medium 1. + // gaz[i][j] = 1.0 / (epsilon1 + (sigma1 * dt) / epsz); + + // Medium 2. + // if (i>=80 && i<=90 && j>=80 && j<=120) { + // gaz[i][j] = 1.0 / (epsilon2 + (sigma2 * dt) / epsz); + // } + + // if (i>=120 && i<=130 && j>=80 && j<=120) { + // gaz[i][j] = 1.0 / (epsilon2 + (sigma2 * dt) / epsz); + // } + + // // Dielectric border. Medium 0. + // if (i >= (src_row - 20) && i <= (src_row + 20) && (j == (src_col - 10) || j == (src_col + 10)) || j >= (src_col - 10) && j <= (src_col + 10) && (i == (src_row - 10))) + // { + // gaz[i][j] = 1.0 / (epsilon0 + (sigma0 * dt) / epsz); + // } + } + } +} + +// One time layer calculation. +void FdtdPml2D::Calculation() { + // Dz field calculation. + for (int i = 1; i < rows; i++) { + for (int j = 1; j < cols; j++) { + dz[i][j] = gi3[i] * gj3[j] * dz[i][j] + gi2[i] * gj2[j] * 0.5 * (hy[i][j] - hy[i - 1][j] - hx[i][j] + hx[i][j - 1]); + } + } + + // Ez field calculation. + for (int i = 1; i < rows; i++) { + for (int j = 1; j < cols; j++) { + ez[i][j] = gaz[i][j] * dz[i][j]; + } + } + + // Put Gaussian beam source. + double source = -2.0 * ((time_step - t0) / tau) * std::exp(-1.0 * std::pow((time_step - t0) / tau, 2)); + ez[src_row][src_col] = source; + + for (int i = 0; i < rows - 1; i++) { + for (int j = 0; j < cols - 1; j++) { + // Hx field calculation. + hx[i][j] = fj3[j] * hx[i][j] + fj2[j] * 0.5 * (ez[i][j] - ez[i][j + 1]); + // Hy field calculation. + hy[i][j] = fi3[i] * hy[i][j] + fi2[i] * 0.5 * (ez[i + 1][j] - ez[i][j]); + } + } +} + +size_t FdtdPml2D::GetStep() { + return 1; +} + +size_t FdtdPml2D::GetCurrentTimeStep() { + return time_step; +} + +FdtdPml2D::FdtdPml2D() {} + +FdtdPml2D::FdtdPml2D(std::vector> &eps, + std::vector> &mu, + std::vector> &sigma, + size_t new_src_position_row, + size_t new_src_position_col) { + SetParams(eps, mu, sigma, new_src_position_row, new_src_position_col); +} + +void FdtdPml2D::CalcNextLayer(std::vector &vectX, + std::vector &vectY, + std::vector &vectEz, + std::vector &vectHy, + std::vector &vectHx, + std::vector &vectEnergy, + double &max, double &min) { + Calculation(); + + size_t step = GetStep(); + + // for (int xx = 1; xx < rows - 1; xx += step) + for (int xx = pml_width; xx < rows - pml_width; xx += step) { + // vectX.push_back(xx); + // vectY.push_back(xx); + // for (int yy = 1; yy < cols - 1; yy += step) + for (int yy = pml_width; yy < cols - pml_width; yy += step) { + // Energy + // double energy = yy1[xx][yy] * yy1[xx][yy] * Ez1[xx][yy] * Ez1[xx][yy] + + // Hy1[xx][yy] * Hy1[xx][yy] + Hx1[xx][yy] * Hx1[xx][yy]; + double energy = 1; + + vectX.push_back(xx - pml_width); + vectY.push_back(yy - pml_width); + vectEz.push_back(ez[yy][xx]); + vectHy.push_back(hy[yy][xx]); + vectHx.push_back(hx[yy][xx]); + vectEnergy.push_back(energy); + + if (xx > (pml_width + 10) && xx < (rows - pml_width - 10) && yy > (pml_width + 10) && yy < (rows - pml_width - 10)) { + if (ez[yy][xx] > max) { + max = ez[yy][xx]; + } + if (ez[yy][xx] < min) { + min = ez[yy][xx]; + } + } + } + } + + time_step++; +} \ No newline at end of file diff --git a/src/fdtd/2d-pml copy/fdtd-pml-2d.h b/src/fdtd/2d-pml copy/fdtd-pml-2d.h new file mode 100644 index 0000000..b7611df --- /dev/null +++ b/src/fdtd/2d-pml copy/fdtd-pml-2d.h @@ -0,0 +1,126 @@ +// 2D FDTD with PML. + +// Based on code from book: +// Electromagnetic simulation using the fdtd method with python. +// Chapter 3. + +#ifndef FDTD_PML_2D +#define FDTD_PML_2D + +#include +#include + +#include // std::fill_n +#include +#include + + +class FdtdPml2D { + static const size_t pml_width = 40; + + static const size_t rows = 220 + pml_width * 2; + static const size_t cols = 220 + pml_width * 2; + + // Set source position. + size_t src_row = rows / 2; + size_t src_col = cols / 2; + + // Light speed. + const double light_spd = 2.99792458e8; + + const double pi = 3.1415926535897932385; + + // Permeability of free space. + const double mu0 = 4.0 * pi * 1.0e-7; + + // Permittivity of free space. + const double epsz = 8.8e-12; + + // Permittivity. + double epsilon0 = 1; + double epsilon1 = 1; + double epsilon2 = 4.2; + + // Conductivity. + double sigma0 = 5000; + double sigma1 = 0.00001; + double sigma2 = 0.001; + + // Field arrays. + double dz[rows][cols]; + double ez[rows][cols]; + double hx[rows][cols]; + double hy[rows][cols]; + + double gaz[rows][cols]; + + // Space grid step. + const double ddx = 1.0e-3; + const double ddy = ddx; + + // Courant factor. + const double cfl_factor = 0.99; + + // Time step corressponding Courant factor. + const double dt = cfl_factor / (light_spd * sqrt(std::pow(1 / ddx, 2) + std::pow(1 / ddy, 2))); + + size_t time_step = 0; + + // Source params. + // Gaussian beam. + const size_t t0 = 20; + + // Beam width. + const size_t tau = 20; + + // PML params. + double gi2[rows]; + double gi3[rows]; + double fi1[rows]; + double fi2[rows]; + double fi3[rows]; + + double gj2[rows]; + double gj3[rows]; + double fj1[rows]; + double fj2[rows]; + double fj3[rows]; + + public: + static size_t GetRows() { + return rows - pml_width * 2; + } + + static size_t GetCols() { + return cols - pml_width * 2; + } + + size_t GetStep(); + size_t GetCurrentTimeStep(); + + void SetParams(std::vector> &eps, + std::vector> &mu, + std::vector> &sigma, + size_t new_src_position_row, + size_t new_src_position_col); + + void Calculation(); + + FdtdPml2D(); + + FdtdPml2D(std::vector> &eps, + std::vector> &mu, + std::vector> &sigma, + size_t new_src_position_row, + size_t new_src_position_col); + + void CalcNextLayer(std::vector &vectX, + std::vector &vectY, + std::vector &vectEz, + std::vector &vectHy, + std::vector &vectHx, + std::vector &vectEnergy, + double &max, double &min); +}; + +#endif \ No newline at end of file diff --git a/src/fdtd/2d-pml/fdtd-pml-2d.cpp b/src/fdtd/2d-pml/fdtd-pml-2d.cpp index 8baf84a..27f3073 100644 --- a/src/fdtd/2d-pml/fdtd-pml-2d.cpp +++ b/src/fdtd/2d-pml/fdtd-pml-2d.cpp @@ -1,187 +1,187 @@ -#include "fdtd-pml-2d.h" - -void FdtdPml2D::SetParams(std::vector> &eps, - std::vector> &mu, - std::vector> &sigma, - size_t new_src_position_row, - size_t new_src_position_col) { - time_step = 0; - - src_row = new_src_position_row + pml_width; - src_col = new_src_position_col + pml_width; - - // Init field arrays. - std::fill_n(&dz[0][0], rows * cols, 0.0); - std::fill_n(&ez[0][0], rows * cols, 0.0); - std::fill_n(&hx[0][0], rows * cols, 0.0); - std::fill_n(&hy[0][0], rows * cols, 0.0); - std::fill_n(&gaz[0][0], rows * cols, 1.0); - - // Init pml arrays. - std::fill_n(&gi2[0], rows, 1.0); - std::fill_n(&gi3[0], rows, 1.0); - std::fill_n(&fi1[0], rows, 0.0); - std::fill_n(&fi2[0], rows, 1.0); - std::fill_n(&fi3[0], rows, 1.0); - - std::fill_n(&gj2[0], rows, 1.0); - std::fill_n(&gj3[0], rows, 1.0); - std::fill_n(&fj1[0], rows, 0.0); - std::fill_n(&fj2[0], rows, 1.0); - std::fill_n(&fj3[0], rows, 1.0); - - // Fill pml arrays. - for (size_t i = 0; i < pml_width; ++i) { - int xnum = pml_width - i; - double xd = pml_width; - double xxn = xnum / xd; - double xn = 0.33 * std::pow(xxn, 3); - - gi2[i] = 1.0 / (1.0 + xn); - gi2[rows - 1 - i] = 1.0 / (1.0 + xn); - gi3[i] = (1.0 - xn) / (1.0 + xn); - gi3[rows - i - 1] = (1.0 - xn) / (1.0 + xn); - - gj2[i] = 1.0 / (1.0 + xn); - gj2[rows - 1 - i] = 1.0 / (1.0 + xn); - gj3[i] = (1.0 - xn) / (1.0 + xn); - gj3[rows - i - 1] = (1.0 - xn) / (1.0 + xn); - - xxn = (xnum - 0.5) / xd; - xn = 0.33 * std::pow(xxn, 3); - - fi1[i] = xn; - fi1[rows - 2 - i] = xn; - fi2[i] = 1.0 / (1.0 + xn); - fi2[rows - 2 - i] = 1.0 / (1.0 + xn); - fi3[i] = (1.0 - xn) / (1.0 - xn); - fi3[rows - 2 - i] = (1.0 - xn) / (1.0 + xn); - - fj1[i] = xn; - fj1[rows - 2 - i] = xn; - fj2[i] = 1.0 / (1.0 + xn); - fj2[rows - 2 - i] = 1.0 / (1.0 + xn); - fj3[i] = (1.0 - xn) / (1.0 - xn); - fj3[rows - 2 - i] = (1.0 - xn) / (1.0 + xn); - } - - // Fill medium array. - for (size_t i = 0; i < rows; ++i) { - for (size_t j = 0; j < cols; ++j) { - if (i > pml_width && i < (rows - pml_width) && j > pml_width && j < (cols - pml_width)) { - gaz[i][j] = 1.0 / (eps[i - pml_width][j - pml_width] + (sigma[i - pml_width][j - pml_width] * dt) / epsz); - } else { - gaz[i][j] = 1.0 / (epsilon1 + (sigma1 * dt) / epsz); - } - - // Medium 1. - // gaz[i][j] = 1.0 / (epsilon1 + (sigma1 * dt) / epsz); - - // Medium 2. - // if (i>=80 && i<=90 && j>=80 && j<=120) { - // gaz[i][j] = 1.0 / (epsilon2 + (sigma2 * dt) / epsz); - // } - - // if (i>=120 && i<=130 && j>=80 && j<=120) { - // gaz[i][j] = 1.0 / (epsilon2 + (sigma2 * dt) / epsz); - // } - - // // Dielectric border. Medium 0. - // if (i >= (src_row - 20) && i <= (src_row + 20) && (j == (src_col - 10) || j == (src_col + 10)) || j >= (src_col - 10) && j <= (src_col + 10) && (i == (src_row - 10))) - // { - // gaz[i][j] = 1.0 / (epsilon0 + (sigma0 * dt) / epsz); - // } - } - } -} - -// One time layer calculation. -void FdtdPml2D::Calculation() { - // Dz field calculation. - for (int i = 1; i < rows; i++) { - for (int j = 1; j < cols; j++) { - dz[i][j] = gi3[i] * gj3[j] * dz[i][j] + gi2[i] * gj2[j] * 0.5 * (hy[i][j] - hy[i - 1][j] - hx[i][j] + hx[i][j - 1]); - } - } - - // Ez field calculation. - for (int i = 1; i < rows; i++) { - for (int j = 1; j < cols; j++) { - ez[i][j] = gaz[i][j] * dz[i][j]; - } - } - - // Put Gaussian beam source. - double source = -2.0 * ((time_step - t0) / tau) * std::exp(-1.0 * std::pow((time_step - t0) / tau, 2)); - ez[src_row][src_col] = source; - - for (int i = 0; i < rows - 1; i++) { - for (int j = 0; j < cols - 1; j++) { - // Hx field calculation. - hx[i][j] = fj3[j] * hx[i][j] + fj2[j] * 0.5 * (ez[i][j] - ez[i][j + 1]); - // Hy field calculation. - hy[i][j] = fi3[i] * hy[i][j] + fi2[i] * 0.5 * (ez[i + 1][j] - ez[i][j]); - } - } -} - -size_t FdtdPml2D::GetStep() { - return 1; -} - -size_t FdtdPml2D::GetCurrentTimeStep() { - return time_step; -} - -FdtdPml2D::FdtdPml2D() {} - -FdtdPml2D::FdtdPml2D(std::vector> &eps, - std::vector> &mu, - std::vector> &sigma, - size_t new_src_position_row, - size_t new_src_position_col) { - SetParams(eps, mu, sigma, new_src_position_row, new_src_position_col); -} - -void FdtdPml2D::CalcNextLayer(std::vector &vectX, - std::vector &vectY, - std::vector &vectEz, - std::vector &vectHy, - std::vector &vectHx, - std::vector &vectEnergy, - double &max, double &min) { - Calculation(); - - size_t step = GetStep(); - - // for (int xx = 1; xx < rows - 1; xx += step) - for (int xx = pml_width; xx < rows - pml_width; xx += step) { - // vectX.push_back(xx); - // vectY.push_back(xx); - // for (int yy = 1; yy < cols - 1; yy += step) - for (int yy = pml_width; yy < cols - pml_width; yy += step) { - // Energy - // double energy = yy1[xx][yy] * yy1[xx][yy] * Ez1[xx][yy] * Ez1[xx][yy] + - // Hy1[xx][yy] * Hy1[xx][yy] + Hx1[xx][yy] * Hx1[xx][yy]; - double energy = 1; - - vectX.push_back(xx - pml_width); - vectY.push_back(yy - pml_width); - vectEz.push_back(ez[yy][xx]); - vectHy.push_back(hy[yy][xx]); - vectHx.push_back(hx[yy][xx]); - vectEnergy.push_back(energy); - - if (xx > (pml_width + 10) && xx < (rows - pml_width - 10) && yy > (pml_width + 10) && yy < (rows - pml_width - 10)) { - if (ez[yy][xx] > max) { - max = ez[yy][xx]; - } - if (ez[yy][xx] < min) { - min = ez[yy][xx]; - } - } - } - } - - time_step++; -} \ No newline at end of file +// #include "fdtd-pml-2d.h" + +// void FdtdPml2D::SetParams(std::vector> &eps, +// std::vector> &mu, +// std::vector> &sigma, +// size_t new_src_position_row, +// size_t new_src_position_col) { +// time_step = 0; + +// src_row = new_src_position_row + pml_width; +// src_col = new_src_position_col + pml_width; + +// // Init field arrays. +// std::fill_n(&dz[0][0], rows * cols, 0.0); +// std::fill_n(&ez[0][0], rows * cols, 0.0); +// std::fill_n(&hx[0][0], rows * cols, 0.0); +// std::fill_n(&hy[0][0], rows * cols, 0.0); +// std::fill_n(&gaz[0][0], rows * cols, 1.0); + +// // Init pml arrays. +// std::fill_n(&gi2[0], rows, 1.0); +// std::fill_n(&gi3[0], rows, 1.0); +// std::fill_n(&fi1[0], rows, 0.0); +// std::fill_n(&fi2[0], rows, 1.0); +// std::fill_n(&fi3[0], rows, 1.0); + +// std::fill_n(&gj2[0], rows, 1.0); +// std::fill_n(&gj3[0], rows, 1.0); +// std::fill_n(&fj1[0], rows, 0.0); +// std::fill_n(&fj2[0], rows, 1.0); +// std::fill_n(&fj3[0], rows, 1.0); + +// // Fill pml arrays. +// for (size_t i = 0; i < pml_width; ++i) { +// int xnum = pml_width - i; +// double xd = pml_width; +// double xxn = xnum / xd; +// double xn = 0.33 * std::pow(xxn, 3); + +// gi2[i] = 1.0 / (1.0 + xn); +// gi2[rows - 1 - i] = 1.0 / (1.0 + xn); +// gi3[i] = (1.0 - xn) / (1.0 + xn); +// gi3[rows - i - 1] = (1.0 - xn) / (1.0 + xn); + +// gj2[i] = 1.0 / (1.0 + xn); +// gj2[rows - 1 - i] = 1.0 / (1.0 + xn); +// gj3[i] = (1.0 - xn) / (1.0 + xn); +// gj3[rows - i - 1] = (1.0 - xn) / (1.0 + xn); + +// xxn = (xnum - 0.5) / xd; +// xn = 0.33 * std::pow(xxn, 3); + +// fi1[i] = xn; +// fi1[rows - 2 - i] = xn; +// fi2[i] = 1.0 / (1.0 + xn); +// fi2[rows - 2 - i] = 1.0 / (1.0 + xn); +// fi3[i] = (1.0 - xn) / (1.0 - xn); +// fi3[rows - 2 - i] = (1.0 - xn) / (1.0 + xn); + +// fj1[i] = xn; +// fj1[rows - 2 - i] = xn; +// fj2[i] = 1.0 / (1.0 + xn); +// fj2[rows - 2 - i] = 1.0 / (1.0 + xn); +// fj3[i] = (1.0 - xn) / (1.0 - xn); +// fj3[rows - 2 - i] = (1.0 - xn) / (1.0 + xn); +// } + +// // Fill medium array. +// for (size_t i = 0; i < rows; ++i) { +// for (size_t j = 0; j < cols; ++j) { +// if (i > pml_width && i < (rows - pml_width) && j > pml_width && j < (cols - pml_width)) { +// gaz[i][j] = 1.0 / (eps[i - pml_width][j - pml_width] + (sigma[i - pml_width][j - pml_width] * dt) / epsz); +// } else { +// gaz[i][j] = 1.0 / (epsilon1 + (sigma1 * dt) / epsz); +// } + +// // Medium 1. +// // gaz[i][j] = 1.0 / (epsilon1 + (sigma1 * dt) / epsz); + +// // Medium 2. +// // if (i>=80 && i<=90 && j>=80 && j<=120) { +// // gaz[i][j] = 1.0 / (epsilon2 + (sigma2 * dt) / epsz); +// // } + +// // if (i>=120 && i<=130 && j>=80 && j<=120) { +// // gaz[i][j] = 1.0 / (epsilon2 + (sigma2 * dt) / epsz); +// // } + +// // // Dielectric border. Medium 0. +// // if (i >= (src_row - 20) && i <= (src_row + 20) && (j == (src_col - 10) || j == (src_col + 10)) || j >= (src_col - 10) && j <= (src_col + 10) && (i == (src_row - 10))) +// // { +// // gaz[i][j] = 1.0 / (epsilon0 + (sigma0 * dt) / epsz); +// // } +// } +// } +// } + +// // One time layer calculation. +// void FdtdPml2D::Calculation() { +// // Dz field calculation. +// for (int i = 1; i < rows; i++) { +// for (int j = 1; j < cols; j++) { +// dz[i][j] = gi3[i] * gj3[j] * dz[i][j] + gi2[i] * gj2[j] * 0.5 * (hy[i][j] - hy[i - 1][j] - hx[i][j] + hx[i][j - 1]); +// } +// } + +// // Ez field calculation. +// for (int i = 1; i < rows; i++) { +// for (int j = 1; j < cols; j++) { +// ez[i][j] = gaz[i][j] * dz[i][j]; +// } +// } + +// // Put Gaussian beam source. +// double source = -2.0 * ((time_step - t0) / tau) * std::exp(-1.0 * std::pow((time_step - t0) / tau, 2)); +// ez[src_row][src_col] = source; + +// for (int i = 0; i < rows - 1; i++) { +// for (int j = 0; j < cols - 1; j++) { +// // Hx field calculation. +// hx[i][j] = fj3[j] * hx[i][j] + fj2[j] * 0.5 * (ez[i][j] - ez[i][j + 1]); +// // Hy field calculation. +// hy[i][j] = fi3[i] * hy[i][j] + fi2[i] * 0.5 * (ez[i + 1][j] - ez[i][j]); +// } +// } +// } + +// size_t FdtdPml2D::GetStep() { +// return 1; +// } + +// size_t FdtdPml2D::GetCurrentTimeStep() { +// return time_step; +// } + +// FdtdPml2D::FdtdPml2D() {} + +// FdtdPml2D::FdtdPml2D(std::vector> &eps, +// std::vector> &mu, +// std::vector> &sigma, +// size_t new_src_position_row, +// size_t new_src_position_col) { +// SetParams(eps, mu, sigma, new_src_position_row, new_src_position_col); +// } + +// void FdtdPml2D::CalcNextLayer(std::vector &vectX, +// std::vector &vectY, +// std::vector &vectEz, +// std::vector &vectHy, +// std::vector &vectHx, +// std::vector &vectEnergy, +// double &max, double &min) { +// Calculation(); + +// size_t step = GetStep(); + +// // for (int xx = 1; xx < rows - 1; xx += step) +// for (int xx = pml_width; xx < rows - pml_width; xx += step) { +// // vectX.push_back(xx); +// // vectY.push_back(xx); +// // for (int yy = 1; yy < cols - 1; yy += step) +// for (int yy = pml_width; yy < cols - pml_width; yy += step) { +// // Energy +// // double energy = yy1[xx][yy] * yy1[xx][yy] * Ez1[xx][yy] * Ez1[xx][yy] + +// // Hy1[xx][yy] * Hy1[xx][yy] + Hx1[xx][yy] * Hx1[xx][yy]; +// double energy = 1; + +// vectX.push_back(xx - pml_width); +// vectY.push_back(yy - pml_width); +// vectEz.push_back(ez[yy][xx]); +// vectHy.push_back(hy[yy][xx]); +// vectHx.push_back(hx[yy][xx]); +// vectEnergy.push_back(energy); + +// if (xx > (pml_width + 10) && xx < (rows - pml_width - 10) && yy > (pml_width + 10) && yy < (rows - pml_width - 10)) { +// if (ez[yy][xx] > max) { +// max = ez[yy][xx]; +// } +// if (ez[yy][xx] < min) { +// min = ez[yy][xx]; +// } +// } +// } +// } + +// time_step++; +// } \ No newline at end of file diff --git a/src/fdtd/2d-pml/fdtd-pml-2d.h b/src/fdtd/2d-pml/fdtd-pml-2d.h index 4122316..b6a062a 100644 --- a/src/fdtd/2d-pml/fdtd-pml-2d.h +++ b/src/fdtd/2d-pml/fdtd-pml-2d.h @@ -12,115 +12,926 @@ #include // std::fill_n #include +#include -class FdtdPml2D { - // Grid sizes. - static const size_t pml_width = 40; - - static const size_t rows = 220 + pml_width * 2; - static const size_t cols = 220 + pml_width * 2; - - // Set source position. - size_t src_row = rows / 2; - size_t src_col = cols / 2; - - // Light speed. - const double light_spd = 2.99792458e8; - - const double pi = 3.1415926535897932385; - - // Permeability of free space. - const double mu0 = 4.0 * pi * 1.0e-7; - - // Permittivity of free space. - const double epsz = 8.8e-12; - - // Permittivity. - double epsilon0 = 1; - double epsilon1 = 1; - double epsilon2 = 4.2; - - // Conductivity. - double sigma0 = 5000; - double sigma1 = 0.00001; - double sigma2 = 0.001; - // Field arrays. - double dz[rows][cols]; - double ez[rows][cols]; - double hx[rows][cols]; - double hy[rows][cols]; - - double gaz[rows][cols]; - - // Space grid step. - const double ddx = 1.0e-3; - const double ddy = ddx; +#include +#include +#include - // Courant factor. - const double cfl_factor = 0.99; +#define MEDIACONSTANT (2) +#define NUMBEROFITERATIONCONSTANT (1400) // was 300 +#define BIGLINESIZE (8192) - // Time step corressponding Courant factor. - const double dt = cfl_factor / (light_spd * sqrt(std::pow(1 / ddx, 2) + std::pow(1 / ddy, 2))); - size_t time_step = 0; +class FdtdPml2D { - // Source params. - // Gaussian beam. - const size_t t0 = 20; - // Beam width. - const size_t tau = 20; +/////////////////////////////////////////////////////// + // char ch; + double eps[MEDIACONSTANT] = {1.0, 1.0}; // index=0 is for vacuum, index=1 is for the metallic cylinder + double sig[MEDIACONSTANT] = {0.0, 1.0e+7}; + double mur[MEDIACONSTANT] = {1.0, 1.0}; + double sim[MEDIACONSTANT] = {0.0, 0.0}; + double ca[MEDIACONSTANT]; + double cb[MEDIACONSTANT]; + double da[MEDIACONSTANT]; + double db[MEDIACONSTANT]; + double source[NUMBEROFITERATIONCONSTANT]; + int n = 0; + int i,j, icenter, jcenter; + double eaf,haf, diam,rad,dist2; + double temporaryi,temporaryj,temporary; + double cc,muz,epsz,pi,freq,lambda,omega,dx,dt; + double rmax, orderbc; + static const size_t ie = 220; //number of grid cells in x-direction + static const size_t je = 220; //number of grid cells in y-direction + int ib, jb, is, js, nmax, iebc, jebc, ibbc, jbbc, iefbc, jefbc, ibfbc, jbfbc ; + int media; + double rtau, tau, delay; + double **ex; //fields in main grid + double **ey; + double **hz; + double **exbcf; //fields in front PML region + double **eybcf; + double **hzxbcf; + double **hzybcf; + double **exbcb; //fields in back PML region + double **eybcb; + double **hzxbcb; + double **hzybcb; + double **exbcl; //fields in left PML region + double **eybcl; + double **hzxbcl; + double **hzybcl; + double **exbcr; //fields in right PML region + double **eybcr; + double **hzxbcr; + double **hzybcr; + double **caex; // main grid coefficents + double **cbex; + double **caey; + double **cbey; + double **dahz; + double **dbhz; + double **caexbcf ; // pml coefficients + double **cbexbcf ; + double **caexbcb ; + double **cbexbcb ; + double **caexbcl ; + double **cbexbcl ; + double **caexbcr ; + double **cbexbcr ; + double **caeybcf ; + double **cbeybcf ; + double **caeybcb ; + double **cbeybcb ; + double **caeybcl ; + double **cbeybcl ; + double **caeybcr ; + double **cbeybcr ; + double **dahzxbcf; + double **dbhzxbcf; + double **dahzxbcb; + double **dbhzxbcb; + double **dahzxbcl; + double **dbhzxbcl; + double **dahzxbcr; + double **dbhzxbcr; + double **dahzybcf; + double **dbhzybcf; + double **dahzybcb; + double **dbhzybcb; + double **dahzybcl; + double **dbhzybcl; + double **dahzybcr; + double **dbhzybcr; + double delbc,sigmam,bcfactor; + double ca1,cb1,y1,y2, sigmay,sigmays,da1,db1; + double x1,x2,sigmax,sigmaxs; + double minimumValue, maximumValue; + char filename[BIGLINESIZE] ; + FILE *filePointer ; + double scaleValue; + int iValue,plottingInterval,centery,centerx ; - // PML params. - double gi2[rows]; - double gi3[rows]; - double fi1[rows]; - double fi2[rows]; - double fi3[rows]; - double gj2[rows]; - double gj3[rows]; - double fj1[rows]; - double fj2[rows]; - double fj3[rows]; public: - // Getters. - static size_t GetRows() { - return rows - pml_width * 2; - } - static size_t GetCols() { - return cols - pml_width * 2; + FdtdPml2D() {} + + size_t GetCurrentTimeStep() { + return n; + } + + struct Output { + // std::array, nx> Ez; + // std::array, nx> Hz; + // std::array Ez; + std::array Hz; + std::array X; + std::array Y; + size_t rows; + size_t cols; + // double maxEz; + // double minEz; + double maxHz; + double minHz; + }; + + struct Output GetValues() { + struct Output output; + + // output.Ez = Ez; + // output.Hz = Hz; + output.rows = ie; + output.cols = je; + + // size_t flatten_array_size = ie * nyje + for(size_t i = 0; i < ie; i += 1) { + for(size_t j = 0; j < je; j += 1) { + // output.Ez[i*je + j] = Ez[i][j]; + output.Hz[i*je + j] = hz[i][j]; + output.X[i*je + j] = i; + output.Y[i*je + j] = j; + } } - size_t GetStep(); - size_t GetCurrentTimeStep(); + // output.maxEz = *std::max_element(std::begin(output.Ez), std::end(output.Ez)); + // output.minEz = *std::min_element(std::begin(output.Ez), std::end(output.Ez)); + + output.maxHz = *std::max_element(std::begin(output.Hz), std::end(output.Hz)); + output.minHz = *std::min_element(std::begin(output.Hz), std::end(output.Hz)); + + return output; + } - void SetParams(std::vector> &eps, - std::vector> &mu, - std::vector> &sigma, - size_t new_src_position_row, - size_t new_src_position_col); - void Calculation(); - FdtdPml2D(); + void CalcNextLayer() { + + + //*********************************************************************** + // Update electric fields (EX and EY) in main grid + //*********************************************************************** + + // Note the 4 edges, ey left (i=0), ey right (i=ie), ex bottom (j=0) and ex top (j=je) are evaluated in the pml section + + for (i = 0; i < ie; i++) { + for (j = 1; j < je; j++) { // dont do ex at j=0 or j=je, it will be done in the PML section + ex[i][j] = caex[i][j] * ex[i][j] + cbex[i][j] * ( hz[i][j] - hz[i][j-1] ); + } /* jForLoop */ + } /* iForLoop */ + + for (i = 1; i < ie; i++) { // dont do ey at i=0 or i=ie, it will be done in the PML section + for (j = 0; j < je; j++) { + ey[i][j] = caey[i][j] * ey[i][j] + cbey[i][j] * ( hz[i-1][j] - hz[i][j] ); + } /* jForLoop */ + } /* iForLoop */ + + + //*********************************************************************** + // Update EX in PML regions + //*********************************************************************** + + + // FRONT + + for (i = 0; i < iefbc; i++) { + for (j = 1; j < jebc; j++) { // don't Evaluate exbcf at j=0, as it is the PEC + // note: sign change in the second term from main grid!! ... (due to the Exponential time stepping algorithm?) + exbcf[i][j] = caexbcf[i][j] * exbcf[i][j] - cbexbcf[i][j] * (hzxbcf[i][j-1] + hzybcf[i][j-1] - hzxbcf[i][j] - hzybcf[i][j]); + } /* jForLoop */ + } /* iForLoop */ + + for (j=0, i = 0; i < ie; i++) { // fill in the edge for ex at j=0 main grid (ties in the pml with the main grid) + ex[i][j] = caex[i][j] * ex[i][j] - cbex[i][j] * (hzxbcf[iebc+i][jebc-1] + hzybcf[iebc+i][jebc-1] - hz[i][j]); + } /* iForLoop */ + + + // BACK + + for (i = 0; i < iefbc; i++) { + for (j = 1; j < jebc; j++) { // don't Evaluate exbcb at j=jebc, as it is the PEC, also dont eval at j=0 as this point is the same as j=je on the main grid + exbcb[i][j] = caexbcb[i][j] * exbcb[i][j] - cbexbcb[i][j] * (hzxbcb[i][j-1] + hzybcb[i][j-1] - hzxbcb[i][j] - hzybcb[i][j]); + } /* jForLoop */ + } /* iForLoop */ + + for (j=je, i = 0; i < ie; i++) { // fill in the edge for ex at j=je main grid (ties in the pml with the main grid) + ex[i][j] = caex[i][j] * ex[i][j] - cbex[i][j] * (hz[i][j-1] - hzxbcb[iebc+i][0] - hzybcb[iebc+i][0]); + } /* iForLoop */ + + + // LEFT + + for (i = 0; i < iebc; i++) { + // don't Evaluate exbcl at j=0, j=0 is a special case, it needs data from the "front grid" see below + // likewise, don't Evaluate exbcl at j=je, j=je is a special case, it needs data from the "back grid" see below + for (j = 1; j < je; j++) { + exbcl[i][j] = caexbcl[i][j] * exbcl[i][j] - cbexbcl[i][j] * (hzxbcl[i][j-1] + hzybcl[i][j-1] - hzxbcl[i][j] - hzybcl[i][j]); + } /* jForLoop */ + } /* iForLoop */ + + for (j=0, i = 0; i < iebc; i++) { // exbcl at j=0 case, uses data from the "front grid" + exbcl[i][j] = caexbcl[i][j] * exbcl[i][j] - cbexbcl[i][j] * (hzxbcf[i][jebc-1] + hzybcf[i][jebc-1] - hzxbcl[i][j] - hzybcl[i][j]); + } /* iForLoop */ + + for (j=je, i = 0; i < iebc; i++) { // exbcl at j=je case, uses data from the "back grid" + exbcl[i][j] = caexbcl[i][j] * exbcl[i][j] - cbexbcl[i][j] * (hzxbcl[i][j-1] + hzybcl[i][j-1] - hzxbcb[i][0] - hzybcb[i][0]); + } /* iForLoop */ + + + // RIGHT + + for (i = 0; i < iebc; i++) { + // don't Evaluate exbcr at j=0, j=0 is a special case, it needs data from the "front grid" see below + // likewise, don't Evaluate exbcr at j=je, j=je is a special case, it needs data from the "back grid" see below + for (j = 1; j < je; j++) { + exbcr[i][j] = caexbcr[i][j] * exbcr[i][j] - cbexbcr[i][j] * (hzxbcr[i][j-1] + hzybcr[i][j-1] - hzxbcr[i][j] - hzybcr[i][j]); + } /* jForLoop */ + } /* iForLoop */ + + for (j=0, i = 0; i < iebc; i++) { // exbcr at j=0 case, uses data from the "front grid" (on the right side) + exbcr[i][j] = caexbcr[i][j] * exbcr[i][j] - cbexbcr[i][j] * (hzxbcf[iebc+ie + i][jebc-1] + hzybcf[iebc+ie + i][jebc-1] - hzxbcr[i][j] - hzybcr[i][j]); + } /* iForLoop */ + + for (j=je, i = 0; i < iebc; i++) { // exbcr at j=je case, uses data from the "back grid" + exbcr[i][j] = caexbcr[i][j] * exbcr[i][j] - cbexbcr[i][j] * (hzxbcr[i][j-1] + hzybcr[i][j-1] - hzxbcb[iebc+ie + i][0] - hzybcb[iebc+ie + i][0]); + } /* iForLoop */ + + + //*********************************************************************** + // Update EY in PML regions + //*********************************************************************** + + // FRONT + + for (i = 1; i < iefbc; i++) { // don't Evaluate eybcf at i=0 or iefbc, as it is the PEC + for (j = 0; j < jebc; j++) { + // note: sign change in the second term from main grid!! + eybcf[i][j] = caeybcf[i][j] * eybcf[i][j] - cbeybcf[i][j] * (hzxbcf[i][j] + hzybcf[i][j] - hzxbcf[i-1][j] - hzybcf[i-1][j]); + } /* jForLoop */ + } /* iForLoop */ + + + // BACK + + for (i = 1; i < iefbc; i++) { // don't Evaluate eybcb at i=0 or iefbc, as it is the PEC + for (j = 0; j < jebc; j++) { + eybcb[i][j] = caeybcb[i][j] * eybcb[i][j] - cbeybcb[i][j] * (hzxbcb[i][j] + hzybcb[i][j] - hzxbcb[i-1][j] - hzybcb[i-1][j]); + } /* jForLoop */ + } /* iForLoop */ + + + // LEFT + + for (i = 1; i < iebc; i++) { // don't Evaluate eybcb at i=0, as it is the PEC + for (j = 0; j < je; j++) { + eybcl[i][j] = caeybcl[i][j] * eybcl[i][j] - cbeybcl[i][j] * (hzxbcl[i][j] + hzybcl[i][j] - hzxbcl[i-1][j] - hzybcl[i-1][j] ); + } /* jForLoop */ + } /* iForLoop */ + + for (i=0, j = 0; j < je; j++) { // fill in the edge for ey at i=0 main grid (ties in the pml with the main grid) + ey[i][j] = caey[i][j] * ey[i][j] - cbey[i][j] * (hz[i][j] - hzxbcl[iebc-1][j] - hzybcl[iebc-1][j]); + } /* jForLoop */ + + + // RIGHT + + for (i = 1; i < iebc; i++) { // don't Evaluate eybcr at i=iebc, as it is the PEC, also dont eval at i=0 as this point is the same as i=ie on the main grid + for (j = 0; j < je; j++) { + eybcr[i][j] = caeybcr[i][j] * eybcr[i][j] - cbeybcr[i][j] * (hzxbcr[i][j] + hzybcr[i][j] - hzxbcr[i-1][j] - hzybcr[i-1][j] ); + } /* jForLoop */ + } /* iForLoop */ + + for (i=ie, j = 0; j < je; j++) { // fill in the edge for ey at i=ie main grid (ties in the pml with the main grid) + ey[i][j] = caey[i][j] * ey[i][j] - cbey[i][j] * (hzxbcr[0][j] + hzybcr[0][j] - hz[i-1][j] ); + } /* jForLoop */ + + + + //*********************************************************************** + // Update magnetic fields (HZ) in main grid + //*********************************************************************** + + + for (i = 0; i < ie; i++) { + for (j = 0; j < je; j++) { + hz[i][j] = dahz[i][j] * hz[i][j] + dbhz[i][j] * ( ex[i][j+1] - ex[i][j] + ey[i][j] - ey[i+1][j] ); + } /* jForLoop */ + } /* iForLoop */ + + + // rtau = 160.0e-12; + // tau = rtau / dt; + // delay = 3 * tau; + // for (i = 0; i < nmax; i++) { + // source[i] = 0.0; + // } /* iForLoop */ + + // int nn; + // for (nn = 0; nn < (int )(7.0 * tau); nn++) { + // temporary = (double )nn - delay; + // source[nn] = sin( omega * (temporary) * dt) * exp(-( (temporary * temporary)/(tau * tau) ) ); + // } /* forLoop */ + + const size_t t0 = 20; + + // Beam width. + const size_t tau = 20; + double src = -2.0 * ((n - t0) / tau) * std::exp(-1.0 * std::pow((n - t0) / tau, 2)); + + hz[is][js] = src; + // hz[is][js] = source[n]; + + + //*********************************************************************** + // Update HZX in PML regions + //*********************************************************************** + + // FRONT + + for (i = 0; i < iefbc; i++) { + for (j = 0; j < jebc; j++) { + // note: sign change in the second term from main grid!! + hzxbcf[i][j] = dahzxbcf[i][j] * hzxbcf[i][j] - dbhzxbcf[i][j] * (eybcf[i+1][j] - eybcf[i][j] ); + } /* jForLoop */ + } /* iForLoop */ + + // BACK + + for (i = 0; i < iefbc; i++) { + for (j = 0; j < jebc; j++) { + hzxbcb[i][j] = dahzxbcb[i][j] * hzxbcb[i][j] - dbhzxbcb[i][j] * (eybcb[i+1][j] - eybcb[i][j] ); + } /* jForLoop */ + } /* iForLoop */ + + // LEFT + + for (i = 0; i < (iebc-1); i++) { // don't evaluate hzxbcl at i=iebc-1 as it needs ey from main grid, see below + for (j = 0; j < je; j++) { + hzxbcl[i][j] = dahzxbcl[i][j] * hzxbcl[i][j] - dbhzxbcl[i][j] * (eybcl[i+1][j] - eybcl[i][j]); + } /* jForLoop */ + } /* iForLoop */ + + for (i=(iebc-1), j = 0; j < je; j++) { // fix-up hzxbcl at i=iebc-1 + hzxbcl[i][j] = dahzxbcl[i][j] * hzxbcl[i][j] - dbhzxbcl[i][j] * (ey[0][j] - eybcl[i][j] ); + } /* jForLoop */ + + // RIGHT + + for (i = 1; i < iebc; i++) { // don't evaluate hzxbcl at i=0 as it needs ey from main grid, see below + for (j = 0; j < je; j++) { + hzxbcr[i][j] = dahzxbcr[i][j] * hzxbcr[i][j] - dbhzxbcr[i][j] * (eybcr[i+1][j] - eybcr[i][j] ); + } /* jForLoop */ + } /* iForLoop */ + + for (i=0, j = 0; j < je; j++) { // fix-up hzxbcl at i=0 + hzxbcr[i][j] = dahzxbcr[i][j] * hzxbcr[i][j] - dbhzxbcr[i][j] * (eybcr[i+1][j] - ey[ie][j] ); + } /* jForLoop */ + + + + //*********************************************************************** + // Update HZY in PML regions + //*********************************************************************** + + // FRONT + + for (i = 0; i < iefbc; i++) { + for (j = 0; j < (jebc-1); j++) { // don't evaluate hzxbcf at j=jebc-1 as it needs data from main,left,right grids, see below + // note: sign change in the second term from main grid!! + hzybcf[i][j] = dahzybcf[i][j] * hzybcf[i][j] - dbhzybcf[i][j] * (exbcf[i][j] - exbcf[i][j+1] ); + } /* jForLoop */ + } /* iForLoop */ + + for (j = (jebc-1), i = 0; i < iebc; i++) { // fix-up hzybcf at j=jebc-1, with left grid + hzybcf[i][j] = dahzybcf[i][j] * hzybcf[i][j] - dbhzybcf[i][j] * (exbcf[i][j] - exbcl[i][0] ); + } /* iForLoop */ + + for (j = (jebc-1), i = 0; i < ie; i++) { // fix-up hzybcf at j=jebc-1, with main grid + hzybcf[iebc+i][j] = dahzybcf[iebc+i][j] * hzybcf[iebc+i][j] - dbhzybcf[iebc+i][j] * (exbcf[iebc+i][j] - ex[i][0]); + } /* iForLoop */ + + for (j = (jebc-1), i = 0; i < iebc; i++) { // fix-up hzybcf at j=jebc-1, with right grid + hzybcf[iebc+ie + i][j] = dahzybcf[iebc+ie + i][j] * hzybcf[iebc+ie + i][j] - dbhzybcf[iebc+ie + i][j] * (exbcf[iebc+ie + i][j] - exbcr[i][0]); + } /* iForLoop */ + + // BACK + + for (i = 0; i < iefbc; i++) { + for (j = 1; j < jebc; j++) { // don't evaluate hzxbcb at j=0 as it needs data from main,left,right grids, see below + hzybcb[i][j] = dahzybcb[i][j] * hzybcb[i][j] - dbhzybcb[i][j] * (exbcb[i][j] - exbcb[i][j+1] ); + } /* jForLoop */ + } /* iForLoop */ + + for (j = 0, i = 0; i < iebc; i++) { // fix-up hzybcb at j=0, with left grid + hzybcb[i][j] = dahzybcb[i][j] * hzybcb[i][j] - dbhzybcb[i][j] * (exbcl[i][je] - exbcb[i][j+1] ); + } /* iForLoop */ + + for (j = 0, i = 0; i < ie; i++) { // fix-up hzybcb at j=0, with main grid + hzybcb[iebc+i][j] = dahzybcb[iebc+i][j] * hzybcb[iebc+i][j] - dbhzybcb[iebc+i][j] * (ex[i][je] - exbcb[iebc+i][j+1]); + } /* iForLoop */ + + for (j = 0, i = 0; i < iebc; i++) { // fix-up hzybcb at j=0, with right grid + hzybcb[iebc+ie + i][j] = dahzybcb[iebc+ie + i][j] * hzybcb[iebc+ie + i][j] - dbhzybcb[iebc+ie + i][j] * (exbcr[i][je] - exbcb[iebc+ie + i][j+1] ); + } /* iForLoop */ + + // LEFT + + for (i = 0; i < iebc; i++) { + for (j = 0; j < je; j++) { + hzybcl[i][j] = dahzybcl[i][j] * hzybcl[i][j] - dbhzybcl[i][j] * (exbcl[i][j] - exbcl[i][j+1] ); + } /* jForLoop */ + } /* iForLoop */ + + // RIGHT + + for (i = 0; i < iebc; i++) { + for (j = 0; j < je; j++) { + hzybcr[i][j] = dahzybcr[i][j] * hzybcr[i][j] - dbhzybcr[i][j] * (exbcr[i][j] - exbcr[i][j+1] ); + } /* jForLoop */ + } /* iForLoop */ + + n++; + } - FdtdPml2D(std::vector> &eps, - std::vector> &mu, - std::vector> &sigma, - size_t new_src_position_row, - size_t new_src_position_col); +// standard C memory allocation for 2-D array +double **AllocateMemory (int imax, int jmax, double initialValue) +{ + int i,j; + double **pointer; + pointer = (double **)malloc(imax * sizeof(double *)); + for (i = 0; i < imax; i++) { + pointer[i] = (double *)malloc(jmax * sizeof(double)); + for (j = 0; j < jmax; j++) { + pointer[i][j] = initialValue; + } /* jForLoop */ + } /* iForLoop */ + return(pointer) ; +} + + + void InitializeFdtd () + { + + + + //*********************************************************************** + // Printing/Plotting variables + //*********************************************************************** + minimumValue = -0.1; + maximumValue = 0.1; + plottingInterval = 0; + centery = 25 ; + centerx = 15 ; + + //*********************************************************************** + // Fundamental constants + //*********************************************************************** + pi = (acos(-1.0)); + cc = 2.99792458e8; //speed of light in free space (meters/second) + muz = 4.0 * pi * 1.0e-7; //permeability of free space + epsz = 1.0 / (cc * cc * muz); //permittivity of free space + + freq = 5.0e+9; //center frequency of source excitation (Hz) + lambda = cc / freq; //center wavelength of source excitation + omega = 2.0 * pi * freq; //center frequency in radians + + + //*********************************************************************** + // Grid parameters + //*********************************************************************** + + // ie = 100; //number of grid cells in x-direction + // je = 50; //number of grid cells in y-direction + + ib = ie + 1; // one extra is needed for fields on the boundaries (ie Ex on top boundary, Ey on right boundary) + jb = je + 1; // ditto + + is = 15; //location of z-directed hard source + js = je / 2; //location of z-directed hard source + + dx = 3.0e-3; //space increment of square lattice (meters) + dt = dx / (2.0 * cc); //time step, seconds, courant limit, Taflove1995 page 177 + + nmax = NUMBEROFITERATIONCONSTANT; //total number of time steps + + iebc = 8; //thickness of left and right PML region + jebc = 8; //thickness of front and back PML region + rmax = 0.00001; // R(0) reflection coefficient (in %) Nikolova part4 p.25 + orderbc = 2; // m, grading order, optimal values: 2 <= m <= 6, Nikolova part4 p.29 + ibbc = iebc + 1; + jbbc = jebc + 1; + iefbc = ie + 2 * iebc; // for front and bottom (width of region) + jefbc = je + 2 * jebc; // not used + ibfbc = iefbc + 1; // one extra for Ey on right boundary + jbfbc = jefbc + 1; // not used + + + //*********************************************************************** + // Material parameters + //*********************************************************************** + + media = MEDIACONSTANT; // number of different medias, ie 2: vacuum, metallicCylinder + + + //*********************************************************************** + // Wave excitation + //*********************************************************************** + + rtau = 160.0e-12; + tau = rtau / dt; + delay = 3 * tau; + for (i = 0; i < nmax; i++) { + source[i] = 0.0; + } /* iForLoop */ + + int nn; + for (nn = 0; nn < (int )(7.0 * tau); nn++) { + temporary = (double )nn - delay; + source[nn] = sin( omega * (temporary) * dt) * exp(-( (temporary * temporary)/(tau * tau) ) ); + } /* forLoop */ + + + //*********************************************************************** + // Field arrays + //*********************************************************************** + + ex = AllocateMemory(ie,jb, 0.0 ); //fields in main grid + ey = AllocateMemory(ib,je, 0.0 ); + hz = AllocateMemory(ie,je, 0.0 ); + + exbcf = AllocateMemory(iefbc,jebc, 0.0 ); //fields in front PML region + eybcf = AllocateMemory(ibfbc,jebc, 0.0 ); + hzxbcf = AllocateMemory(iefbc,jebc, 0.0 ); + hzybcf = AllocateMemory(iefbc,jebc, 0.0 ); + + exbcb = AllocateMemory(iefbc,jbbc, 0.0 ); //fields in back PML region + eybcb = AllocateMemory(ibfbc,jebc, 0.0 ); + hzxbcb = AllocateMemory(iefbc,jebc, 0.0 ); + hzybcb = AllocateMemory(iefbc,jebc, 0.0 ); + + exbcl = AllocateMemory(iebc,jb, 0.0 ); //fields in left PML region + eybcl = AllocateMemory(iebc,je, 0.0 ); + hzxbcl = AllocateMemory(iebc,je, 0.0 ); + hzybcl = AllocateMemory(iebc,je, 0.0 ); + + exbcr = AllocateMemory(iebc,jb, 0.0 ); //fields in right PML region + eybcr = AllocateMemory(ibbc,je, 0.0 ); + hzxbcr = AllocateMemory(iebc,je, 0.0 ); + hzybcr = AllocateMemory(iebc,je, 0.0 ); + + + //*********************************************************************** + // Updating coefficients + //*********************************************************************** + + for (i = 0; i < media; i++) { + eaf = dt * sig[i] / (2.0 * epsz * eps[i] ); // Taflove1995 p.67 + ca[i] = (1.0 - eaf) / (1.0 + eaf); // ditto + cb[i] = dt / epsz / eps[i] / dx / (1.0 + eaf); // ditto + haf = dt * sim[i] / (2.0 * muz * mur[i]); // ditto + da[i] = (1.0 - haf) / (1.0 + haf); // ditto + db[i] = dt / muz / mur[i] / dx / (1.0 + haf); // ditto + } /* iForLoop */ + + + //*********************************************************************** + // Geometry specification (main grid) + //*********************************************************************** + + // Initialize entire main grid to free space + + caex = AllocateMemory(ie,jb, ca[0]); + cbex = AllocateMemory(ie,jb, cb[0]); + + caey = AllocateMemory(ib,je, ca[0] ); + cbey = AllocateMemory(ib,je, cb[0] ); + + dahz = AllocateMemory(ie,je, da[0] ); + dbhz = AllocateMemory(ie,je, db[0] ); + + // Add metal cylinder + + diam = 20; // diameter of cylinder: 6 cm + rad = diam / 2.0; // radius of cylinder: 3 cm + icenter = (4 * ie) / 5; // i-coordinate of cylinder's center + jcenter = je / 2; // j-coordinate of cylinder's center + for (i = 0; i < ie; i++) { + for (j = 0; j < je; j++) { + temporaryi = (double )(i - icenter) ; + temporaryj = (double )(j - jcenter) ; + dist2 = (temporaryi + 0.5) * (temporaryi + 0.5) + (temporaryj) * (temporaryj); + if (dist2 <= (rad * rad)) { + caex[i][j] = ca[1]; + cbex[i][j] = cb[1]; + } /* if */ + // This looks tricky! Why can't caey/cbey use the same 'if' statement as caex/cbex above ?? + dist2 = (temporaryj + 0.5) * (temporaryj + 0.5) + (temporaryi) * (temporaryi); + if (dist2 <= (rad * rad)) { + caey[i][j] = ca[1]; + cbey[i][j] = cb[1]; + } /* if */ + } /* jForLoop */ + } /* iForLoop */ + + + //*********************************************************************** + // Fill the PML regions + //*********************************************************************** + + caexbcf = AllocateMemory(iefbc,jebc, 0.0 ); + cbexbcf = AllocateMemory(iefbc,jebc, 0.0 ); + caexbcb = AllocateMemory(iefbc,jbbc, 0.0 ); + cbexbcb = AllocateMemory(iefbc,jbbc, 0.0 ); + caexbcl = AllocateMemory(iebc,jb, 0.0 ); + cbexbcl = AllocateMemory(iebc,jb, 0.0 ); + caexbcr = AllocateMemory(iebc,jb, 0.0 ); + cbexbcr = AllocateMemory(iebc,jb, 0.0 ); + caeybcf = AllocateMemory(ibfbc,jebc, 0.0 ); + cbeybcf = AllocateMemory(ibfbc,jebc, 0.0 ); + caeybcb = AllocateMemory(ibfbc,jebc, 0.0 ); + cbeybcb = AllocateMemory(ibfbc,jebc, 0.0 ); + caeybcl = AllocateMemory(iebc,je, 0.0 ); + cbeybcl = AllocateMemory(iebc,je, 0.0 ); + caeybcr = AllocateMemory(ibbc,je, 0.0 ); + cbeybcr = AllocateMemory(ibbc,je, 0.0 ); + dahzxbcf = AllocateMemory(iefbc,jebc, 0.0 ); + dbhzxbcf = AllocateMemory(iefbc,jebc, 0.0 ); + dahzxbcb = AllocateMemory(iefbc,jebc, 0.0 ); + dbhzxbcb = AllocateMemory(iefbc,jebc, 0.0 ); + dahzxbcl = AllocateMemory(iebc,je, 0.0 ); + dbhzxbcl = AllocateMemory(iebc,je, 0.0 ); + dahzxbcr = AllocateMemory(iebc,je, 0.0 ); + dbhzxbcr = AllocateMemory(iebc,je, 0.0 ); + dahzybcf = AllocateMemory(iefbc,jebc, 0.0 ); + dbhzybcf = AllocateMemory(iefbc,jebc, 0.0 ); + dahzybcb = AllocateMemory(iefbc,jebc, 0.0 ); + dbhzybcb = AllocateMemory(iefbc,jebc, 0.0 ); + dahzybcl = AllocateMemory(iebc,je, 0.0 ); + dbhzybcl = AllocateMemory(iebc,je, 0.0 ); + dahzybcr = AllocateMemory(iebc,je, 0.0 ); + dbhzybcr = AllocateMemory(iebc,je, 0.0 ); + + delbc = (double )iebc * dx; // width of PML region (in mm) + + // SigmaMaximum, using polynomial grading (Nikolova part 4, p.30), rmax=reflectionMax in percent + sigmam =-log(rmax/100.0) * epsz * cc * (orderbc + 1) / (2 * delbc); + + // bcfactor comes from the polynomial grading equation: sigma_x = sigmaxMaximum * (x/d)^m, where d=width of PML, m=gradingOrder, (Nikolova part4, p.28) + // IMPORTANT: The conductivity (sigma) must use the "average" value at each mesh point as follows: + // sigma_x = sigma_Maximum/dx * Integral_from_x0_to_x1 of (x/d)^m dx, where x0=currentx-0.5, x1=currentx+0.5 (Nikolova part 4, p.32) + // integrating gives: sigma_x = (sigmaMaximum / (dx * d^m * m+1)) * ( x1^(m+1) - x0^(m+1) ) (Nikolova part 4, p.32) + // the first part is "bcfactor", so, sigma_x = bcfactor * ( x1^(m+1) - x0^(m+1) ) (Nikolova part 4, p.32) + // note: it's not exactly clear what the term eps[0] is for. It's probably to cover the case in which eps[0] is not equal to one (ie the main grid area next to the pml boundary is not vacuum) + bcfactor = eps[0] * sigmam / ( dx * (pow(delbc,orderbc)) * (orderbc + 1)); + + // FRONT region + + for (i = 0; i < iefbc; i++) { // IS THIS EVER USED? -- coef for the pec ex at j=0 set so that ex_t+1 = ex_t, but ex is never evaluated at j=0. + caexbcf[i][0] = 1.0; + cbexbcf[i][0] = 0.0; + } /* iForLoop */ + + for (j = 1; j < jebc; j++) { // calculate the coefs for the PML layer (except for the boundary at the main grid, which is a special case (see below)) + y1 = ((double )(jebc - j) + 0.5) * dx; // upper bounds for point j (re-adujsted for C indexes!) + y2 = ((double )(jebc - j) - 0.5) * dx; // lower bounds for point j + sigmay = bcfactor * (pow(y1,(orderbc+1)) - pow(y2,(orderbc+1)) ); // polynomial grading + ca1 = exp (-sigmay * dt / (epsz * eps[0])); // exponential time step, Taflove1995 p.77,78 + cb1 = (1.0 - ca1) / (sigmay * dx); // ditto, but note sign change from Taflove1995 + for (i = 0; i < iefbc; i++) { + caexbcf[i][j] = ca1; + cbexbcf[i][j] = cb1; + } /* iForLoop */ + } /* jForLoop */ + + sigmay = bcfactor * pow( (0.5 * dx), (orderbc+1) ); // calculate for the front edge of the pml at j=0 in the main grid (half vacuum (sigma=0) and half pml) + ca1 = exp( -sigmay * dt / (epsz * eps[0]) ); + cb1 = ( 1 - ca1) / (sigmay * dx); + for (i = 0; i < ie; i++) { + caex[i][0] = ca1; + cbex[i][0] = cb1; + } /* iForLoop */ + + for (i = 0; i < iebc; i++) { // this continues the front edge into the left and right grids + caexbcl[i][0] = ca1; + cbexbcl[i][0] = cb1; + caexbcr[i][0] = ca1; + cbexbcr[i][0] = cb1; + } /* iForLoop */ + + for (j = 0; j < jebc; j++) { // for ey and hz (which are offset spacially 1/2 dx from ex) + y1 = ((double )(jebc - j) + 0.0) * dx; // upper bounds for point j + y2 = ((double )(jebc - j) - 1.0) * dx; // lower bounds for point j + sigmay = bcfactor * (pow(y1,(orderbc+1)) - pow(y2,(orderbc+1)) ); // polynomial grading + sigmays = sigmay * (muz / (epsz*eps[0]) ); // Taflove1995 p.182 (for no reflection: sigmaM = sigmaE * mu0/eps0) + da1 = exp (-sigmays * dt / muz); // exponential time step, Taflove1995 p.77,78 + db1 = (1.0 - da1) / (sigmays * dx); + for (i = 0; i < iefbc; i++) { + dahzybcf[i][j] = da1; + dbhzybcf[i][j] = db1; + dahzxbcf[i][j] = da[0]; // important note: hzx is Perpendicular to the front pml and so is not attenuated (sigma=0) (looks like vacuum) + dbhzxbcf[i][j] = db[0]; // ditto + } /* iForLoop */ + for (i = 0; i < ibfbc; i++) { + caeybcf[i][j] = ca[0]; // important note: ey is Perpendicular to the front pml and so is not attenuated (sigma=0) (looks like vacuum) + cbeybcf[i][j] = cb[0]; // ditto + } /* iForLoop */ + } /* jForLoop */ + + + // BACK region + + for (j=jebc, i = 0; i < iefbc; i++) { // IS THIS EVER USED? -- coef for the pec ex at j=jebc set to ex_t+1 = ex_t + caexbcb[i][j] = 1.0; + cbexbcb[i][j] = 0.0; + } /* iForLoop */ + + for (j = 1; j < jebc; j++) { // calculate the coefs for the PML layer (except for the boundary at the main grid, which is a special case (see below)) + y1 = ((double )(j) + 0.5) * dx; // upper bounds for point j (re-adujsted for C indexes!) + y2 = ((double )(j) - 0.5) * dx; // lower bounds for point j + sigmay = bcfactor * (pow(y1,(orderbc+1)) - pow(y2,(orderbc+1)) ); // polynomial grading + ca1 = exp (-sigmay * dt / (epsz * eps[0])); // exponential time step + cb1 = (1.0 - ca1) / (sigmay * dx); // ditto, but note sign change from Taflove + for (i = 0; i < iefbc; i++) { + caexbcb[i][j] = ca1; + cbexbcb[i][j] = cb1; + } /* iForLoop */ + } /* jForLoop */ + + sigmay = bcfactor * pow( (0.5 * dx), (orderbc+1) ); // calculate for the front edge of the pml at j=0 in the main grid (half vacuum (sigma=0) and half pml) + ca1 = exp( -sigmay * dt / (epsz * eps[0]) ); + cb1 = ( 1 - ca1) / (sigmay * dx); + for (j=je, i = 0; i < ie; i++) { + caex[i][j] = ca1; + cbex[i][j] = cb1; + } /* iForLoop */ + + for (j=je, i = 0; i < iebc; i++) { // this continues the back edge into the left and right grids + caexbcl[i][j] = ca1; + cbexbcl[i][j] = cb1; + caexbcr[i][j] = ca1; + cbexbcr[i][j] = cb1; + } /* iForLoop */ + + for (j = 0; j < jebc; j++) { // for ey and hz (which are offset spacially 1/2 dx from ex) + y1 = ((double )(j) + 1.0) * dx; // upper bounds for point j + y2 = ((double )(j) + 0.0) * dx; // lower bounds for point j + sigmay = bcfactor * (pow(y1,(orderbc+1)) - pow(y2,(orderbc+1)) ); // polynomial grading + sigmays = sigmay * (muz / (epsz*eps[0]) ); + da1 = exp (-sigmays * dt / muz); + db1 = (1.0 - da1) / (sigmays * dx); + for (i = 0; i < iefbc; i++) { + dahzybcb[i][j] = da1; + dbhzybcb[i][j] = db1; + dahzxbcb[i][j] = da[0]; // important note: hzx is Perpendicular to the back pml and so is not attenuated (sigma=0) (looks like vacuum) + dbhzxbcb[i][j] = db[0]; // ditto + } /* iForLoop */ + for (i = 0; i < ibfbc; i++) { + caeybcb[i][j] = ca[0]; // important note: ey is Perpendicular to the back pml and so is not attenuated (sigma=0) (looks like vacuum) + cbeybcb[i][j] = cb[0]; // ditto + } /* iForLoop */ + } /* jForLoop */ + + // LEFT region + + + for (i=0, j = 0; j < je; j++) { // IS THIS EVER USED? -- coef for the pec ey at i=0 set to ey_t+1 = ey_t + caeybcl[i][j] = 1.0; + cbeybcl[i][j] = 0.0; + } /* jForLoop */ + + for (i = 1; i < iebc; i++) { // calculate the coefs for the PML layer (except for the boundary at the main grid, which is a special case (see below)) + x1 = ((double )(iebc - i) + 0.5) * dx; // upper bounds for point i (re-adujsted for C indexes!) + x2 = ((double )(iebc - i) - 0.5) * dx; // lower bounds for point i + sigmax = bcfactor * ( pow(x1,(orderbc+1)) - pow(x2,(orderbc+1)) ); + ca1 = exp(-sigmax * dt / (epsz * eps[0]) ); + cb1 = (1.0 - ca1) / (sigmax * dx); + for (j = 0; j < je; j++) { + caeybcl[i][j] = ca1; + cbeybcl[i][j] = cb1; + } /* jForLoop */ + for (j = 0; j < jebc; j++) { // fill in the front left and back left corners for ey + caeybcf[i][j] = ca1; + cbeybcf[i][j] = cb1; + caeybcb[i][j] = ca1; + cbeybcb[i][j] = cb1; + } /* jForLoop */ + } /* iForLoop */ + + sigmax = bcfactor * pow( (0.5 * dx), (orderbc+1) ); // calculate for the left edge of the pml at x=0 in the main grid (half vacuum (sigma=0) and half pml) + ca1 = exp(-sigmax * dt / (epsz * eps[0]) ); + cb1 = (1.0 - ca1) / (sigmax * dx); + for (i=0, j = 0; j < je; j++) { + caey[i][j] = ca1; + cbey[i][j] = cb1; + } /* jForLoop */ + for (i=iebc, j = 0; j < jebc; j++) { // continue the left edge into the front and back grids + caeybcf[i][j] = ca1; + cbeybcf[i][j] = cb1; + caeybcb[i][j] = ca1; + cbeybcb[i][j] = cb1; + } /* jForLoop */ + + + for (i = 0; i < iebc; i++) { // for ex and hz (which are offset spacially 1/2 dx from ey) + x1 = ((double )(iebc - i) + 0.0) * dx; // upper bounds for point i (re-adujsted for C indexes!) + x2 = ((double )(iebc - i) - 1.0) * dx; // lower bounds for point i + sigmax = bcfactor * ( pow(x1,(orderbc+1)) - pow(x2,(orderbc+1)) ); + sigmaxs = sigmax * (muz / (epsz * eps[0]) ); + da1 = exp(-sigmaxs * dt / muz); + db1 = (1 - da1) / (sigmaxs * dx); + for (j = 0; j < je; j++) { + dahzxbcl[i][j] = da1; + dbhzxbcl[i][j] = db1; + dahzybcl[i][j] = da[0]; // important note: hzy is Perpendicular to the left pml and so is not attenuated (sigma=0) (looks like vacuum) + dbhzybcl[i][j] = db[0]; + } /* jForLoop */ + for (j = 0; j < jebc; j++) { // fill in the front left and back left corners for hzx + dahzxbcf[i][j] = da1; + dbhzxbcf[i][j] = db1; + dahzxbcb[i][j] = da1; + dbhzxbcb[i][j] = db1; + } /* jForLoop */ + for (j = 1; j < je; j++) { // important note: ex is Perpendicular to the left pml and so is not attenuated (sigma=0) (looks like vacuum) + caexbcl[i][j] =ca[0]; + cbexbcl[i][j] =cb[0]; + } /* jForLoop */ + } /* iForLoop */ + + + + // RIGHT region + + for (i=iebc, j = 0; j < je; j++) { // IS THIS EVER USED? -- coef for the pec ey at i=iebc set to ey_t+1 = ey_t + caeybcr[i][j] = 1.0; + cbeybcr[i][j] = 0.0; + } /* jForLoop */ + + for (i = 1; i < iebc; i++) { // calculate the coefs for the PML layer (except for the boundary at the main grid, which is a special case (see below)) + x1 = ((double )(i) + 0.5) * dx; // upper bounds for point i (re-adujsted for C indexes!) + x2 = ((double )(i) - 0.5) * dx; // lower bounds for point i + sigmax = bcfactor * ( pow(x1,(orderbc+1)) - pow(x2,(orderbc+1)) ); + ca1 = exp(-sigmax * dt / (epsz * eps[0]) ); + cb1 = (1.0 - ca1) / (sigmax * dx); + for (j = 0; j < je; j++) { + caeybcr[i][j] = ca1; + cbeybcr[i][j] = cb1; + } /* jForLoop */ + for (j = 0; j < jebc; j++) { // fill in the front right and back right corners for ey + caeybcf[i+iebc+ie][j] = ca1; + cbeybcf[i+iebc+ie][j] = cb1; + caeybcb[i+iebc+ie][j] = ca1; + cbeybcb[i+iebc+ie][j] = cb1; + } /* jForLoop */ + } /* iForLoop */ + + + sigmax = bcfactor * pow( (0.5 * dx), (orderbc+1) ); // calculate for the right edge of the pml at x=ic in the main grid (half vacuum (sigma=0) and half pml) + ca1 = exp(-sigmax * dt / (epsz * eps[0]) ); + cb1 = (1.0 - ca1) / (sigmax * dx); + for (i=ie, j = 0; j < je; j++) { + caey[i][j] = ca1; + cbey[i][j] = cb1; + } /* jForLoop */ + for (i=iebc+ie, j = 0; j < jebc; j++) { // continue the right edge into the front and back grids + caeybcf[i][j] = ca1; + cbeybcf[i][j] = cb1; + caeybcb[i][j] = ca1; + cbeybcb[i][j] = cb1; + } /* jForLoop */ + + + for (i = 0; i < iebc; i++) { // for ex and hz (which are offset spacially 1/2 dx from ey) + x1 = ((double )(i) + 1.0) * dx; // upper bounds for point i (re-adujsted for C indexes!) + x2 = ((double )(i) + 0.0) * dx; // lower bounds for point i + sigmax = bcfactor * ( pow(x1,(orderbc+1)) - pow(x2,(orderbc+1)) ); + sigmaxs = sigmax * (muz / (epsz * eps[0]) ); + da1 = exp(-sigmaxs * dt / muz); + db1 = (1 - da1) / (sigmaxs * dx); + for (j = 0; j < je; j++) { + dahzxbcr[i][j] = da1; + dbhzxbcr[i][j] = db1; + dahzybcr[i][j] = da[0]; // important note: hzy is Perpendicular to the right pml and so is not attenuated (sigma=0) (looks like vacuum) + dbhzybcr[i][j] = db[0]; + } /* jForLoop */ + for (j = 0; j < jebc; j++) { // fill in the front right and back right corners for hzx + dahzxbcf[i+ie+iebc][j] = da1; + dbhzxbcf[i+ie+iebc][j] = db1; + dahzxbcb[i+ie+iebc][j] = da1; + dbhzxbcb[i+ie+iebc][j] = db1; + } /* jForLoop */ + for (j = 1; j < je; j++) { // important note: ex is Perpendicular to the right pml and so is not attenuated (sigma=0) (looks like vacuum) + caexbcr[i][j] =ca[0]; + cbexbcr[i][j] =cb[0]; + } /* jForLoop */ + } /* iForLoop */ + - void CalcNextLayer(std::vector &vectX, - std::vector &vectY, - std::vector &vectEz, - std::vector &vectHy, - std::vector &vectHx, - std::vector &vectEnergy, - double &max, double &min); + } + }; #endif \ No newline at end of file diff --git a/src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.cpp b/src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.cpp new file mode 100644 index 0000000..21c6467 --- /dev/null +++ b/src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.cpp @@ -0,0 +1,909 @@ +#include "fdtd-2d-upml-tf-sf.h" + +TFSF::TFSF() { + SetParams(); +} + +size_t TFSF::GetCurrentTimeStep() { + return T; +} + +struct TFSF::Output TFSF::GetValues() { + struct Output output; + + // output.Ez = Ez; + // output.Hz = Hz; + output.rows = nx; + output.cols = ny; + + // size_t flatten_array_size = nx * ny; + + for(size_t i = 0; i < nx; i += 1) { + for(size_t j = 0; j < ny; j += 1) { + output.Ez[i*ny + j] = Ez[i][j]; + output.Hz[i*ny + j] = Hz[i][j]; + output.X[i*ny + j] = i; + output.Y[i*ny + j] = j; + } + } + + output.maxEz = *std::max_element(std::begin(output.Ez), std::end(output.Ez)); + output.minEz = *std::min_element(std::begin(output.Ez), std::end(output.Ez)); + + output.maxHz = *std::max_element(std::begin(output.Hz), std::end(output.Hz)); + output.minHz = *std::min_element(std::begin(output.Hz), std::end(output.Hz)); + + return output; +} + + +void TFSF::CalcNextLayer() { + // Calculate incident 1D plain waves for TF/SF implementation + // TE mode (Hz) + + // H + for(size_t i=0; i < nx_b; i += 1) { + Ey_1D[i] = k_Ey_a[i]*Ey_1D[i] - k_Ey_b[i]*(Hz_1D[i+1]-Hz_1D[i]); // C + } + + // std::cout << "Ey_1d\n"; + + Hz_1D[0] = E0*std::sin(2 * pi * frequency * T * dt); + + // H + for(size_t i=1; i < nx_b; i += 1) { + Hz_1D[i] = k_Hz_a[i] * Hz_1D[i] - k_Hz_b[i]*(Ey_1D[i]-Ey_1D[i-1]); // C + } + + + // TM mode (Ez) + for(size_t i=0; i < nx_b; i += 1) { + Hy_1D[i] = k_Hy_a[i]*Hy_1D[i] + k_Hy_b[i]*(Ez_1D[i+1]-Ez_1D[i]); // C + } + + + Ez_1D[0] = E0*std::sin(2*pi*frequency*T*dt); + + for(size_t i=0; i < nx_b-1; i += 1) { + Fz_1D_r[i] = Fz_1D[i+1]; // C + } + + for(size_t i=1; i < nx_b; i += 1) { + Fz_1D[i] = k_Fz_a[i]*Fz_1D[i] + k_Fz_b[i]*(Hy_1D[i] - Hy_1D[i-1]); // C + } + + for(size_t i=1; i < nx_b; i += 1) { + Ez_1D[i] = k_Ez_a*Ez_1D[i] + k_Ez_b*( Fz_1D[i] - Fz_1D_r[i-1]); // C + } + + + // Calculate Ez (TM mode) and Hz (TE mode) total fields + // TE: Wz -> Hz + + // h + // C + for(size_t i=1; i < nx-1; i += 1) { + for(size_t j=1; j < ny-1; j += 1) { + Fz_r1[i-1][j-1] = Wz[i][j]; + } + } + + // h + + for(size_t i=1; i < nx-1; i += 1) { + for(size_t j=1; j < ny-1; j += 1) { + + + // C + Wz[i][j] = k_Fz_1_new[i-1][j-1] * Wz[i][j] + + k_Fz_2_new[i-1][j-1] * ((Ex[i][j]-Ex[i][j-1])/dy - + (Ey[i][j]-Ey[i-1][j])/dx ); + + + // C + Hz[i][j] = k_Hz_1_new[i-1][j-1]*Hz[i][j] + + k_Hz_2_new[i-1][j-1] * ( Wz[i][j]-Fz_r1[i-1][j-1] )/M1[i-1][j-1]; + } + } + + + + // TM: Fz -> Tz -> Ez + // h + for(size_t i=1; i < nx-1; i += 1) { + for(size_t j=1; j < ny-1; j += 1) { + Fz_r[i-1][j-1] = Fz[i][j]; + Tz_r[i-1][j-1] = Tz[i][j]; + } + } + + + // h + for(size_t i=1; i < nx-1; i += 1) { + for(size_t j=1; j < ny-1; j += 1) { + Fz[i][j] = k_Fz_1_new[i-1][j-1] * Fz[i][j] + k_Fz_2_new[i-1][j-1]*( (Hy[i][j] - + Hy[i-1][j])/dx - (Hx[i][j] - + Hx[i][j-1])/dy); + } + } + + // h + for(size_t i=0; i < nx-2; i += 1) { + for(size_t j=0; j < ny-2; j += 1) { + Tz[i+1][j+1] = K_a_new[i][j] * Tz[i+1][j+1] + + K_b_new[i][j]*( Fz[i+1][j+1] - Fz_r[i][j]); + } + } + + // h + for(size_t i=1; i < nx-1; i += 1) { + for(size_t j=1; j < ny-1; j += 1) { + Ez[i][j] = k_Ez_1_new[i-1][j-1]*Ez[i][j] + + k_Ez_2_new[i-1][j-1]*( Tz[i][j] - Tz_r[i-1][j-1] )/ M0[i-1][j-1]; + } + } + + // Calculate scattered field Ez and Hz in TF/SF + // TE mode + size_t i = nx_a; + + //h + for(size_t j=ny_a; j <= ny_b; j += 1) { + Hz[i][j] = Hz[i][j] + dt / (mu_0*Material[Index[i][j]][1]*dx) * Ey_1D[0]; + } + + //h + i = nx_b; + for(size_t j=ny_a; j <= ny_b; j += 1) { + Hz[i][j] = Hz[i][j] - dt/(mu_0*Material[Index[i][j]][1]*dx) + *Ey_1D[nx_b-nx_a+1]; // C + } + + // TM mode + // h + i = nx_a; + for(size_t j=ny_a; j <= ny_b; j += 1) { + Ez[i][j] = Ez[i][j] - dt / (epsilon_0*Material[Index[i][j]][0]*dx)*Hy_1D[0]; + } + + // h + i = nx_b; + for(size_t j=ny_a; j <= ny_b; j += 1) { + Ez[i][j] = Ez[i][j] + dt/(epsilon_0*Material[Index[i][j]][0]*dx)*Hy_1D[nx_b-nx_a+1]; + } + + // Calculate Hx and Ex total fields + // TE mode + // h + for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny-1; j += 1) { + Gx_r1[i][j] = Mx[i][j]; + } + } + + // h + // C + for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny-1; j += 1) { + Mx[i][j] = k_Gx_1_new[i][j] * Mx[i][j] + + k_Gx_2_new[i][j]*( Hz[i][j+1]-Hz[i][j] ); + } + } + + // h + // + for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny-1; j += 1) { + Ex[i][j] = K_a1[IndexX[i][j]] * (K_b1[IndexX[i][j]] * + Ex[i][j] + k_Ex_1_new[i][j]*Mx[i][j]-k_Ex_2_new[i][j]*Gx_r1[i][j]); + + // if(T == 9 && i == 1 && j == 1 ) { + // std::cout << "a: " << k_Fz_1_new[i-1][j-1] * Wz[i][j] << "\n"; + // std::cout << "b: " << ((Ex[i][j]-Ex[i][j-1])/dy - + // (Ey[i][j]-Ey[i-1][j])/dx ) << "\n"; + // std::cout << "ex: " << (K_b1[IndexX[i][j]] * + // Ex[i][j] + k_Ex_1_new[i][j]*Mx[i][j]-k_Ex_2_new[i][j]*Gx_r1[i][j]) << "\n"; + // } + } + } + + // TM mode + // h + for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny-1; j += 1) { + Gx_r[i][j] = Gx[i][j]; + } + } + + // // h + for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny-1; j += 1) { + Gx[i][j] = k_Gx_1_new[i][j] * Gx[i][j] - + k_Gx_2_new[i][j] * (Ez[i][j+1] - Ez[i][j]); + } + } + + // // h + for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny-1; j += 1) { + Hx[i][j] = Hx[i][j] + (k_Hx_1_new[i][j] * Gx[i][j] - k_Hx_2_new[i][j]*Gx_r[i][j]) / M2[i][j]; + } + } + + + // // Calculate Hy and Ey total fields + // TE mode + // h + for(size_t i=0; i < nx-1; i += 1) { + for(size_t j=0; j < ny; j += 1) { + Gy_r1[i][j] = My[i][j]; + } + } + + // h + // C + for(size_t i=0; i < nx-1; i += 1) { + for(size_t j=0; j < ny; j += 1) { + My[i][j] = k_Gy_1_new[i][j]*My[i][j] - + k_Gy_2_new[i][j] * (Hz[i+1][j]-Hz[i][j]); + } + } + + + // h + // C + for(size_t i=0; i < nx-1; i += 1) { + for(size_t j=0; j < ny; j += 1) { + Ey[i][j] = K_a1[IndexY[i][j]] * (K_b1[IndexY[i][j]] * + Ey[i][j] + k_Ey_1_new[i][j] * My[i][j] - k_Ey_2_new[i][j] * Gy_r1[i][j]); + } + } + + + + + // ??? + // TM mode + // h + // C + for(size_t i=0; i < nx-1; i += 1) { + for(size_t j=0; j < ny; j += 1) { + Gy_r[i][j] = Gy[i][j]; + } + } + + // h + // C + for(size_t i=0; i < nx-1; i += 1) { + for(size_t j=0; j < ny; j += 1) { + Gy[i][j] = k_Gy_1_new[i][j] * Gy[i][j] + + k_Gy_2_new[i][j] * (Ez[i+1][j]-Ez[i][j]); + } + } + + // h + // C + for(size_t i=0; i < nx-1; i += 1) { + for(size_t j=0; j < ny; j += 1) { + Hy[i][j] = Hy[i][j] + + (k_Hy_1_new[i][j] * Gy[i][j] - k_Hy_2_new[i][j] * Gy_r[i][j]) / M3[i][j]; + } + } + + + // // Calculate scattered field Hx and Ex in TF/SF + // // TE mode + // h + // C + size_t j = ny_a-1; + for(size_t i=nx_a; i <= nx_b; i += 1) { + Ex[i][j] = Ex[i][j] - + 2 * dt / dy * K_a1[Index[i][j]] * Hz_1D[i-nx_a+1]; + + } + + // h + j = ny_b; + for(size_t i=nx_a; i <= nx_b; i += 1) { + Ex[i][j] = Ex[i][j] + + 2 * dt / dy * K_a1[Index[i][j]] * Hz_1D[i-nx_a+1]; + + + } + + // TM mode + //h + j = ny_a-1; + for(size_t i=nx_a; i <= nx_b; i += 1) { + Hx[i][j] = Hx[i][j] + + dt / (mu_0 * dy * Material[Index[i][j]][1]) * Ez_1D[i-nx_a+1]; + } + + // h + j = ny_b; + for(size_t i=nx_a; i <= nx_b; i += 1) { + Hx[i][j] = Hx[i][j] - + dt / (mu_0 * dy * Material[Index[i][j]][1]) * Ez_1D[i-nx_a+1]; + } + + // Calculate scattered field Hy and Ey in TF/SF + // TE mode + // h + i = nx_a-1; + for(size_t j=ny_a; j <= ny_b; j += 1) { + Ey[i][j] = Ey[i][j] + + 2 * dt / dx * K_a1[Index[i][j]] * Hz_1D[1]; + } // K_a1???? + + // h + i = nx_b; + for(size_t j=ny_a; j <= ny_b; j += 1) { + Ey[i][j] = Ey[i][j] - + 2 * dt / dx * K_a1[Index[i][j]] * Hz_1D[nx_b-nx_a+1]; + + + }// K_a1???? + + // TM mode + // h + i = nx_a-1; + for(size_t j=ny_a; j <= ny_b; j += 1) { + Hy[i][j] = Hy[i][j] - + dt / (mu_0 * dx * Material[Index[i][j]][1]) * Ez_1D[1]; + } + + + // h + // C + i = nx_b; + for(size_t j=ny_a; j <= ny_b; j += 1) { + Hy[i][j] = Hy[i][j] + + dt / (mu_0 * dx * Material[Index[i][j]][1]) * Ez_1D[nx_b-nx_a+1]; + + // if(T == 9 ) { + // std::cout << j << ": " << dt / (mu_0 * dx * Material[Index[i][j]][1]) * Ez_1D[nx_b-nx_a+1] << " "; + + // } + // std::cout << "\n"; + } + + + if(T == 9) { + // for(size_t i=0; i < Fz_r.size(); i += 1) { + // for(size_t j=0; j < Fz_r[0].size(); j += 1) { + // std::cout << Fz_r[i][j] << " "; + // } + // std::cout << "\n"; + // } + // << "Fz_r1[0][j]" << "\n"; + std::ofstream myfile; + myfile.open("Wz.txt"); + if (myfile.is_open()) + { + // myfile << "This is a line.\n"; + // myfile << "This is another line.\n"; + + for(size_t i=0; i < Hz.size(); i += 1) { + for(size_t j=0; j < Hz[0].size(); j += 1) { + + myfile << Hz[i][j] << " "; + } + myfile << '\n'; + } + +// std::cout << "filfe" << '\n'; + myfile.close(); + } + else std::cout << "Unable to open file"; + } + + + T++; +} + + +void TFSF::SetParams() { + + // Free space FDTD coefficients + k_Fz_1.fill(1.0); + k_Fz_2.fill(dt); + + k_Ez_1.fill(1.0); + k_Ez_2.fill(1.0/epsilon_0); + k_Gx_1.fill(1.0); + k_Gx_2.fill(dt/dy); + k_Hx_1.fill(1.0/mu_0); + k_Hx_2.fill(1.0/mu_0); + k_Gy_1.fill(1.0); + k_Gy_2.fill(dt/dx); + k_Hy_1.fill(1.0/mu_0); + k_Hy_2.fill(1.0/mu_0); + k_Hz_1.fill(1.0); + k_Hz_2.fill(1.0/mu_0); + k_Ex_1.fill(2.0); + k_Ex_2.fill(2.0); + k_Ey_1.fill(2.0); + k_Ey_2.fill(2.0); + + + + k_Ez_a = (2.0*epsilon_0 * Material[0][0] - Material[0][2]*dt)/(2.0*epsilon_0*Material[0][0] + Material[0][2]*dt); + k_Ez_b = 2.0/(2.0*epsilon_0*Material[0][0] + Material[0][2]*dt); + + + // // Free space plane wave coefficients + k_Hy_a.fill(1.0); + k_Hy_b.fill(dt/(mu_0*Material[0][1]*dx)); + k_Fz_a.fill(1.0); + k_Fz_b.fill(dt/dx); + + k_Hz_a.fill(1.0); + k_Hz_b.fill(dt/(mu_0*Material[0][1]*dx)); + k_Ey_a.fill(1.0); + k_Ey_b.fill(dt/(epsilon_0*Material[0][0]*dx)); + + + + eta = std::sqrt(mu_0*Material[0][1]/epsilon_0/Material[0][0]); + + // Field transformation coefficients in upml_width areas for TM ans TE modes + // Along x-axis + sigma_max = -(m+1)*std::log(R_err)/(2*eta*upml_width*dx); + + for(size_t i=0; i < upml_width; i += 1) { + + size_t end; + + double sigma_x = sigma_max * std::pow((double)(upml_width-i)/upml_width, m); + double ka_x = 1+(ka_max-1)*std::pow((double)(upml_width-i)/upml_width, m); + + k_Ez_1[i] = (2*epsilon_0*ka_x-sigma_x*dt)/(2*epsilon_0*ka_x+sigma_x*dt); + end = k_Ez_1.size() - 1; + k_Ez_1[end-i] = k_Ez_1[i]; + + k_Ez_2[i] = 2.0/(2*epsilon_0*ka_x + sigma_x*dt); + end = k_Ez_2.size() - 1; + k_Ez_2[end-i] = k_Ez_2[i]; + + k_Hx_1[i] = (2*epsilon_0*ka_x+sigma_x*dt)/(2*epsilon_0*mu_0); + end = k_Hx_1.size() - 1; + k_Hx_1[end-i] = k_Hx_1[i]; + + k_Hx_2[i] = (2*epsilon_0*ka_x-sigma_x*dt)/(2*epsilon_0*mu_0); + end = k_Hx_2.size() - 1; + k_Hx_2[end-i] = k_Hx_2[i]; + + k_Hz_1[i] = (2*epsilon_0*ka_x-sigma_x*dt)/(2*epsilon_0*ka_x+sigma_x*dt); + end = k_Hz_1.size() - 1; + k_Hz_1[end-i] = k_Hz_1[i]; + + k_Hz_2[i] = 2*epsilon_0/(2*epsilon_0*ka_x+sigma_x*dt)/mu_0; + end = k_Hz_2.size() - 1; + k_Hz_2[end-i] = k_Hz_2[i]; + + k_Ex_1[i] = (2*epsilon_0*ka_x+sigma_x*dt)/epsilon_0; + end = k_Ex_1.size() - 1; + k_Ex_1[end-i] = k_Ex_1[i]; + + k_Ex_2[i] = (2*epsilon_0*ka_x-sigma_x*dt)/epsilon_0; + end = k_Ex_2.size() - 1; + k_Ex_2[end-i] = k_Ex_2[i]; + + + sigma_x = sigma_max * std::pow((double)(upml_width-i-0.5)/upml_width, m); + ka_x = 1+(ka_max-1)*std::pow((double)(upml_width-i-0.5)/upml_width, m); + + k_Gy_1[i] = (2*epsilon_0*ka_x-sigma_x*dt)/(2*epsilon_0*ka_x+sigma_x*dt); + end = k_Gy_1.size() - 1; + k_Gy_1[end-i] = k_Gy_1[i]; + + k_Gy_2[i] = 2*epsilon_0*dt/(2*epsilon_0*ka_x+sigma_x*dt)/dx; + end = k_Gy_2.size() - 1; + k_Gy_2[end-i] = k_Gy_2[i]; + } + + + + +// % Along y-axis + + sigma_max = -(m+1)*std::log(R_err)/(2*eta*upml_width*dy); + + for(size_t i=0; i < upml_width; i += 1) { + + size_t end; + + double sigma_y = sigma_max * std::pow((double)(upml_width-i)/upml_width, m); + double ka_y = 1+(ka_max-1)*std::pow((double)(upml_width-i)/upml_width, m); + + k_Fz_1[i] = (2*epsilon_0*ka_y-sigma_y*dt)/(2*epsilon_0*ka_y+sigma_y*dt); + end = k_Fz_1.size() - 1; + k_Fz_1[end-i] = k_Fz_1[i]; + + k_Fz_2[i] = 2*epsilon_0*dt/(2*epsilon_0*ka_y+sigma_y*dt); + end = k_Fz_2.size() - 1; + k_Fz_2[end-i] = k_Fz_2[i]; + + k_Hy_1[i] = (2*epsilon_0*ka_y+sigma_y*dt)/(2*epsilon_0*mu_0); + end = k_Hy_1.size() - 1; + k_Hy_1[end-i] = k_Hy_1[i]; + + k_Hy_2[i] = (2*epsilon_0*ka_y-sigma_y*dt)/(2*epsilon_0*mu_0); + end = k_Hy_2.size() - 1; + k_Hy_2[end-i] = k_Hy_2[i]; + + k_Ey_1[i] = (2*epsilon_0*ka_y+sigma_y*dt)/epsilon_0; + end = k_Ey_1.size() - 1; + k_Ey_1[end-i] = k_Ey_1[i]; + + k_Ey_2[i] = (2*epsilon_0*ka_y-sigma_y*dt)/epsilon_0; + end = k_Ey_2.size() - 1; + k_Ey_2[end-i] = k_Ey_2[i]; + + + sigma_y = sigma_max * std::pow((double)(upml_width-i-0.5)/upml_width, m); + ka_y = 1+(ka_max-1)*std::pow((double)(upml_width-i-0.5)/upml_width, m); + + k_Gx_1[i] = (2*epsilon_0*ka_y-sigma_y*dt)/(2*epsilon_0*ka_y+sigma_y*dt); + end = k_Gx_1.size() - 1; + k_Gx_1[end-i] = k_Gx_1[i]; + + k_Gx_2[i] = 2*epsilon_0*dt/(2*epsilon_0*ka_y+sigma_y*dt)/dy; + end = k_Gx_2.size() - 1; + k_Gx_2[end-i] = k_Gx_2[i]; + } + + +// for (auto & row : arr) { +// for (auto & col : row) { +// std::cout << col << ' '; +// } +// std::cout << std::endl; +// } + + +// % Vectorize transformation coefficients +size_t rows; +size_t cols; + +//nx-2 ny-2 +rows = k_Fz_1_new.size(); +cols = k_Fz_1_new[0].size(); +std::fill_n(&k_Fz_1_new[0][0], rows * cols, k_Fz_1[upml_width]); // C +std::fill_n(&k_Fz_2_new[0][0], rows * cols, k_Fz_2[upml_width]); + +// left right border +for(size_t i=0; i < rows; i += 1) { + for(size_t j=0; j < upml_width-1; j += 1) { + // left + k_Fz_1_new[i][j] = k_Fz_1[j+1]; // C + k_Fz_2_new[i][j] = k_Fz_2[j+1]; + + // right + k_Fz_1_new[i][cols-j-1] = k_Fz_1[j+1]; // C + k_Fz_2_new[i][cols-j-1] = k_Fz_2[j+1]; + } +} + + + +rows = k_Ez_1_new.size(); +cols = k_Ez_1_new[0].size(); +std::fill_n(&k_Ez_1_new[0][0], rows * cols, k_Ez_1[upml_width]); +std::fill_n(&k_Ez_2_new[0][0], rows * cols, k_Ez_2[upml_width]); + +// up down border +for(size_t i=0; i < upml_width-1; i += 1) { + for(size_t j=0; j < cols; j += 1) { + + // up + k_Ez_1_new[i][j] = k_Ez_1[i+1]; + k_Ez_2_new[i][j] = k_Ez_2[i+1]; + + //down + k_Ez_1_new[rows-i-1][j] = k_Ez_1[i+1]; + k_Ez_2_new[rows-i-1][j] = k_Ez_2[i+1]; + } +} + + +rows = k_Gx_1_new.size(); +cols = k_Gx_1_new[0].size(); +std::fill_n(&k_Gx_1_new[0][0], rows * cols, k_Gx_1[upml_width]); +std::fill_n(&k_Gx_2_new[0][0], rows * cols, k_Gx_2[upml_width]); + +// left right border +for(size_t i=0; i < rows; i += 1) { + for(size_t j=0; j < upml_width; j += 1) { + // left + k_Gx_1_new[i][j] = k_Gx_1[j]; + k_Gx_2_new[i][j] = k_Gx_2[j]; + + // right + k_Gx_1_new[i][cols-j-1] = k_Gx_1[j]; + k_Gx_2_new[i][cols-j-1] = k_Gx_2[j]; + } +} + + +rows = k_Hx_1_new.size(); +cols = k_Hx_1_new[0].size(); +std::fill_n(&k_Hx_1_new[0][0], rows * cols, k_Hx_1[upml_width]); +std::fill_n(&k_Hx_2_new[0][0], rows * cols, k_Hx_2[upml_width]); +std::fill_n(&k_Ex_1_new[0][0], rows * cols, k_Ex_1[upml_width]); +std::fill_n(&k_Ex_2_new[0][0], rows * cols, k_Ex_2[upml_width]); + +// up down border +for(size_t i=0; i < upml_width; i += 1) { + for(size_t j=0; j < cols; j += 1) { + // up + k_Hx_1_new[i][j] = k_Hx_1[i]; + k_Hx_2_new[i][j] = k_Hx_2[i]; + + k_Ex_1_new[i][j] = k_Ex_1[i]; + k_Ex_2_new[i][j] = k_Ex_2[i]; + + //down + k_Hx_1_new[rows-i-1][j] = k_Hx_1[i]; + k_Hx_2_new[rows-i-1][j] = k_Hx_2[i]; + k_Ex_1_new[rows-i-1][j] = k_Ex_1[i]; + k_Ex_2_new[rows-i-1][j] = k_Ex_2[i]; + } +} + + + + +rows = k_Gy_1_new.size(); +cols = k_Gy_1_new[0].size(); +std::fill_n(&k_Gy_1_new[0][0], rows * cols, k_Gy_1[upml_width]); +std::fill_n(&k_Gy_2_new[0][0], rows * cols, k_Gy_2[upml_width]); + +// up down border +for(size_t i=0; i < upml_width; i += 1) { + for(size_t j=0; j < cols; j += 1) { + + // up + k_Gy_1_new[i][j] = k_Gy_1[i]; + k_Gy_2_new[i][j] = k_Gy_2[i]; + + //down + k_Gy_1_new[rows-i-1][j] = k_Gy_1[i]; + k_Gy_2_new[rows-i-1][j] = k_Gy_2[i]; + } +} + + + +rows = k_Hy_1_new.size(); +cols = k_Hy_1_new[0].size(); +std::fill_n(&k_Hy_1_new[0][0], rows * cols, k_Hy_1[upml_width]); +std::fill_n(&k_Hy_2_new[0][0], rows * cols, k_Hy_2[upml_width]); +std::fill_n(&k_Ey_1_new[0][0], rows * cols, k_Ey_1[upml_width]); +std::fill_n(&k_Ey_2_new[0][0], rows * cols, k_Ey_2[upml_width]); + +// left right border +for(size_t i=0; i < rows; i += 1) { + for(size_t j=0; j < upml_width; j += 1) { + // left + k_Hy_1_new[i][j] = k_Hy_1[j]; + k_Hy_2_new[i][j] = k_Hy_2[j]; + + k_Ey_1_new[i][j] = k_Ey_1[j]; + k_Ey_2_new[i][j] = k_Ey_2[j]; + + // right + k_Hy_1_new[i][cols-j-1] = k_Hy_1[j]; + k_Hy_2_new[i][cols-j-1] = k_Hy_2[j]; + + k_Ey_1_new[i][cols-j-1] = k_Ey_1[j]; + k_Ey_2_new[i][cols-j-1] = k_Ey_2[j]; + + } +} + + +rows = k_Hz_1_new.size(); +cols = k_Hz_1_new[0].size(); +std::fill_n(&k_Hz_1_new[0][0], rows * cols, k_Hz_1[upml_width]); +std::fill_n(&k_Hz_2_new[0][0], rows * cols, k_Hz_2[upml_width]); + +// up down border +for(size_t i=0; i < upml_width-1; i += 1) { + for(size_t j=0; j < cols; j += 1) { + + // up + k_Hz_1_new[i][j] = k_Hz_1[i+1]; + k_Hz_2_new[i][j] = k_Hz_2[i+1]; + + //down + k_Hz_1_new[rows-i-1][j] = k_Hz_1[i+1]; + k_Hz_2_new[rows-i-1][j] = k_Hz_2[i+1]; + } +} + + + +std::fill_n(&Index[0][0], nx * ny, 0.0); // C +std::fill_n(&IndexX[0][0], nx * ny, 0.0); // C +std::fill_n(&IndexY[0][0], nx * ny, 0.0); // C +// % User material config coordinates in cells +// % [i, cells_i, j, cells_j] +// mat_conf = [ +// [160, 9, 240, 20], +// [180, 15, 240, 20], +// [200, 10, 240, 20], +// [220, 10, 240, 20], +// [240, 10, 240, 20], +// [260, 10, 240, 20], +// [280, 10, 240, 20], +// [300, 10, 240, 20], +// [320, 10, 240, 20] +// ]; + +// mat_count = size(mat_conf, 1); +// % ///////////////////////////////////////////////////// + + +// Fill geometry matrix for 1/2 cylinder (no vectorization here), but for TE +// mode we need three different Index for Ex, Ey and Hz due to leapfrog +for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny; j += 1) { + for(size_t mat_idx=0; mat_idx < mat_count; mat_idx += 1) { + + //conf = mat_conf(mat_idx); + size_t x = mat_conf[mat_idx][0]; + size_t x_cells = mat_conf[mat_idx][1]; + size_t y = mat_conf[mat_idx][2]; + size_t y_cells = mat_conf[mat_idx][3]; + + if ((i > x) && (i < (x + x_cells)) + && (j > y) && (j < (y + y_cells))) { + Index[i][j] = 1; + } + if ((i+0.5 > x) && (i+0.5 < (x + x_cells)) + && (j > y) && (j < (y + y_cells))) { + IndexX[i][j] = 1; + } + if ((i > x) && (i < (x + x_cells)) + && (j+0.5 > y) && (j+0.5 < (y + y_cells))) { + IndexY[i][j] = 1; + } + } + } +} + + +for(size_t i=0; i < upml_width; i += 1) { + + size_t end; + + double sigma_x = sigma_max * std::pow((double)(upml_width-i)/upml_width, m); + double ka_x = 1+(ka_max-1)*std::pow((double)(upml_width-i)/upml_width, m); + + end = k_Fz_a.size() - 1; + k_Fz_a[end-i] = (2*epsilon_0*ka_x-sigma_x*dt)/(2*epsilon_0*ka_x+sigma_x*dt); + + end = k_Fz_b.size() - 1; + k_Fz_b[end-i] = 2*epsilon_0*dt/(2*epsilon_0*ka_x+sigma_x*dt)/dx; + + end = k_Hz_a.size() - 1; + k_Hz_a[end-i] = (2*epsilon_0*ka_x - sigma_x*dt)/(2*epsilon_0*ka_x + sigma_x*dt); + + end = k_Hz_b.size() - 1; + k_Hz_b[end-i] = 2*epsilon_0*dt/(2*epsilon_0*ka_x+sigma_x*dt)/(mu_0*Material[0][1]*dx); + + sigma_x = sigma_max * std::pow((double)(upml_width-i-0.5)/upml_width, m); + ka_x = 1+(ka_max-1)*std::pow((double)(upml_width-i-0.5)/upml_width, m); + + end = k_Hy_a.size() - 1; + k_Hy_a[end-i] = (2*epsilon_0*ka_x - sigma_x*dt)/(2*epsilon_0*ka_x + sigma_x*dt); + + end = k_Hy_b.size() - 1; + k_Hy_b[end-i] = 2*epsilon_0*dt/(2*epsilon_0*ka_x + sigma_x*dt)/(mu_0*Material[0][1]*dx); + + end = k_Ey_a.size() - 1; + k_Ey_a[end-i] = (2*epsilon_0*ka_x - sigma_x*dt)/(2*epsilon_0*ka_x + sigma_x*dt); + + end = k_Ey_b.size() - 1; + k_Ey_b[end-i] = 2*dt/(2*epsilon_0*ka_x + sigma_x*dt)/(Material[0][0]*dx); +} + + + +// end + +// %% Another transformation coefficients + +for(size_t i=0; i < number_of_materials; i += 1) { + K_a[i] = (2*epsilon_0*Material[i][0] - + Material[i][2]*dt) / (2*epsilon_0*Material[i][0] + + Material[i][2]*dt); + + K_b[i] = 2*epsilon_0*Material[i][0] / + (2*epsilon_0*Material[i][0] + Material[i][2]*dt); + + K_a1[i] = 1/(2*epsilon_0*Material[i][0] + Material[i][2]*dt); + + K_b1[i] = 2*epsilon_0*Material[i][0] - Material[i][2]*dt; +} + + +rows = K_a_new.size(); +cols = K_a_new[0].size(); +std::fill_n(&K_a_new[0][0], rows * cols, 1); +std::fill_n(&K_b_new[0][0], rows * cols, 1); + + +// % Replace [1,2] -> [1,-1] +for(size_t i=0; i < nx-2; i += 1) { + for(size_t j=0; j < ny-2; j += 1) { + + if(K_a_new[i][j] == 2) { + K_a_new[i][j] = K_a[1]; + } + + if(K_b_new[i][j] == 2) { + K_b_new[i][j] = K_b[1]; + } + } +} + + +for(size_t i=0; i < nx-2; i += 1) { + for(size_t j=0; j < ny-2; j += 1) { + M0[i][j] = Material[Index[i+1][j+1]][0]; + M1[i][j] = Material[Index[i+1][j+1]][1]; + } +} + +for(size_t i=0; i < nx; i += 1) { + for(size_t j=0; j < ny-1; j += 1) { + M2[i][j] = Material[Index[i][j]][1]; + } +} + +for(size_t i=0; i < nx-1; i += 1) { + for(size_t j=0; j < ny; j += 1) { + M3[i][j] = Material[Index[i][j]][1]; + } +} + + + +for(size_t j=0; j < K_b.size(); j += 1) { + // std::cout << K_b[j] << " "; + } +// std::cout << "\n" << "K_b[j]" << "\n"; + +std::ofstream myfile; + myfile.open("k_Fz_1_new3.txt"); + if (myfile.is_open()) + { + myfile << "This is a line.\n"; + myfile << "This is another line.\n"; + + // for(size_t i=0; i < nx; i += 1) { + // for(size_t j=0; j < ny; j += 1) { + + // std::cout << Index[i][j] << " "; + // } + // std::cout << '\n'; + // } + +// std::cout << "filfe" << '\n'; + myfile.close(); + } + else std::cout << "Unable to open file"; +for(size_t i=0; i < M3.size(); i += 1) { + for(size_t j=0; j < M3[0].size(); j += 1) { + + // std::cout << M3[i][j] << " "; + } + // std::cout << '\n'; +} + + + + // std::cout << "ka_max" << ka_max << std::endl; + } diff --git a/src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.h b/src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.h new file mode 100644 index 0000000..87b812e --- /dev/null +++ b/src/fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.h @@ -0,0 +1,320 @@ +// 2D FDTD with UPML and TF/SF. + +// Based on code from book: +// k_Fz_1_new +#ifndef TF_SF_H +#define TF_SF_H + +#include +#include + +#include // std::fill_n +#include +#include +#include +#include + + +class TFSF { + // Matlab array index start + const size_t mmm = 0; + + // Physical constants + // Absolute vacuum permittivity + const double epsilon_0 = 8.854*1e-12; + + const double pi = 3.1415926535897932385; + + // Absolute vacuum permeability + const double mu_0 = 4 * pi * 1e-7; + + + // Light speed in vacuum + const double light_speed = 2.99792458e8; + + // Main parameters + // Calculation area length per x and y axes + static constexpr double area_width = 2.5; + static constexpr double area_height = 2.5; + + // Uniform grid points for x and y axes + // static const size_t nx = 500; + // static const size_t ny = 500; + // static const size_t nx = 30; + // static const size_t ny = 30; + static const size_t nx = 220; + static const size_t ny = 220; + + // Perfect match layer (upml_width) thickness in uniform grid cells + static const size_t upml_width = 10; + // static const size_t upml_width = 3; + + // End time [s] + double t_end = 15e-9; + + // Excitation source amplitude and frequency [Hz] + double E0 = 1.0; + double frequency = 2.0e+9; // 2 GHz + + + //Width of alinea between total field area and calculation area border + //(scattered field interface width) [m] + static constexpr double len_tf_dx = 0.5; + static constexpr double len_tf_dy = 0.5; + + // Grid calculation + static constexpr double dx = area_width / nx; //C + static constexpr double dy = area_width / ny; //C + + // % C++ ignore + // X[nx]; + // Y[ny]; + + // % C++ ignore + // for i = 1:nx + // X(i) = (i-1)*dx - area_width/2; + // end + + // % C++ ignore + // for i = 1:ny + // Y(i) = (i-1)*dy - area_height/2; + // end + + + // Time step from CFL condition + const double CFL_FACTOR = 0.99; + // const double dt = ( 1/light_speed/sqrt( 1/(dx^2) + 1/(dy^2) ) )*CFL_FACTOR; + const double dt = CFL_FACTOR / (light_speed * sqrt(std::pow(1 / dx, 2) + std::pow(1 / dy, 2))); // C + + + // number_of_iterations = round(t_end/dt); + + // Geometry matrix + // Fill with 0 in code below + std::array, nx> Index; // C + std::array, nx> IndexX; // C + std::array, nx> IndexY; // C + + // Calculate size of total field area in TF/SF formalism + // static constexpr size_t nx_a = round(len_tf_dx/dx); + // static constexpr size_t nx_b = round((area_width-len_tf_dx)/dx); + // static constexpr size_t ny_a = round(len_tf_dy/dy); + // static constexpr size_t ny_b = round((area_height-len_tf_dy)/dy); + // static constexpr size_t nx_a = 100; + // static constexpr size_t nx_b = 400; + // static constexpr size_t ny_a = 100; + // static constexpr size_t ny_b = 400; + + // static constexpr size_t nx_a = 6; + // static constexpr size_t nx_b = 24; + // static constexpr size_t ny_a = 6; + // static constexpr size_t ny_b = 24; + + // static constexpr size_t nx_a = 60; + // static constexpr size_t nx_b = 160; + // static constexpr size_t ny_a = 60; + // static constexpr size_t ny_b = 160; + + // static constexpr size_t nx_a = 15; + // static constexpr size_t nx_b = 205; + // static constexpr size_t ny_a = 15; + // static constexpr size_t ny_b = 205; + + static constexpr size_t nx_a = 10; + static constexpr size_t nx_b = 209; + static constexpr size_t ny_a = 10; + static constexpr size_t ny_b = 209; + + // Pre-allocate 1D fields for TF/SF interface + // TM components + std::array Ez_1D; // C + std::array Fz_1D; // C + std::array Hy_1D; // C + + std::array k_Fz_a; // C + std::array k_Fz_b; // C + std::array k_Hy_a; // C + std::array k_Hy_b; // C + + // TE components + std::array Hz_1D; // C + std::array Ey_1D; // C + + + std::array k_Hz_a; // C half + std::array k_Hz_b; // C half + std::array k_Ey_a; // C + std::array k_Ey_b; // C + + const double ka_max = 1; + const double m = 4; + const double R_err = 1e-16; + + + // auxiliary + std::array, nx-2> k_Fz_1_new; // C + std::array, nx-2> k_Fz_2_new; // C + std::array, nx-2> k_Ez_1_new; // C + std::array, nx-2> k_Ez_2_new; // C + std::array, nx> k_Gx_1_new; // C + std::array, nx> k_Gx_2_new; // C + std::array, nx> k_Hx_1_new; // C + std::array, nx> k_Hx_2_new; // C + std::array, nx> k_Ex_1_new; // C + std::array, nx> k_Ex_2_new; // C + + std::array, nx-1> k_Gy_1_new; // C + std::array, nx-1> k_Gy_2_new; // C half + + std::array, nx-1> k_Ey_1_new; // C + std::array, nx-1> k_Ey_2_new; // C half + std::array, nx-1> k_Hy_1_new; // C half + std::array, nx-1> k_Hy_2_new; // C half + + std::array, nx-2> k_Hz_1_new; // C + std::array, nx-2> k_Hz_2_new; // C half + + std::array, nx-2> M0; // c + std::array, nx-2> M1; // c + std::array, nx> M2; // c + std::array, nx-1> M3; // c + + + std::array Fz_1D_r; // C + std::array, nx-2> Fz_r1; + + std::array, nx-2> Fz_r; // + std::array, nx-2> Tz_r; + + std::array, nx-1> Gy_r; + std::array, nx> Gx_r; + std::array, nx> Gx_r1; + std::array, nx-1> Gy_r1; + + std::array, nx-2> K_a_new; + std::array, nx-2> K_b_new; + + // Allocate 2D arrays + // TM physical and auxiliary fields + std::array, nx> Fz; + std::array, nx> Tz; + std::array, nx> Gx; + std::array, nx-1> Gy; + std::array, nx> Ez; + std::array, nx> Hx; + std::array, nx-1> Hy; + + // TE physical and auxiliary fields + std::array, nx> Wz; + std::array, nx> Mx; + std::array, nx-1> My; + std::array, nx> Hz; + std::array, nx> Ex; + std::array, nx-1> Ey; + + // Allocate UPML FDTD 1D coefficient arrays + // TM coefficients + std::array k_Fz_1; // C + std::array k_Fz_2; // C + std::array k_Ez_1; // C + std::array k_Ez_2; // C + std::array k_Gx_1; // C + std::array k_Gx_2; // C + std::array k_Hx_1; // C + std::array k_Hx_2; // C + + std::array k_Gy_1; // C + std::array k_Gy_2; // C half + + std::array k_Hy_1; // C half + std::array k_Hy_2; // C half + + // TM coefficients + std::array k_Hz_1; // C half + std::array k_Hz_2; // C half + std::array k_Ex_1; // C + std::array k_Ex_2; // C + std::array k_Ey_1; // C + std::array k_Ey_2; // C half + + + // Number of materials (include vacuum backrgound) + static const size_t number_of_materials = 2; + //Materials matrix + // Material = zeros(number_of_materials,3); + + // % eps, mu, sigma + // % Background relative permittivity, relative permeability and absolute + // % conductivity + std::array, number_of_materials> Material = {{ + {1.0, 1.0, 0.0}, //vacuum + {3, 3, 9.0e+70} + }}; + + std::array K_a; // C + std::array K_a1; // C + std::array K_b; // C + std::array K_b1; // C + + // tmp here start + size_t mat_count = 3; + +// mat_conf = [ +// 80, 15, 100, 20; +// 100, 10, 100, 20; +// 120, 10, 100, 20; +// 140, 10, 100, 20; +// ]; + + // std::array, 4> mat_conf = {{ + // {80, 15, 100, 20}, + // {100, 10, 100, 20}, + // {120, 10, 100, 20}, + // {140, 10, 100, 20} + + // }}; + + std::array, 3> mat_conf = {{ + {100, 25, 10, 85}, + {100, 25, 105, 10}, + {100, 25, 125, 80}, + // {120, 10, 100, 20}, + // {140, 10, 100, 20} + + }}; + + // tmp here end + + double k_Ez_a; + double k_Ez_b; + + double eta; + double sigma_max; + + size_t T = 0; + +public: + TFSF(); + void SetParams(); + void CalcNextLayer(); + size_t GetCurrentTimeStep(); + + struct Output { + // std::array, nx> Ez; + // std::array, nx> Hz; + std::array Ez; + std::array Hz; + std::array X; + std::array Y; + size_t rows; + size_t cols; + double maxEz; + double minEz; + double maxHz; + double minHz; + }; + + struct Output GetValues(); +}; + +#endif \ No newline at end of file diff --git a/src/index.cpp b/src/index.cpp index de60b02..f6857b6 100644 --- a/src/index.cpp +++ b/src/index.cpp @@ -4,12 +4,15 @@ #include "./transform-to-js/fdtd-1d/fdtd-1d.h" #include "./transform-to-js/fdtd-2d/fdtd-2d.h" +#include "./transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.h" // Callback method when module is registered with Node.js. Napi::Object Init(Napi::Env env, Napi::Object exports) { Napi::Object new_exports = exports; new_exports = Fdtd1D::Init(env, new_exports); new_exports = Fdtd2D::Init(env, new_exports); + new_exports = Fdtd2DTFSF::Init(env, new_exports); + return Fdtd1D::Init(env, new_exports); } diff --git a/src/transform-to-js/fdtd-2d copy/fdtd-2d.cpp b/src/transform-to-js/fdtd-2d copy/fdtd-2d.cpp new file mode 100644 index 0000000..ac2fb67 --- /dev/null +++ b/src/transform-to-js/fdtd-2d copy/fdtd-2d.cpp @@ -0,0 +1,323 @@ +// #include +// #include +// #include "fdtd-2d.h" + +// Napi::Object Fdtd2D::Init(Napi::Env env, Napi::Object exports) { +// // This method is used to hook the accessor and method callbacks. +// Napi::Function func = +// DefineClass(env, +// "Fdtd2D", +// { +// InstanceMethod("getNextTimeLayer", &Fdtd2D::GetNextTimeLayer), +// }); + +// Napi::FunctionReference* constructor = new Napi::FunctionReference(); + +// // Create a persistent reference to the class constructor. This will allow +// // a function called on a class prototype and a function +// // called on instance of a class to be distinguished from each other. +// *constructor = Napi::Persistent(func); +// exports.Set("Fdtd2D", func); + +// // Store the constructor as the add-on instance data. This will allow this +// // add-on to support multiple instances of itself running on multiple worker +// // threads, as well as multiple instances of itself running in different +// // contexts on the same thread. +// // +// // By default, the value set on the environment here will be destroyed when +// // the add-on is unloaded using the `delete` operator, but it is also +// // possible to supply a custom deleter. +// env.SetInstanceData(constructor); + +// return exports; +// } + +// Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info) +// : Napi::ObjectWrap(info), fdtd() { +// Napi::Env env = info.Env(); + +// // Object argument: +// // { +// // lambda, +// // beamsize, +// // isReload, +// // materialVector, // material matrix(flatten) +// // eps, +// // mu, +// // sigma, +// // dataReturnType, // ('Ez' = 0 | 'Hy' = 1 |'Hx' = 2 |'Energy' = 3) +// // srcPosition +// // } + +// if (info.Length() != 1) { +// Napi::TypeError::New(env, "Wrong number of arguments") +// .ThrowAsJavaScriptException(); +// return; +// } + +// if (!info[0].IsObject()) { +// Napi::TypeError::New(env, "Wrong argument! Should be an object.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// Napi::Object argObj = info[0].As().ToObject(); + +// // lambda +// if (!argObj.Has("lambda") || !argObj.Get("lambda").IsNumber()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'lambda' property and it should be number.") +// .ThrowAsJavaScriptException(); +// return; +// } +// double lambda = (double)argObj.Get("lambda").As(); + +// // beamsize +// if (!argObj.Has("beamsize") || !argObj.Get("beamsize").IsNumber()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'beamsize' property and it should be number.") +// .ThrowAsJavaScriptException(); +// return; +// } +// double beamsize = (double)argObj.Get("beamsize").As(); + +// // Reload params checker. +// if (!argObj.Has("isReload") || !argObj.Get("isReload").IsBoolean()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'isReload' property and it should be boolean.") +// .ThrowAsJavaScriptException(); +// return; +// } +// bool reload_check = static_cast(argObj.Get("isReload").As()); + +// // Material matrix(flatten). +// if (!argObj.Has("materialVector") || !argObj.Get("materialVector").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'materialVector' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // Material matrix transformation JS -> C++. +// const Napi::Array material_matrix_js = argObj.Get("materialVector").As(); +// int material_matrix_size = std::sqrt(material_matrix_js.Length()); + +// // Eps vector. +// if (!argObj.Has("eps") || !argObj.Get("eps").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'eps' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // eps transformation JS -> C++. +// const Napi::Array eps_js = argObj.Get("eps").As(); + +// // mu vector. +// if (!argObj.Has("mu") || !argObj.Get("mu").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'mu' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // mu transformation JS -> C++. +// const Napi::Array mu_js = argObj.Get("mu").As(); + +// // sigma vector. +// if (!argObj.Has("sigma") || !argObj.Get("sigma").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'sigma' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // sigma transformation JS -> C++. +// const Napi::Array sigma_js = argObj.Get("sigma").As(); + +// // Data return type('Ez' = 0 | 'Hy' = 1 |'Hx' = 2 |'Energy' = 3) +// if (!argObj.Has("dataReturnType") || !argObj.Get("dataReturnType").IsNumber()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'dataReturnType' property and it should be number(integer).") +// .ThrowAsJavaScriptException(); +// return; +// } + +// int data_return_type = static_cast(argObj.Get("dataReturnType").As()); +// this->data_return_type = data_return_type; + +// std::cout << "before ";// << material_matrix_size; + +// // Temporary matrix. +// std::vector> temp_matrix; + +// std::cout << "--" << material_matrix_size << "--"; + +// // Transform input flatten matrix into 2-dimensional matrix. +// for (int i = 0; i < material_matrix_size; i++) { +// temp_matrix.push_back(std::vector()); + +// for (int j = 0; j < material_matrix_size; j++) { +// temp_matrix[i].push_back( +// (int)material_matrix_js[i * material_matrix_size + j] +// .As()); +// std::cout << ", " << i; +// } +// } + +// std::cout << "after"; + +// const size_t rows = FdtdPml2D::GetRows(); +// const size_t cols = FdtdPml2D::GetCols(); + +// // Matrix size coefficient. +// size_t coeff = rows / material_matrix_size; + +// // Initialization eps, mu, sigma matrixes. +// std::vector> eps_matrix; +// std::vector> mu_matrix; +// std::vector> sigma_matrix; +// for (int i = 0; i < rows; i++) { +// eps_matrix.push_back(std::vector()); +// mu_matrix.push_back(std::vector()); +// sigma_matrix.push_back(std::vector()); +// for (int j = 0; j < cols; j++) { +// eps_matrix[i].push_back(0); +// mu_matrix[i].push_back(0); +// sigma_matrix[i].push_back(0); +// } +// } + +// // Filling eps, mu, sigma matrixes. +// for (int i = 0; i < material_matrix_size; i++) { +// for (int j = 0; j < material_matrix_size; j++) { +// for (int k = 0; k < coeff; k++) { +// for (int f = 0; f < coeff; f++) { +// int index = temp_matrix[i][j]; + +// // Rotate matrix on 90 degree for correctness in numerical method. +// // eps_matrix[j * coeff + f][i * coeff + k] = static_cast(eps_js[index].As()); +// // mu_matrix[j * coeff + f][i * coeff + k] = static_cast(mu_js[index].As()); +// // sigma_matrix[j * coeff + f][i * coeff + k] = static_cast(sigma_js[index].As()); + +// // Without rotate. +// eps_matrix[i * coeff + k][j * coeff + f] = static_cast(eps_js[index].As()); +// mu_matrix[i * coeff + k][j * coeff + f] = static_cast(mu_js[index].As()); +// sigma_matrix[i * coeff + k][j * coeff + f] = static_cast(sigma_js[index].As()); +// } +// } +// } +// } + +// // srcPositionArray. +// if (!argObj.Has("srcPosition") || !argObj.Get("srcPosition").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'srcPosition' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // srcPosition transformation JS -> C++. +// const Napi::Array src_position_array = argObj.Get("srcPosition").As(); +// double relative_src_position_x = (float)src_position_array[(uint32_t)0].As(); +// double relative_src_position_y = (float)src_position_array[(uint32_t)1].As(); + +// size_t src_position_row = static_cast(relative_src_position_y * rows); +// size_t src_position_col = static_cast(relative_src_position_x * cols); + +// fdtd.SetParams(eps_matrix, mu_matrix, sigma_matrix, src_position_row, src_position_col); +// } + +// // Fdtd method in 1D case. +// Napi::Value Fdtd2D::GetNextTimeLayer(const Napi::CallbackInfo& info) { +// Napi::Env env = info.Env(); + +// std::vector vect_X = {}; +// std::vector vect_Y = {}; +// std::vector vect_Ez = {}; +// std::vector vect_Hy = {}; +// std::vector vect_Hx = {}; +// std::vector vect_Energy = {}; + +// double max = 0.001; +// double min = -0.001; + +// this->fdtd.CalcNextLayer(vect_X, vect_Y, vect_Ez, vect_Hy, vect_Hx, vect_Energy, max, min); + +// const size_t rows = FdtdPml2D::GetRows(); +// const size_t cols = FdtdPml2D::GetCols(); + +// // Matrix sizes. +// size_t client_rows = rows / fdtd.GetStep(); +// size_t client_cols = cols / fdtd.GetStep(); + +// const size_t js_arrays_size = client_rows * client_cols; + +// // Creating JS arrays to store C++ arrays. +// Napi::Array js_data_X = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Y = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Ez = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Hy = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Hx = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Energy = Napi::Array::New(env, js_arrays_size); + +// // Filling JS arrays with C++ arrays data. +// for (size_t i = 0; i < js_arrays_size; i++) { +// js_data_X[i] = Napi::Number::New(env, vect_X[i]); +// js_data_Y[i] = Napi::Number::New(env, vect_Y[i]); +// js_data_Ez[i] = Napi::Number::New(env, vect_Ez[i]); +// js_data_Hy[i] = Napi::Number::New(env, vect_Hy[i]); +// js_data_Hx[i] = Napi::Number::New(env, vect_Hx[i]); +// js_data_Energy[i] = Napi::Number::New(env, vect_Energy[i]); +// } + +// // Creating JS object to return. +// Napi::Object data = Napi::Array::New(env); +// data.Set("dataX", js_data_X); +// data.Set("dataY", js_data_Y); +// data.Set("rows", client_rows); +// data.Set("cols", client_cols); +// data.Set("timeStep", fdtd.GetCurrentTimeStep()); + +// switch (this->data_return_type) { +// case 0: +// data.Set("dataEz", js_data_Ez); +// // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); +// // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); +// break; +// case 1: +// data.Set("dataHy", js_data_Hy); +// // max = *std::max_element(std::begin(vect_Hy), std::end(vect_Hy)); +// // min = *std::min_element(std::begin(vect_Hy), std::end(vect_Hy)); +// break; +// case 2: +// data.Set("dataHx", js_data_Hx); +// // max = *std::max_element(std::begin(vect_Hx), std::end(vect_Hx)); +// // min = *std::min_element(std::begin(vect_Hx), std::end(vect_Hx)); +// break; +// case 3: +// data.Set("dataEnergy", js_data_Energy); +// // max = *std::max_element(std::begin(vect_Energy), std::end(vect_Energy)); +// // min = *std::min_element(std::begin(vect_Energy), std::end(vect_Energy)); +// break; + +// default: +// // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); +// // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); +// break; +// } +// // Fill max and min values. +// data.Set("max", max); +// data.Set("min", min); + +// return data; +// } \ No newline at end of file diff --git a/src/transform-to-js/fdtd-2d copy/fdtd-2d.h b/src/transform-to-js/fdtd-2d copy/fdtd-2d.h new file mode 100644 index 0000000..3a41e6f --- /dev/null +++ b/src/transform-to-js/fdtd-2d copy/fdtd-2d.h @@ -0,0 +1,20 @@ +// #ifndef FDTD_2D_H +// #define FDTD_2D_H + +// #include + +// #include "../../fdtd/2d-pml/fdtd-pml-2d.h" + +// class Fdtd2D : public Napi::ObjectWrap { +// public: +// static Napi::Object Init(Napi::Env env, Napi::Object exports); +// Fdtd2D(const Napi::CallbackInfo& info); + +// private: +// Napi::Value GetNextTimeLayer(const Napi::CallbackInfo& info); +// FdtdPml2D fdtd; +// // Data return type('Ez' = 0 | 'Hy' = 1 |'Hx' = 2 |'Energy' = 3) +// int data_return_type = 0; +// }; + +// #endif \ No newline at end of file diff --git a/src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.cpp b/src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.cpp new file mode 100644 index 0000000..d191aaa --- /dev/null +++ b/src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.cpp @@ -0,0 +1,352 @@ +#include +#include +#include "fdtd-2d-tf-sf.h" + +Napi::Object Fdtd2DTFSF::Init(Napi::Env env, Napi::Object exports) { + // This method is used to hook the accessor and method callbacks. + Napi::Function func = + DefineClass(env, + "Fdtd2DTFSF", + { + InstanceMethod("getNextTimeLayer", &Fdtd2DTFSF::GetNextTimeLayer), + } + ); + + Napi::FunctionReference* constructor = new Napi::FunctionReference(); + + // Create a persistent reference to the class constructor. This will allow + // a function called on a class prototype and a function + // called on instance of a class to be distinguished from each other. + *constructor = Napi::Persistent(func); + exports.Set("Fdtd2DTFSF", func); + + // Store the constructor as the add-on instance data. This will allow this + // add-on to support multiple instances of itself running on multiple worker + // threads, as well as multiple instances of itself running in different + // contexts on the same thread. + // + // By default, the value set on the environment here will be destroyed when + // the add-on is unloaded using the `delete` operator, but it is also + // possible to supply a custom deleter. + env.SetInstanceData(constructor); + + return exports; +} + +Fdtd2DTFSF::Fdtd2DTFSF(const Napi::CallbackInfo& info) + : Napi::ObjectWrap(info), fdtd() { + Napi::Env env = info.Env(); + + // Object argument: + // { + // lambda, + // beamsize, + // isReload, + // materialVector, // material matrix(flatten) + // eps, + // mu, + // sigma, + // dataReturnType, // ('Ez' = 0 | 'Hy' = 1 |'Hx' = 2 |'Energy' = 3) + // srcPosition + // } + +// if (info.Length() != 1) { +// Napi::TypeError::New(env, "Wrong number of arguments") +// .ThrowAsJavaScriptException(); +// return; +// } + +// if (!info[0].IsObject()) { +// Napi::TypeError::New(env, "Wrong argument! Should be an object.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// Napi::Object argObj = info[0].As().ToObject(); + +// // lambda +// if (!argObj.Has("lambda") || !argObj.Get("lambda").IsNumber()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'lambda' property and it should be number.") +// .ThrowAsJavaScriptException(); +// return; +// } +// double lambda = (double)argObj.Get("lambda").As(); + +// // beamsize +// if (!argObj.Has("beamsize") || !argObj.Get("beamsize").IsNumber()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'beamsize' property and it should be number.") +// .ThrowAsJavaScriptException(); +// return; +// } +// double beamsize = (double)argObj.Get("beamsize").As(); + +// // Reload params checker. +// if (!argObj.Has("isReload") || !argObj.Get("isReload").IsBoolean()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'isReload' property and it should be boolean.") +// .ThrowAsJavaScriptException(); +// return; +// } +// bool reload_check = static_cast(argObj.Get("isReload").As()); + +// // Material matrix(flatten). +// if (!argObj.Has("materialVector") || !argObj.Get("materialVector").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'materialVector' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // Material matrix transformation JS -> C++. +// const Napi::Array material_matrix_js = argObj.Get("materialVector").As(); +// int material_matrix_size = std::sqrt(material_matrix_js.Length()); + +// // Eps vector. +// if (!argObj.Has("eps") || !argObj.Get("eps").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'eps' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // eps transformation JS -> C++. +// const Napi::Array eps_js = argObj.Get("eps").As(); + +// // mu vector. +// if (!argObj.Has("mu") || !argObj.Get("mu").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'mu' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // mu transformation JS -> C++. +// const Napi::Array mu_js = argObj.Get("mu").As(); + +// // sigma vector. +// if (!argObj.Has("sigma") || !argObj.Get("sigma").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'sigma' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // sigma transformation JS -> C++. +// const Napi::Array sigma_js = argObj.Get("sigma").As(); + +// // Data return type('Ez' = 0 | 'Hy' = 1 |'Hx' = 2 |'Energy' = 3) +// if (!argObj.Has("dataReturnType") || !argObj.Get("dataReturnType").IsNumber()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'dataReturnType' property and it should be number(integer).") +// .ThrowAsJavaScriptException(); +// return; +// } + +// int data_return_type = static_cast(argObj.Get("dataReturnType").As()); +// this->data_return_type = data_return_type; + +// std::cout << "before ";// << material_matrix_size; + +// // Temporary matrix. +// std::vector> temp_matrix; + +// std::cout << "--" << material_matrix_size << "--"; + +// // Transform input flatten matrix into 2-dimensional matrix. +// for (int i = 0; i < material_matrix_size; i++) { +// temp_matrix.push_back(std::vector()); + +// for (int j = 0; j < material_matrix_size; j++) { +// temp_matrix[i].push_back( +// (int)material_matrix_js[i * material_matrix_size + j] +// .As()); +// std::cout << ", " << i; +// } +// } + +// std::cout << "after"; + +// const size_t rows = FdtdPml2D::GetRows(); +// const size_t cols = FdtdPml2D::GetCols(); + +// // Matrix size coefficient. +// size_t coeff = rows / material_matrix_size; + +// // Initialization eps, mu, sigma matrixes. +// std::vector> eps_matrix; +// std::vector> mu_matrix; +// std::vector> sigma_matrix; +// for (int i = 0; i < rows; i++) { +// eps_matrix.push_back(std::vector()); +// mu_matrix.push_back(std::vector()); +// sigma_matrix.push_back(std::vector()); +// for (int j = 0; j < cols; j++) { +// eps_matrix[i].push_back(0); +// mu_matrix[i].push_back(0); +// sigma_matrix[i].push_back(0); +// } +// } + +// // Filling eps, mu, sigma matrixes. +// for (int i = 0; i < material_matrix_size; i++) { +// for (int j = 0; j < material_matrix_size; j++) { +// for (int k = 0; k < coeff; k++) { +// for (int f = 0; f < coeff; f++) { +// int index = temp_matrix[i][j]; + +// // Rotate matrix on 90 degree for correctness in numerical method. +// // eps_matrix[j * coeff + f][i * coeff + k] = static_cast(eps_js[index].As()); +// // mu_matrix[j * coeff + f][i * coeff + k] = static_cast(mu_js[index].As()); +// // sigma_matrix[j * coeff + f][i * coeff + k] = static_cast(sigma_js[index].As()); + +// // Without rotate. +// eps_matrix[i * coeff + k][j * coeff + f] = static_cast(eps_js[index].As()); +// mu_matrix[i * coeff + k][j * coeff + f] = static_cast(mu_js[index].As()); +// sigma_matrix[i * coeff + k][j * coeff + f] = static_cast(sigma_js[index].As()); +// } +// } +// } +// } + +// // srcPositionArray. +// if (!argObj.Has("srcPosition") || !argObj.Get("srcPosition").IsArray()) { +// Napi::TypeError::New( +// env, +// "Object should contain 'srcPosition' property and it should be array.") +// .ThrowAsJavaScriptException(); +// return; +// } + +// // srcPosition transformation JS -> C++. +// const Napi::Array src_position_array = argObj.Get("srcPosition").As(); +// double relative_src_position_x = (float)src_position_array[(uint32_t)0].As(); +// double relative_src_position_y = (float)src_position_array[(uint32_t)1].As(); + +// size_t src_position_row = static_cast(relative_src_position_y * rows); +// size_t src_position_col = static_cast(relative_src_position_x * cols); + +// fdtd.SetParams(eps_matrix, mu_matrix, sigma_matrix, src_position_row, src_position_col); +} + +// Fdtd method in 1D case. +Napi::Value Fdtd2DTFSF::GetNextTimeLayer(const Napi::CallbackInfo& info) { + Napi::Env env = info.Env(); + +// std::vector vect_X = {}; +// std::vector vect_Y = {}; +// std::vector vect_Ez = {}; +// std::vector vect_Hy = {}; +// std::vector vect_Hx = {}; +// std::vector vect_Energy = {}; + +// double max = 0.001; +// double min = -0.001; + + // this->fdtd.CalcNextLayer(vect_X, vect_Y, vect_Ez, vect_Hy, vect_Hx, vect_Energy, max, min); + this->fdtd.CalcNextLayer(); + + struct TFSF::Output fdtd_output = this->fdtd.GetValues(); + size_t rows = fdtd_output.rows; + size_t cols = fdtd_output.cols; + + // std::cout << rows << '\n'; + // std::cout << cols << '\n'; + +// for(size_t i=0; i < fields.Ez.size(); i += 1) { +// for(size_t j=0; j < fields.Ez[0].size(); j += 1) { + +// std::cout << fields.Ez[i][j] << " "; +// } +// std::cout << '\n'; +// } + +// const size_t rows = FdtdPml2D::GetRows(); +// const size_t cols = FdtdPml2D::GetCols(); + +// // Matrix sizes. +// size_t client_rows = rows / fdtd.GetStep(); +// size_t client_cols = cols / fdtd.GetStep(); + +// const size_t js_arrays_size = client_rows * client_cols; + const size_t js_arrays_size = rows * cols; + + // Creating JS arrays to store C++ arrays. + Napi::Array js_data_X = Napi::Array::New(env, js_arrays_size); + Napi::Array js_data_Y = Napi::Array::New(env, js_arrays_size); + Napi::Array js_data_Ez = Napi::Array::New(env, js_arrays_size); + Napi::Array js_data_Hz = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Hy = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Hx = Napi::Array::New(env, js_arrays_size); +// Napi::Array js_data_Energy = Napi::Array::New(env, js_arrays_size); + +// // Filling JS arrays with C++ arrays data. + for (size_t i = 0; i < js_arrays_size; i++) { + js_data_X[i] = Napi::Number::New(env, fdtd_output.X[i]); + js_data_Y[i] = Napi::Number::New(env, fdtd_output.Y[i]); + js_data_Ez[i] = Napi::Number::New(env, fdtd_output.Ez[i]); + js_data_Hz[i] = Napi::Number::New(env, fdtd_output.Hz[i]); +// js_data_Hy[i] = Napi::Number::New(env, vect_Hy[i]); +// js_data_Hx[i] = Napi::Number::New(env, vect_Hx[i]); +// js_data_Energy[i] = Napi::Number::New(env, vect_Energy[i]); + } + + // Creating JS object to return. + Napi::Object data = Napi::Array::New(env); + data.Set("dataX", js_data_X); + data.Set("dataY", js_data_Y); + data.Set("rows", rows); + data.Set("cols", cols); + data.Set("dataEz", js_data_Ez); + data.Set("dataHz", js_data_Hz); + data.Set("timeStep", fdtd.GetCurrentTimeStep()); + + // data.Set("maxEz", fdtd_output.maxEz); + // data.Set("minEz", fdtd_output.minEz); + data.Set("max", fdtd_output.maxEz); + data.Set("min", fdtd_output.minEz); + data.Set("maxHz", fdtd_output.maxHz); + data.Set("minHz", fdtd_output.minHz); + + // switch (this->data_return_type) { + // case 0: + // data.Set("dataEz", js_data_Ez); + // // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); + // // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); + // break; + // case 1: + // data.Set("dataHy", js_data_Hy); + // // max = *std::max_element(std::begin(vect_Hy), std::end(vect_Hy)); + // // min = *std::min_element(std::begin(vect_Hy), std::end(vect_Hy)); + // break; + // case 2: + // data.Set("dataHx", js_data_Hx); + // // max = *std::max_element(std::begin(vect_Hx), std::end(vect_Hx)); + // // min = *std::min_element(std::begin(vect_Hx), std::end(vect_Hx)); + // break; + // case 3: + // data.Set("dataEnergy", js_data_Energy); + // // max = *std::max_element(std::begin(vect_Energy), std::end(vect_Energy)); + // // min = *std::min_element(std::begin(vect_Energy), std::end(vect_Energy)); + // break; + + // default: + // // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); + // // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); + // break; + // } + // // Fill max and min values. + // data.Set("max", max); + // data.Set("min", min); + + return data; +} \ No newline at end of file diff --git a/src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.h b/src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.h new file mode 100644 index 0000000..580eaff --- /dev/null +++ b/src/transform-to-js/fdtd-2d-tf-sf/fdtd-2d-tf-sf.h @@ -0,0 +1,20 @@ +#ifndef FDTD_2D_TF_SF_H +#define FDTD_2D_TF_SF_H + +#include + +#include "../../fdtd/2d-upml-tf-sf/fdtd-2d-upml-tf-sf.h" + +class Fdtd2DTFSF : public Napi::ObjectWrap { + public: + static Napi::Object Init(Napi::Env env, Napi::Object exports); + Fdtd2DTFSF(const Napi::CallbackInfo& info); + + private: + Napi::Value GetNextTimeLayer(const Napi::CallbackInfo& info); + TFSF fdtd; + // Data return type('Ez' = 0 | 'Hy' = 1 |'Hx' = 2 |'Energy' = 3) + // int data_return_type = 0; +}; + +#endif \ No newline at end of file diff --git a/src/transform-to-js/fdtd-2d/fdtd-2d.cpp b/src/transform-to-js/fdtd-2d/fdtd-2d.cpp index cab10cf..c452e93 100644 --- a/src/transform-to-js/fdtd-2d/fdtd-2d.cpp +++ b/src/transform-to-js/fdtd-2d/fdtd-2d.cpp @@ -1,3 +1,5 @@ +#include +#include #include "fdtd-2d.h" Napi::Object Fdtd2D::Init(Napi::Env env, Napi::Object exports) { @@ -102,7 +104,7 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info) // Material matrix transformation JS -> C++. const Napi::Array material_matrix_js = argObj.Get("materialVector").As(); - int material_matrix_size = material_matrix_js.Length() / 2; + int material_matrix_size = std::sqrt(material_matrix_js.Length()); // Eps vector. if (!argObj.Has("eps") || !argObj.Get("eps").IsArray()) { @@ -140,9 +142,6 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info) // sigma transformation JS -> C++. const Napi::Array sigma_js = argObj.Get("sigma").As(); - // Temporary matrix. - std::vector> temp_matrix; - // Data return type('Ez' = 0 | 'Hy' = 1 |'Hx' = 2 |'Energy' = 3) if (!argObj.Has("dataReturnType") || !argObj.Get("dataReturnType").IsNumber()) { Napi::TypeError::New( @@ -155,18 +154,31 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info) int data_return_type = static_cast(argObj.Get("dataReturnType").As()); this->data_return_type = data_return_type; + // std::cout << "before ";// << material_matrix_size; + + // Temporary matrix. + std::vector> temp_matrix; + + // std::cout << "--" << material_matrix_size << "--"; + // Transform input flatten matrix into 2-dimensional matrix. for (int i = 0; i < material_matrix_size; i++) { temp_matrix.push_back(std::vector()); + for (int j = 0; j < material_matrix_size; j++) { temp_matrix[i].push_back( (int)material_matrix_js[i * material_matrix_size + j] .As()); + // std::cout << ", " << i; } } - const size_t rows = FdtdPml2D::GetRows(); - const size_t cols = FdtdPml2D::GetCols(); + // std::cout << "after"; + + // const size_t rows = FdtdPml2D::GetRows(); + // const size_t cols = FdtdPml2D::GetCols(); + const size_t rows = 220; + const size_t cols = 220; // Matrix size coefficient. size_t coeff = rows / material_matrix_size; @@ -224,7 +236,8 @@ Fdtd2D::Fdtd2D(const Napi::CallbackInfo& info) size_t src_position_row = static_cast(relative_src_position_y * rows); size_t src_position_col = static_cast(relative_src_position_x * cols); - fdtd.SetParams(eps_matrix, mu_matrix, sigma_matrix, src_position_row, src_position_col); + // fdtd.SetParams(eps_matrix, mu_matrix, sigma_matrix, src_position_row, src_position_col); + fdtd.InitializeFdtd(); } // Fdtd method in 1D case. @@ -233,41 +246,53 @@ Napi::Value Fdtd2D::GetNextTimeLayer(const Napi::CallbackInfo& info) { std::vector vect_X = {}; std::vector vect_Y = {}; - std::vector vect_Ez = {}; - std::vector vect_Hy = {}; - std::vector vect_Hx = {}; - std::vector vect_Energy = {}; + std::vector vect_Hz = {}; + // std::vector vect_Ez = {}; + // std::vector vect_Hy = {}; + // std::vector vect_Hx = {}; + // std::vector vect_Energy = {}; double max = 0.001; double min = -0.001; - this->fdtd.CalcNextLayer(vect_X, vect_Y, vect_Ez, vect_Hy, vect_Hx, vect_Energy, max, min); + // this->fdtd.CalcNextLayer(vect_X, vect_Y, vect_Ez, vect_Hy, vect_Hx, vect_Energy, max, min); + this->fdtd.CalcNextLayer(); + + struct FdtdPml2D::Output fdtd_output = this->fdtd.GetValues(); + size_t rows = fdtd_output.rows; + size_t cols = fdtd_output.cols; - const size_t rows = FdtdPml2D::GetRows(); - const size_t cols = FdtdPml2D::GetCols(); + // const size_t rows = FdtdPml2D::GetRows(); + // const size_t cols = FdtdPml2D::GetCols(); + // const size_t rows = 100; + // const size_t cols = 100; // Matrix sizes. - size_t client_rows = rows / fdtd.GetStep(); - size_t client_cols = cols / fdtd.GetStep(); + // size_t client_rows = rows / fdtd.GetStep(); + // size_t client_cols = cols / fdtd.GetStep(); + size_t client_rows = rows; + size_t client_cols = cols; const size_t js_arrays_size = client_rows * client_cols; // Creating JS arrays to store C++ arrays. Napi::Array js_data_X = Napi::Array::New(env, js_arrays_size); Napi::Array js_data_Y = Napi::Array::New(env, js_arrays_size); - Napi::Array js_data_Ez = Napi::Array::New(env, js_arrays_size); - Napi::Array js_data_Hy = Napi::Array::New(env, js_arrays_size); - Napi::Array js_data_Hx = Napi::Array::New(env, js_arrays_size); - Napi::Array js_data_Energy = Napi::Array::New(env, js_arrays_size); + Napi::Array js_data_Hz = Napi::Array::New(env, js_arrays_size); + // Napi::Array js_data_Ez = Napi::Array::New(env, js_arrays_size); + // Napi::Array js_data_Hy = Napi::Array::New(env, js_arrays_size); + // Napi::Array js_data_Hx = Napi::Array::New(env, js_arrays_size); + // Napi::Array js_data_Energy = Napi::Array::New(env, js_arrays_size); // Filling JS arrays with C++ arrays data. for (size_t i = 0; i < js_arrays_size; i++) { - js_data_X[i] = Napi::Number::New(env, vect_X[i]); - js_data_Y[i] = Napi::Number::New(env, vect_Y[i]); - js_data_Ez[i] = Napi::Number::New(env, vect_Ez[i]); - js_data_Hy[i] = Napi::Number::New(env, vect_Hy[i]); - js_data_Hx[i] = Napi::Number::New(env, vect_Hx[i]); - js_data_Energy[i] = Napi::Number::New(env, vect_Energy[i]); + js_data_X[i] = Napi::Number::New(env, fdtd_output.X[i]); + js_data_Y[i] = Napi::Number::New(env, fdtd_output.Y[i]); + js_data_Hz[i] = Napi::Number::New(env, fdtd_output.Hz[i]); + // js_data_Ez[i] = Napi::Number::New(env, vect_Ez[i]); + // js_data_Hy[i] = Napi::Number::New(env, vect_Hy[i]); + // js_data_Hx[i] = Napi::Number::New(env, vect_Hx[i]); + // js_data_Energy[i] = Napi::Number::New(env, vect_Energy[i]); } // Creating JS object to return. @@ -276,38 +301,41 @@ Napi::Value Fdtd2D::GetNextTimeLayer(const Napi::CallbackInfo& info) { data.Set("dataY", js_data_Y); data.Set("rows", client_rows); data.Set("cols", client_cols); + data.Set("dataEz", js_data_Hz); + data.Set("max", fdtd_output.maxHz); + data.Set("min", fdtd_output.minHz); data.Set("timeStep", fdtd.GetCurrentTimeStep()); - switch (this->data_return_type) { - case 0: - data.Set("dataEz", js_data_Ez); - // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); - // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); - break; - case 1: - data.Set("dataHy", js_data_Hy); - // max = *std::max_element(std::begin(vect_Hy), std::end(vect_Hy)); - // min = *std::min_element(std::begin(vect_Hy), std::end(vect_Hy)); - break; - case 2: - data.Set("dataHx", js_data_Hx); - // max = *std::max_element(std::begin(vect_Hx), std::end(vect_Hx)); - // min = *std::min_element(std::begin(vect_Hx), std::end(vect_Hx)); - break; - case 3: - data.Set("dataEnergy", js_data_Energy); - // max = *std::max_element(std::begin(vect_Energy), std::end(vect_Energy)); - // min = *std::min_element(std::begin(vect_Energy), std::end(vect_Energy)); - break; - - default: - // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); - // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); - break; - } - // Fill max and min values. - data.Set("max", max); - data.Set("min", min); - + // switch (this->data_return_type) { + // case 0: + // data.Set("dataEz", js_data_Ez); + // // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); + // // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); + // break; + // case 1: + // data.Set("dataHy", js_data_Hy); + // // max = *std::max_element(std::begin(vect_Hy), std::end(vect_Hy)); + // // min = *std::min_element(std::begin(vect_Hy), std::end(vect_Hy)); + // break; + // case 2: + // data.Set("dataHx", js_data_Hx); + // // max = *std::max_element(std::begin(vect_Hx), std::end(vect_Hx)); + // // min = *std::min_element(std::begin(vect_Hx), std::end(vect_Hx)); + // break; + // case 3: + // data.Set("dataEnergy", js_data_Energy); + // // max = *std::max_element(std::begin(vect_Energy), std::end(vect_Energy)); + // // min = *std::min_element(std::begin(vect_Energy), std::end(vect_Energy)); + // break; + + // default: + // // max = *std::max_element(std::begin(vect_Ez), std::end(vect_Ez)); + // // min = *std::min_element(std::begin(vect_Ez), std::end(vect_Ez)); + // break; + // } + // // Fill max and min values. + // data.Set("max", max); + // data.Set("min", min); + return data; } \ No newline at end of file diff --git a/test-addon.ts b/test-addon.ts index fef135f..185a497 100644 --- a/test-addon.ts +++ b/test-addon.ts @@ -33,7 +33,7 @@ const test1D = () => { const test2D = () => { const lambda = 1; const beamsize = 1; - const materialVector = [1, 0, 2, 0]; + const materialVector = [1, 0, 2, 0, 1, 1, 1,1 ,1]; const eps = [1.0, 1.2, 1.1]; const mu = [0.51, 0.5, 0.57]; const sigma = [1.0, 0.001, 1.0]; @@ -54,7 +54,7 @@ const test2D = () => { }); let data; - for (let j = 0; j < 49; ++j) { + for (let j = 0; j < 300; ++j) { data = fdtd.getNextTimeLayer(); } console.log(data); @@ -68,6 +68,22 @@ function testMemoryUsage() { } -test1D(); -// test2D(); +const test2DTFSF = () => { + + + let fdtd = new addon.Fdtd2DTFSF(); + + let data; + for (let j = 0; j < 30; ++j) { + data = fdtd.getNextTimeLayer(); + } + // console.log(fdtd); + console.log(data); +}; + + +// test1D(); +test2D(); + +// test2DTFSF(); // testMemoryUsage(); diff --git a/upml.cpp b/upml.cpp new file mode 100644 index 0000000..e69de29 diff --git a/upml_tf_sf.m b/upml_tf_sf.m new file mode 100644 index 0000000..7357fe4 --- /dev/null +++ b/upml_tf_sf.m @@ -0,0 +1,983 @@ + +%% FDTD 2D with UPML and TF/SF interface example. + + +%% Initialize workspace +clear all; close all; clc; format short; + +% Matlab array index start +mmm = 1; + +%% Fundamental physical constants +% Absolute vacuum permittivity +epsilon_0 = 8.854*1e-12; +% Absulute vacuum permeability +mu_0 = 4*pi*1e-7; +% Light speed in vacuum +light_speed = 1/sqrt(epsilon_0*mu_0); + +%% Main parameters +% Calculation area length per x and y axes +area_width = 2.5; +area_height = 2.5; + +% Uniform grid points for x and y axes +nx = 500; +ny = 500; + +% Perfect match layer (pml_width) thickness in uniform grid cells +pml_width = 10; + +% End time [s] +t_end = 15e-9; + +% Excitation source amplitude and frequency [Hz] +E0 = 1.0; +f = 2.0e+9; % 2 GHz + +% Number of materials (include vacuum backrgound) +number_of_materials = 2; +% Background relative permittivity, relative permeability and absolute +% conductivity +eps_back = 1.0; +mu_back = 1.0; +sig_back = 0.0; + +% Width of alinea between total field area and calculation area border +% (scattered field interface width) [m] +len_tf_dx = 0.5; +len_tf_dy = 0.5; + +% Grid calculation +% X = linspace(0,area_width,nx)-area_width/2; +% Y = linspace(0,area_height,ny)-area_height/2; +% dx = X(nx)-X(nx-1); +% dy = Y(ny)-Y(ny-1); + +dx = area_width / nx; +dy = area_width / ny; +X = []; +Y = []; +for i = 1:nx + X(i) = (i-1)*dx - area_width/2; + Y(i) = (i-1)*dy - area_height/2; +end + + +% Time step from CFL condition +CFL_FACTOR = 0.99; +dt = ( 1/light_speed/sqrt( 1/(dx^2) + 1/(dy^2) ) )*CFL_FACTOR; +number_of_iterations = round(t_end/dt); + +%% Geometry matrix +Index = zeros(nx,ny); +IndexX = zeros(nx,ny); +IndexY = zeros(nx,ny); + + +%% Materials matrix +Material = zeros(number_of_materials,3); +Material(1,1) = eps_back; +Material(1,2) = mu_back; +Material(1,3) = sig_back; + + +%% Simple geometry (here metal round cylinder 1/2) +% Diameter [m] +d0 = 0.4; +% Cylinder materials +Material(2,1) = 3; % relative permittivity +Material(2,2) = 3; % relative permeability +Material(2,3) = 3.0e+70; % absolute conductivity + +% User material config coordinates in cells +% ///////////////////////////////////////////////////// +% [i, cells_i, j, cells_j] +% mat_conf = [ +% 160, 9, 240, 20; +% 180, 15, 240, 20; +% 200, 10, 240, 20; +% 220, 10, 240, 20; +% 240, 10, 240, 20; +% 260, 10, 240, 20; +% 280, 10, 240, 20; +% 300, 10, 240, 20; +% 320, 10, 240, 20; +% ]; + +mat_conf = [ + 200, 100, 200, 100; + + ]; + +mat_count = size(mat_conf, 1); +% ///////////////////////////////////////////////////// + + +% Fill geometry matrix for 1/2 cylinder (no vectorization here), but for TE +% mode we need three different Index for Ex, Ey and Hz due to leapfrog +for I = 1:nx + for J = 1:ny +% half sphere +% if ((J-nx/2)^2 + (I-ny/2)^2 <= 0.25*(d0/dx)^2 && (area_width-I*dx<=J*dy)) +% Index(J,I) = 1; +% end +% if ((J-nx/2)^2 + (I+0.5-ny/2)^2 <= 0.25*(d0/dx)^2 && (area_width-(I+0.5)*dx<=J*dy)) +% IndexX(J,I) = 1; +% end +% if ((J+0.5-nx/2)^2 + (I-ny/2)^2 <= 0.25*(d0/dx)^2 && (area_width-I*dx<=(J+0.5)*dy)) +% IndexY(J,I) = 1; +% end + + +% squares + for mat_idx = mmm:mat_count + conf = mat_conf(mat_idx); + x = mat_conf(mat_idx, 1); + x_cells = mat_conf(mat_idx, 2); + y = mat_conf(mat_idx, 3); + y_cells = mat_conf(mat_idx, 4); + + if ((I > x) && (I < (x + x_cells))... + && (J > y) && (J < (y + y_cells))) + Index(J,I) = 1; + end + if ((I+0.5 > x) && (I+0.5 < (x + x_cells))... + && (J > y) && (J < (y + y_cells))) + IndexX(J,I) = 1; + end + if ((I > x) && (I < (x + x_cells))... + && (J+0.5 > y) && (J+0.5 < (y + y_cells))) + IndexY(J,I) = 1; + end + end + end +end + + +%% Calculate size of total field area in TF/SF formalism +nx_a = round(len_tf_dx/dx); +nx_b = round((area_width-len_tf_dx)/dx); +ny_a = round(len_tf_dy/dy); +ny_b = round((area_height-len_tf_dy)/dy); + +%% Pre-allocate 1D fields for TF/SF interface +% TM components +Ez_1D = zeros(nx_b+1,1); +Fz_1D = zeros(nx_b+1,1); +Hy_1D = zeros(nx_b,1); + +k_Fz_a = zeros(nx_b+1,1); +k_Fz_b = zeros(nx_b+1,1); +k_Hy_a = zeros(nx_b,1); + +k_Hy_b = zeros(nx_b,1); +k_Ez_a = (2.0*epsilon_0*Material(1,1) - Material(1,3)*dt)/(2.0*epsilon_0*Material(1,1) + Material(1,3)*dt); +k_Ez_b = 2.0/(2.0*epsilon_0*Material(1,1) + Material(1,3)*dt); + +% TE components +Hz_1D = zeros(nx_b+1,1); +Ey_1D = zeros(nx_b,1); +k_Hz_a = zeros(nx_b+1,1); +k_Hz_b = zeros(nx_b+1,1); +k_Ey_a = zeros(nx_b,1); +k_Ey_b = zeros(nx_b,1); + +% Free space plane wave coefficients +k_Hy_a(:) = 1.0; +k_Hy_b(:) = dt/(mu_0*Material(1,2)*dx); +k_Fz_a(:) = 1.0; +k_Fz_b(:) = dt/dx; + +k_Hz_a(:) = 1.0; +k_Hz_b(:) = dt/(mu_0*Material(1,2)*dx); +k_Ey_a(:) = 1.0; +k_Ey_b(:) = dt/(epsilon_0*Material(1,1)*dx); + +ka_max = 1; +m = 4; +R_err = 1e-16; +eta = sqrt(mu_0*Material(1,2)/epsilon_0/Material(1,1)); + +sigma_max = -(m+1)*log(R_err)/(2*eta*pml_width*dx); +sigma_x = zeros(1, pml_width); +ka_x = zeros(1, pml_width); + +for ii = 0:(pml_width-1) + sigma_x(ii+mmm) = sigma_max*((pml_width - ii)/ pml_width).^m; + ka_x(ii+mmm) = 1 + (ka_max - 1)*((pml_width - ii)/ pml_width).^m; + + k_Fz_a(end-ii) = (2*epsilon_0*ka_x(ii+mmm) - sigma_x(ii+mmm)*dt)./... + (2*epsilon_0*ka_x(ii+mmm) + sigma_x(ii+mmm)*dt); + + k_Fz_b(end-ii) = 2*epsilon_0*dt./(2*epsilon_0*ka_x(ii+mmm)+sigma_x(ii+mmm)*dt)/dx; + k_Hz_a(end-ii) = (2*epsilon_0*ka_x(ii+mmm) - sigma_x(ii+mmm)*dt)./(2*epsilon_0*ka_x(ii+mmm) + sigma_x(ii+mmm)*dt); + k_Hz_b(end-ii) = 2*epsilon_0*dt./(2*epsilon_0*ka_x(ii+mmm)+sigma_x(ii+mmm)*dt)/(mu_0*Material(1,2)*dx); +end + +for ii = 0:(pml_width-1) + sigma_x(ii+mmm) = sigma_max*((pml_width - ii - 0.5)/ pml_width).^m; + ka_x(ii+mmm) = 1 + (ka_max - 1)*((pml_width - ii - 0.5)/ pml_width).^m; + k_Hy_a(end-ii) = (2*epsilon_0*ka_x(ii+mmm) - sigma_x(ii+mmm)*dt)./(2*epsilon_0*ka_x(ii+mmm) + sigma_x(ii+mmm)*dt); + k_Hy_b(end-ii) = 2*epsilon_0*dt./(2*epsilon_0*ka_x(ii+mmm) + sigma_x(ii+mmm)*dt)/(mu_0*Material(1,2)*dx); + k_Ey_a(end-ii) = (2*epsilon_0*ka_x(ii+mmm) - sigma_x(ii+mmm)*dt)./(2*epsilon_0*ka_x(ii+mmm) + sigma_x(ii+mmm)*dt); + k_Ey_b(end-ii) = 2*dt./(2*epsilon_0*ka_x(ii+mmm) + sigma_x(ii+mmm)*dt)/(Material(1,1)*dx); +end + +%% Another transformation coefficients +for ii = mmm:number_of_materials + K_a(ii) = (2*epsilon_0*Material(ii,1) - ... + Material(ii,3)*dt)./(2*epsilon_0*Material(ii,1) + ... + Material(ii,3)*dt); + K_b(ii) = 2*epsilon_0*Material(ii,1)./ ... + (2*epsilon_0*Material(ii,1) + Material(ii,3)*dt); + + K_a1(ii) = 1./(2*epsilon_0*Material(ii,1) + ... + Material(ii,3)*dt); + K_b1(ii) = 2*epsilon_0*Material(ii,1) - ... + Material(ii,3)*dt; +end + +% size(K_a) +% ????????????????????????? +% ????????????????????????? +% ????????????????????????? +% K_a_new = K_a(Index(mmm+1:nx-1,mmm+1:ny-1)+1); +K_a_new = Index(mmm+1:nx-1,mmm+1:ny-1)+1; +K_b_new = Index(mmm+1:nx-1,mmm+1:ny-1)+1; +% K_old = K_a_new; + +% K_b_new = K_b(K_b_new); +% K_a_new = K_a(K_a_new); +% Replace [1,2] -> [1,-1] +for i = mmm:nx-2 + for j = mmm:ny-2 + if(K_a_new(i,j) == 2) + K_a_new(i,j) = K_a(2); + end + if(K_b_new(i,j) == 2) + K_b_new(i,j) = K_b(2); + end + end +end + + + + + +% ssssssssssssssssssssssssssssssssssssssssssssss +% pppp = 0; +% io = 0; +% jo=0; +% for i = mmm:nx-2 +% for j = mmm:ny-2 +% if(K_a_new(i,j) ~= K_old(i,j)) +% pppp = pppp+1; +% io = i; +% jo = j; +% end +% end +% end +% pppp +% io; +% jo; +% ssssssssssssssssssssssssssssssssssssssssssssss + +% aaa = [1, 1]; +% bbb = [1,2,3,4; 5,6,7,8; 9,10,11,12; 13,14,15,16;]; +% aaa(Index(100:104-1,100:104-1)+1); + +% Index(100:104-1,100:104-1)+2 +% BB = aaa() + + +% size(K_a) +% ????????????????????????? +% K_a_new = zeros(nx-2,ny-2); +% K_b_new = zeros(nx-2,ny-2); +% for i = mmm+1:nx-1 +% for j = mmm+1:ny-1 +% K_a_new = Index(i,j)+1; +% K_b_new = Index(i,j)+1; +% end +% end + + +%% Allocate 2D arrays +% TM physical and auxiliary fields +Fz = zeros(nx,ny); Tz = zeros(nx,ny); Gx = zeros(nx,ny-1); Gy = zeros(nx-1,ny); +Ez = zeros(nx,ny); Hx = zeros(nx,ny-1); Hy = zeros(nx-1,ny); +% TE physical and auxiliary fields +Wz = zeros(nx,ny); Mx = zeros(nx,ny-1); My = zeros(nx-1,ny); +Hz = zeros(nx,ny); Ex = zeros(nx,ny-1); Ey = zeros(nx-1,ny); + +%% Allocate UPML FDTD 1D coefficient arrays +% TM coefficients +k_Fz_1 = zeros(ny,1); k_Fz_2 = zeros(ny,1); +k_Ez_1 = zeros(nx,1); k_Ez_2 = zeros(nx,1); +k_Gx_1 = zeros(ny-1,1); k_Gx_2 = zeros(ny-1,1); +k_Hx_1 = zeros(nx,1); k_Hx_2 = zeros(nx,1); +k_Gy_1 = zeros(nx-1,1); k_Gy_2 = zeros(nx-1,1); +k_Hy_1 = zeros(ny,1); k_Hy_2 = zeros(ny,1); +% TM coefficients +k_Hz_1 = zeros(nx,1); k_Hz_2 = zeros(nx,1); +k_Ex_1 = zeros(nx,1); k_Ex_2 = zeros(nx,1); +k_Ey_1 = zeros(ny,1); k_Ey_2 = zeros(ny,1); + +%% Free space FDTD coefficients +k_Fz_1(:) = 1.0; k_Fz_2(:) = dt; +k_Ez_1(:) = 1.0; k_Ez_2(:) = 1.0/epsilon_0; +k_Gx_1(:) = 1.0; k_Gx_2(:) = dt/dy; +k_Hx_1(:) = 1.0/mu_0; k_Hx_2(:) = 1.0/mu_0; +k_Gy_1(:) = 1.0; k_Gy_2(:) = dt/dx; +k_Hy_1(:) = 1.0/mu_0; k_Hy_2(:) = 1.0/mu_0; +k_Hz_1(:) = 1.0; k_Hz_2(:) = 1.0/mu_0; +k_Ex_1(:) = 2.0; k_Ex_2(:) = 2.0; +k_Ey_1(:) = 2.0; k_Ey_2(:) = 2.0; + +% //////////////////// LABEL //////////////////////// + +%% Field transformation coefficients in pml_width areas for TM ans TE modes +% Along x-axis +sigma_max = -(m+1)*log(R_err)/(2*eta*pml_width*dx); + +for i = mmm:pml_width + sigma_x(i) = sigma_max*((pml_width-i+1)/pml_width).^m; + ka_x(i) = 1+(ka_max-1)*((pml_width-i+1)/pml_width).^m; + + k_Ez_1(i) = (2*epsilon_0*ka_x(i)-sigma_x(i)*dt)./(2*epsilon_0*ka_x(i)+sigma_x(i)*dt); + k_Ez_1(end-(i)+1) = k_Ez_1(i); + k_Ez_2(i) = 2./(2*epsilon_0*ka_x(i) + sigma_x(i)*dt); + k_Ez_2(end-(i)+1) = k_Ez_2(i); + k_Hx_1(i) = (2*epsilon_0*ka_x(i)+sigma_x(i)*dt)/(2*epsilon_0*mu_0); + k_Hx_1(end-(i)+1) = k_Hx_1(i); + k_Hx_2(i) = (2*epsilon_0*ka_x(i)-sigma_x(i)*dt)/(2*epsilon_0*mu_0); + k_Hx_2(end-(i)+1) = k_Hx_2(i); + k_Hz_1(i) = (2*epsilon_0*ka_x(i)-sigma_x(i)*dt)./(2*epsilon_0*ka_x(i)+sigma_x(i)*dt); + k_Hz_1(end-(i)+1) = k_Hz_1(i); + k_Hz_2(i) = 2*epsilon_0./(2*epsilon_0*ka_x(i)+sigma_x(i)*dt)/mu_0; + k_Hz_2(end-(i)+1) = k_Hz_2(i); + k_Ex_1(i) = (2*epsilon_0*ka_x(i)+sigma_x(i)*dt)/epsilon_0; + k_Ex_1(end-(i)+1) = k_Ex_1(i); + k_Ex_2(i) = (2*epsilon_0*ka_x(i)-sigma_x(i)*dt)/epsilon_0; + k_Ex_2(end-(i)+1) = k_Ex_2(i); +end + +for i = mmm:pml_width + sigma_x(i) = sigma_max*((pml_width-(i)+0.5)/pml_width).^m; + ka_x(i) = 1+(ka_max-1)*((pml_width-(i)+0.5)/pml_width).^m; + k_Gy_1(i) = (2*epsilon_0*ka_x(i)-sigma_x(i)*dt)./(2*epsilon_0*ka_x(i)+sigma_x(i)*dt); + k_Gy_1(end-(i)+1) = k_Gy_1(i); + k_Gy_2(i) = 2*epsilon_0*dt./(2*epsilon_0*ka_x(i)+sigma_x(i)*dt)/dx; + k_Gy_2(end-(i)+1) = k_Gy_2(i); +end + + + +% Along y-axis +sigma_max = -(m+1)*log(R_err)/(2*eta*pml_width*dy); +for i = mmm:pml_width + sigma_y(i) = sigma_max*((pml_width-(i)+1)/pml_width).^m; + ka_y(i) = 1+(ka_max-1)*((pml_width-(i)+1)/pml_width).^m; + k_Fz_1(i) = (2*epsilon_0*ka_y(i)-sigma_y(i)*dt)./(2*epsilon_0*ka_y(i)+sigma_y(i)*dt); + k_Fz_1(end-(i)+1) = k_Fz_1(i); + k_Fz_2(i) = 2*epsilon_0*dt./(2*epsilon_0*ka_y(i)+sigma_y(i)*dt); + k_Fz_2(end-(i)+1) = k_Fz_2(i); + k_Hy_1(i) = (2*epsilon_0*ka_y(i)+sigma_y(i)*dt)/(2*epsilon_0*mu_0); + k_Hy_1(end-(i)+1) = k_Hy_1(i); + k_Hy_2(i) = (2*epsilon_0*ka_y(i)-sigma_y(i)*dt)/(2*epsilon_0*mu_0); + k_Hy_2(end-(i)+1) = k_Hy_2(i); + k_Ey_1(i) = (2*epsilon_0*ka_y(i)+sigma_y(i)*dt)/epsilon_0; + k_Ey_1(end-(i)+1) = k_Ey_1(i); + k_Ey_2(i) = (2*epsilon_0*ka_y(i)-sigma_y(i)*dt)/epsilon_0; + k_Ey_2(end-(i)+1) = k_Ey_2(i); +end + +for i = mmm:pml_width + sigma_y(i) = sigma_max*((pml_width-(i)+0.5)/pml_width).^m; + ka_y(i) = 1+(ka_max-1)*((pml_width-(i)+0.5)/pml_width).^m; + k_Gx_1(i) = (2*epsilon_0*ka_y(i)-sigma_y(i)*dt)./(2*epsilon_0*ka_y(i)+sigma_y(i)*dt); + k_Gx_1(end-(i)+1) = k_Gx_1(i); + k_Gx_2(i) = 2*epsilon_0*dt./(2*epsilon_0*ka_y(i)+sigma_y(i)*dt)/dy; + k_Gx_2(end-(i)+1) = k_Gx_2(i); +end + + +% Vectorize transformation coefficients +k_Fz_1_old = k_Fz_1; +k_Fz_1_new = zeros(nx-2,ny-2); +k_Fz_1_new(:) = k_Fz_1(pml_width+1); + +k_Fz_2_old = k_Fz_2; +k_Fz_2_new = zeros(nx-2,ny-2); +k_Fz_2_new(:) = k_Fz_2(pml_width+1); + + +% left right border +for i = mmm:nx-2 + for j = mmm:pml_width-1 + % left + k_Fz_1_new(i,j) = k_Fz_1(j+1); + k_Fz_2_new(i,j) = k_Fz_2(j+1); + + % right + % ?????? index + k_Fz_1_new(i,end-j+1) = k_Fz_1(j+1); + k_Fz_2_new(i,end-j+1) = k_Fz_2(j+1); + end +end + + +k_Ez_1_old = k_Ez_1; +k_Ez_1_new = zeros(nx-2,ny-2); +k_Ez_1_new(:) = k_Ez_1(pml_width+1); + +k_Ez_2_old = k_Ez_2; +k_Ez_2_new = zeros(nx-2,ny-2); +k_Ez_2_new(:) = k_Ez_2(pml_width+1); + + +% up down border +for i = mmm:pml_width-1 + for j = mmm:ny-2 + % up + k_Ez_1_new(i,j) = k_Ez_1(i+1); + k_Ez_2_new(i,j) = k_Ez_2(i+1); + + % down + % ?????? index + k_Ez_1_new(end-i+1,j) = k_Ez_1(i+1); + k_Ez_2_new(end-i+1,j) = k_Ez_2(i+1); + end +end + + +k_Gx_1_old = k_Gx_1; +k_Gx_1_new = zeros(nx,ny-1); +k_Gx_1_new(:) = k_Gx_1(pml_width+1); + +k_Gx_2_old = k_Gx_2; +k_Gx_2_new = zeros(nx,ny-1); +k_Gx_2_new(:) = k_Gx_2(pml_width+1); + + +% left right border +for i = mmm:nx + for j = mmm:pml_width + % left + k_Gx_1_new(i,j) = k_Gx_1(j); + k_Gx_2_new(i,j) = k_Gx_2(j); + + % right + % ?????? index + k_Gx_1_new(i,end-j+1) = k_Gx_1(j); + k_Gx_2_new(i,end-j+1) = k_Gx_2(j); + end +end + + +k_Hx_1_old = k_Hx_1; +k_Hx_1_new = zeros(nx,ny-1); +k_Hx_1_new(:) = k_Hx_1(pml_width+1); + +k_Hx_2_old = k_Hx_2; +k_Hx_2_new = zeros(nx,ny-1); +k_Hx_2_new(:) = k_Hx_2(pml_width+1); + + +k_Ex_1_old = k_Ex_1; +k_Ex_1_new = zeros(nx,ny-1); +k_Ex_1_new(:) = k_Ex_1(pml_width+1); + +k_Ex_2_old = k_Ex_2; +k_Ex_2_new = zeros(nx,ny-1); +k_Ex_2_new(:) = k_Ex_2(pml_width+1); + +% up down border +for i = mmm:pml_width + for j = mmm:ny-1 + % up + k_Hx_1_new(i,j) = k_Hx_1(i); + k_Hx_2_new(i,j) = k_Hx_2(i); + + k_Ex_1_new(i,j) = k_Ex_1(i); + k_Ex_2_new(i,j) = k_Ex_2(i); + + % down + % ?????? index + k_Hx_1_new(end-i+1,j) = k_Hx_1(i); + k_Hx_2_new(end-i+1,j) = k_Hx_2(i); + + k_Ex_1_new(end-i+1,j) = k_Ex_1(i); + k_Ex_2_new(end-i+1,j) = k_Ex_2(i); + end +end + +k_Gy_1_old = k_Gy_1; +k_Gy_1_new = zeros(nx-1,ny); +k_Gy_1_new(:) = k_Gy_1(pml_width+1); + +k_Gy_2_old = k_Gy_2; +k_Gy_2_new = zeros(nx-1,ny); +k_Gy_2_new(:) = k_Gy_2(pml_width+1); + + +% up down border +for i = mmm:pml_width + for j = mmm:ny + % up + k_Gy_1_new(i,j) = k_Gy_1(i); + k_Gy_2_new(i,j) = k_Gy_2(i); + + % down + % ?????? index + k_Gy_1_new(end-i+1,j) = k_Gy_1(i); + k_Gy_2_new(end-i+1,j) = k_Gy_2(i); + end +end + + +k_Hy_1_old = k_Hy_1; +k_Hy_1_new = zeros(nx-1,ny); +k_Hy_1_new(:) = k_Hy_1(pml_width+1); + +k_Hy_2_old = k_Hy_2; +k_Hy_2_new = zeros(nx-1,ny); +k_Hy_2_new(:) = k_Hy_2(pml_width+1); + + +k_Ey_1_old = k_Ey_1; +k_Ey_1_new = zeros(nx-1,ny); +k_Ey_1_new(:) = k_Ey_1(pml_width+1); + +k_Ey_2_old = k_Ey_2; +k_Ey_2_new = zeros(nx-1,ny); +k_Ey_2_new(:) = k_Ey_2(pml_width+1); + + +% left right border +for i = mmm:nx-1 + for j = mmm:pml_width + % left + k_Hy_1_new(i,j) = k_Hy_1(j); + k_Hy_2_new(i,j) = k_Hy_2(j); + + k_Ey_1_new(i,j) = k_Ey_1(j); + k_Ey_2_new(i,j) = k_Ey_2(j); + + % right + % ?????? index + k_Hy_1_new(i,end-j+1) = k_Hy_1(j); + k_Hy_2_new(i,end-j+1) = k_Hy_2(j); + + k_Ey_1_new(i,end-j+1) = k_Ey_1(j); + k_Ey_2_new(i,end-j+1) = k_Ey_2(j); + end +end + +k_Hz_1_old = k_Hz_1; +k_Hz_1_new = zeros(nx-2,ny-2); +k_Hz_1_new(:) = k_Hz_1(pml_width+1); + +k_Hz_2_old = k_Hz_2; +k_Hz_2_new = zeros(nx-2,ny-2); +k_Hz_2_new(:) = k_Hz_2(pml_width+1); + + +% up down border +for i = mmm:pml_width-1 + for j = mmm:ny-2 + % up + k_Hz_1_new(i,j) = k_Hz_1(i+1); + k_Hz_2_new(i,j) = k_Hz_2(i+1); + + % down + % ?????? index + k_Hz_1_new(end-i+1,j) = k_Hz_1(i+1); + k_Hz_2_new(end-i+1,j) = k_Hz_2(i+1); + end +end + + +% k_Fz_1 = repmat(k_Fz_1(2:ny-1)',nx-2,1); +% k_Fz_2 = repmat(k_Fz_2(2:ny-1)',nx-2,1); + +% k_Ez_1 = repmat(k_Ez_1(2:nx-1),1,ny-2); +% k_Ez_2 = repmat(k_Ez_2(2:nx-1),1,ny-2); + +% k_Gx_1 = repmat(k_Gx_1(1:ny-1)',nx,1); +% k_Gx_2 = repmat(k_Gx_2_new(1:ny-1)',nx,1); + +% k_Hx_1 = repmat(k_Hx_1(1:nx),1,ny-1); +% k_Hx_2 = repmat(k_Hx_2(1:nx),1,ny-1); + +% k_Gy_1 = repmat(k_Gy_1(1:nx-1),1,ny); +% k_Gy_2 = repmat(k_Gy_2(1:nx-1),1,ny); + +% k_Hy_1 = repmat(k_Hy_1(1:ny)',nx-1,1); +% k_Hy_2 = repmat(k_Hy_2(1:ny)',nx-1,1); + +% k_Hz_1 = repmat(k_Hz_1(2:nx-1),1,ny-2); +% k_Hz_2 = repmat(k_Hz_2(2:nx-1),1,ny-2); + +% k_Ex_1 = repmat(k_Ex_1(1:nx),1,ny-1); +% k_Ex_2 = repmat(k_Ex_2(1:nx),1,ny-1); + + +% k_Ey_1 = repmat(k_Ey_1(1:ny)',nx-1,1); +% k_Ey_2 = repmat(k_Ey_2(1:ny)',nx-1,1); + + +% M0_old = Material(Index(2:nx-1,2:ny-1)+1,1); + +% M0 = reshape(Material(Index(2:nx-1,2:ny-1)+1,1),nx-2,[]); +% M1 = reshape(Material(Index(2:nx-1,2:ny-1)+1,2),nx-2,[]); +% M2 = reshape(Material(Index(1:nx,1:ny-1)+1,2),nx,[]); +% M3 = reshape(Material(Index(1:nx-1,1:ny)+1,2),[],ny); + + +M0 = zeros(nx-2,ny-2); +M1 = zeros(nx-2,ny-2); + +for i = mmm:nx-2 + for j = mmm:ny-2 + M0(i,j) = Material(Index(i+1,j+1)+1, 1); + M1(i,j) = Material(Index(i+1,j+1)+1, 2); + end +end + + +M2 = zeros(nx,ny-1); + +for i = mmm:nx + for j = mmm:ny-1 + M2(i,j) = Material(Index(i,j)+1, 2); + end +end + +M3 = zeros(nx-1,ny); + +for i = mmm:nx-1 + for j = mmm:ny + M3(i,j) = Material(Index(i,j)+1, 2); + end +end + +% M0 = reshape(Material(Index(2:nx-1,2:ny-1)+1,1),nx-2,ny-2); +% M1 = reshape(Material(Index(2:nx-1,2:ny-1)+1,2),nx-2,ny-2); +% M2 = reshape(Material(Index(1:nx,1:ny-1)+1,2),nx,ny-1); +% M3 = reshape(Material(Index(1:nx-1,1:ny)+1,2),nx-1,ny); + + + + +%% Create fullscreen figure and set a double-buffer +figure('units','normalized','outerposition',[0 0 1 1]); +set(gcf,'doublebuffer','on'); + +tic; + +Fz_1D_r = zeros(nx_b-1, 1); +Fz_r1 = zeros(nx-2, ny-2); +Fz_r = zeros(nx-2, ny-2); +Tz_r = zeros(nx-2, ny-2); + +Gy_r = zeros(nx-1, ny); +Gx_r = zeros(nx, ny-1); +Gx_r1 = zeros(nx, ny-1); +Gy_r1 = zeros(nx-1, ny); + + +%% Main FDTD UPML routine with TF/SF excitation interface calculates TE and +%% TM modes simultaneously +for T = 1:number_of_iterations + %% Calculate incident 1D plain waves for TF/SF implementation + % TE mode (Hz) + for i = mmm:nx_b + Ey_1D(i) = k_Ey_a(i)*Ey_1D(i) ... + - k_Ey_b(i)*(Hz_1D(i+1)-Hz_1D(i)); + end + + + Hz_1D(1) = E0*sin(2*pi*f*(T-1)*dt); + + for i = mmm+1:nx_b + Hz_1D(i) = k_Hz_a(i)*Hz_1D(i) ... + - k_Hz_b(i)*(Ey_1D(i)-Ey_1D(i-1)); + end + + % TM mode (Ez) + for i = mmm:nx_b + Hy_1D(i) = k_Hy_a(i)*Hy_1D(i) + ... + k_Hy_b(i)*(Ez_1D(i+1)-Ez_1D(i)); + end + + Ez_1D(1) = E0*sin(2*pi*f*(T-1)*dt); + + % Fz_1D_r = Fz_1D(2:nx_b); + for i = mmm:nx_b-1 + Fz_1D_r(i) = Fz_1D(i+1); + end + + for i = mmm+1:nx_b + Fz_1D(i) = k_Fz_a(i)*Fz_1D(i) + ... + k_Fz_b(i)*( Hy_1D(i) - Hy_1D(i-1) ); + end + + + for i = mmm+1:nx_b + Ez_1D(i) = k_Ez_a*Ez_1D(i) + k_Ez_b*( Fz_1D(i) - Fz_1D_r(i-1)); + end + + %% Calculate Ez (TM mode) and Hz (TE mode) total fields + % TE: Wz -> Hz + + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Fz_r1(i-1,j-1) = Wz(i,j); + end + end + + + + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Wz(i,j) = k_Fz_1_new(i-1,j-1) * Wz(i,j) + ... + k_Fz_2_new(i-1,j-1) * ((Ex(i,j)-Ex(i,j-1))/dy - ... + (Ey(i,j)-Ey(i-1,j))/dx ); + + Hz(i,j) = k_Hz_1_new(i-1,j-1)*Hz(i,j) + ... + k_Hz_2_new(i-1,j-1) * ( Wz(i,j)-Fz_r1(i-1,j-1) )/M1(i-1,j-1); + end + end + % Wz(2:nx-1,2:ny-1) = k_Fz_1_new.*Wz(2:nx-1,2:ny-1) + ... + % k_Fz_2_new.*( (Ex(2:nx-1,2:ny-1)-Ex(2:nx-1,1:ny-2))/dy - ... + % (Ey(2:nx-1,2:ny-1)-Ey(1:nx-2,2:ny-1))/dx ); + % Hz(2:nx-1,2:ny-1) = k_Hz_1_new.*Hz(2:nx-1,2:ny-1) + ... + % k_Hz_2_new.*( Wz(2:nx-1,2:ny-1)-Fz_r1 )./M1; + + + % TM: Fz -> Tz -> Ez + % //////////////////////////////////////////////// label 2 + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Fz_r(i-1,j-1) = Fz(i,j); + Tz_r(i-1,j-1) = Tz(i,j); + end + end + + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Fz(i,j) = k_Fz_1_new(i-1,j-1) * Fz(i,j) + k_Fz_2_new(i-1,j-1)*( (Hy(i,j) - ... + Hy(i-1,j))/dx - (Hx(i,j) - ... + Hx(i,j-1))/dy ); + end + end + + + for i = mmm:nx-2 + for j = mmm:ny-2 + Tz(i+1,j+1) = K_a_new(i,j)* ... + Tz(i+1,j+1) + ... + K_b_new(i,j)*( Fz(i+1,j+1) - Fz_r(i,j)); + end + end + + + for i = mmm+1:nx-1 + for j = mmm+1:ny-1 + Ez(i,j) = k_Ez_1_new(i-1,j-1)*Ez(i,j) + ... + k_Ez_2_new(i-1,j-1)*( Tz(i,j) - Tz_r(i-1,j-1) )/ M0(i-1,j-1); + end + end + + %% Calculate scattered field Ez and Hz in TF/SF + % TE mode + i = nx_a+1; + for j = ny_a+1:ny_b+1 + Hz(i,j) = Hz(i,j) + ... + dt / (mu_0*Material(Index(i,j)+1, 2)*dx) * Ey_1D(1); + end + + i = nx_b+1; + for j = ny_a+1:ny_b+1 + Hz(i,j) = Hz(i,j) - ... + dt/(mu_0*Material(Index(i,j)+1, 2)*dx)*Ey_1D(nx_b-nx_a+2); + end + + % TM mode + i = nx_a+1; + for j = ny_a+1:ny_b+1 + Ez(i,j) = Ez(i,j) - ... + dt./(epsilon_0*Material(Index(i,j)+1,1)*dx)*Hy_1D(1); + end + + i = nx_b+1; + for j = ny_a+1:ny_b+1 + Ez(i,j) = Ez(i,j) + ... + dt/(epsilon_0*Material(Index(i,j)+1,1)*dx)*Hy_1D(nx_b-nx_a+2); + end + + %% Calculate Hx and Ex total fields + % TE mode + for i = mmm:nx + for j = mmm:ny-1 + Gx_r1(i,j) = Mx(i,j); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Mx(i,j) = k_Gx_1_new(i,j) * Mx(i,j) + ... + k_Gx_2_new(i,j)*( Hz(i,j+1)-Hz(i,j) ); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Ex(i,j) = K_a1(IndexX(i,j)+1)*( K_b1(IndexX(i,j)+1)*... + Ex(i,j) + k_Ex_1_new(i,j)*Mx(i,j)-k_Ex_2_new(i,j)*Gx_r1(i,j) ); + end + end + + % TM mode + for i = mmm:nx + for j = mmm:ny-1 + Gx_r(i,j) = Gx(i,j); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Gx(i,j) = k_Gx_1_new(i,j) * Gx(i,j) - ... + k_Gx_2_new(i,j) * ( Ez(i,j+1)-Ez(i,j) ); + end + end + + for i = mmm:nx + for j = mmm:ny-1 + Hx(i,j) = Hx(i,j) + (k_Hx_1_new(i,j) * Gx(i,j) - k_Hx_2_new(i,j)*Gx_r(i,j)) / M2(i,j); + end + end + + %% Calculate Hy and Ey total fields + % TE mode + for i = mmm:nx-1 + for j = mmm:ny + Gy_r1(i,j) = My(i,j); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + My(i,j) = k_Gy_1_new(i,j)*My(i,j) - ... + k_Gy_2_new(i,j)*( Hz(i+1,j)-Hz(i,j)); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + Ey(i,j) = K_a1(IndexY(i,j)+1)*( K_b1(IndexY(i,j)+1)*... + Ey(i,j) + k_Ey_1_new(i,j)*My(i,j)-k_Ey_2_new(i,j)*Gy_r1(i,j) ); + end + end + + % TM mode + for i = mmm:nx-1 + for j = mmm:ny + Gy_r(i,j) = Gy(i,j); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + Gy(i,j) = k_Gy_1_new(i,j)*Gy(i,j) + ... + k_Gy_2_new(i,j)*(Ez(i+1,j)-Ez(i,j)); + end + end + + for i = mmm:nx-1 + for j = mmm:ny + Hy(i,j) = Hy(i,j) + ... + (k_Hy_1_new(i,j)*Gy(i,j) - k_Hy_2_new(i,j)*Gy_r(i,j))/M3(i,j); + end + end + + %% Calculate scattered field Hx and Ex in TF/SF + % TE mode + j = ny_a; + for i = nx_a+1:nx_b+1 + Ex(i,j) = Ex(i,j) - ... + 2*dt/dy*K_a1(Index(i,j)+1)*Hz_1D(i-nx_a+1); + end + + j = ny_b+1; + for i = nx_a+1:nx_b+1 + Ex(i,j) = Ex(i,j) + ... + 2*dt/dy*K_a1(Index(i,j)+1)*Hz_1D(i-nx_a+1); + end + + % TM mode + j = ny_a; + for i = nx_a+1:nx_b+1 + Hx(i,j) = Hx(i,j) + ... + dt/(mu_0*dy*Material(Index(i,j)+1,2))*Ez_1D(i-nx_a+1); + end + + j = ny_b+1; + for i = nx_a+1:nx_b+1 + Hx(i,j) = Hx(i,j) - ... + dt/(mu_0*dy*Material(Index(i,j)+1,2))*Ez_1D(i-nx_a+1); + end + + %% Calculate scattered field Hy and Ey in TF/SF + % TE mode + i = nx_a; + for j = ny_a+1:ny_b+1 + Ey(i,j) = Ey(i,j) + ... + 2*dt/dx*K_a1(Index(i,j)+1,1)*Hz_1D(mmm+1); + end + + i = nx_b+1; + for j = ny_a+1:ny_b+1 + Ey(i,j) = Ey(i,j) - ... + 2*dt/dx*K_a1(Index(i,j)+1,1)*Hz_1D(nx_b-nx_a+2); + end + + % TM mode + i = nx_a; + for j = ny_a+1:ny_b+1 + Hy(i,j) = Hy(i,j) - ... + dt/(mu_0*dx*Material(Index(i,j)+1,2))*Ez_1D(2); + end + + i = nx_b+1; + for j = j + Hy(i,j) = Hy(i,j) + ... + dt/(mu_0*dx*Material(Index(i,j)+1,2))*Ez_1D(nx_b-nx_a+2); + end + + %% Plot Ez and Hz fields dynamics + if (mod(T,8) == 0) + subplot(1,2,1); + pcolor(Y(pml_width:ny-pml_width),X(pml_width:nx-pml_width),Ez(pml_width:nx-pml_width,pml_width:ny-pml_width)); + shading interp; + caxis([-E0 E0]); + axis image; + colorbar; + title('E_{z}(x,y)', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + xlabel('x, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + ylabel('y, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + subplot(1,2,2); + pcolor(Y(pml_width:ny-pml_width),X(pml_width:nx-pml_width),Hz(pml_width:nx-pml_width,pml_width:ny-pml_width)); + shading interp; + caxis([-E0 E0]); + axis image; + colorbar; + title('H_{z}(x,y)', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + xlabel('x, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + ylabel('y, [m]', 'FontSize',18, 'FontName', 'Arial', 'FontWeight', 'Bold'); + drawnow; + end +end +disp(['Time elapsed - ',num2str(toc/60),' minutes']);