diff --git a/PlotScript/Animations/animation_Ex.gif b/PlotScript/Animations/animation_Ex.gif index 03d3fcd..2e710db 100644 Binary files a/PlotScript/Animations/animation_Ex.gif and b/PlotScript/Animations/animation_Ex.gif differ diff --git a/include/FDTD.h b/include/FDTD.h index ead71e2..429af27 100644 --- a/include/FDTD.h +++ b/include/FDTD.h @@ -6,27 +6,11 @@ #include #include -#include "Structures.h" +#include "Writer.h" +#include "Field.h" using namespace FDTDstruct; -class Field -{ -private: - int Ni, Nj, Nk; - std::vector field; - -public: - Field(const int _Ni, const int _Nj, const int _Nk); - Field& operator= (const Field& other); - - double& operator() (int i, int j, int k); - - int get_Ni() { return Ni; } - int get_Nj() { return Nj; } - int get_Nk() { return Nk; } -}; - class FDTD { private: @@ -59,16 +43,12 @@ class FDTD void set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], double SGm, std::function dist); - double calc_max_value(Field& field); - public: FDTD(Parameters _parameters, double _dt, double _pml_percent); Field& get_field(Component); std::vector& get_current(Component); - std::vector> update_fields(const int time); - - double calc_reflection(Field E_start[3], Field B_start[3], - Field E_final[3], Field B_final[3]); + std::vector update_fields(const int time, bool write_result = false, + Axis write_axis = Axis::X, std::string base_path = ""); }; diff --git a/include/Field.h b/include/Field.h new file mode 100644 index 0000000..f79fdf4 --- /dev/null +++ b/include/Field.h @@ -0,0 +1,20 @@ +#pragma once + +#include + +class Field +{ +private: + int Ni, Nj, Nk; + std::vector field; + +public: + Field(const int _Ni = 1, const int _Nj = 1, const int _Nk = 1); + Field& operator= (const Field& other); + + double& operator() (int i, int j, int k); + + int get_Ni() { return Ni; } + int get_Nj() { return Nj; } + int get_Nk() { return Nk; } +}; diff --git a/include/Writer.h b/include/Writer.h index 9ef7f64..dae03fb 100644 --- a/include/Writer.h +++ b/include/Writer.h @@ -1,172 +1,16 @@ #pragma once +#include "Structures.h" +#include "Field.h" + +#include #include #include -#include "FDTD.h" - -void write_x(Field& this_field, std::ofstream& fout) -{ - for (int k = 0; k < this_field.get_Nk(); ++k) - { - for (int j = 0; j < this_field.get_Nj(); ++j) - { - for (int i = 0; i < this_field.get_Ni(); ++i) - { - fout << this_field(i, j, k); - if (i == this_field.get_Ni() - 1) - { - fout << std::endl; - } - else - { - fout << ";"; - } - } - } - } - fout << std::endl << std::endl; -} -void write_y(Field& this_field, std::ofstream& fout) -{ - for (int k = 0; k < this_field.get_Nk(); ++k) - { - for (int i = 0; i < this_field.get_Ni(); ++i) - { - for (int j = 0; j < this_field.get_Nj(); ++j) - { - fout << this_field(i, j, k); - if (j == this_field.get_Nj() - 1) - { - fout << std::endl; - } - else - { - fout << ";"; - } - } - } - } - fout << std::endl << std::endl; -} -void write_z(Field& this_field, std::ofstream& fout) -{ - for (int i = 0; i < this_field.get_Ni(); ++i) - { - for (int j = 0; j < this_field.get_Nj(); ++j) - { - for (int k = 0; k < this_field.get_Nk(); ++k) - { - fout << this_field(i, j, k); - if (k == this_field.get_Nk() - 1) - { - fout << std::endl; - } - else - { - fout << ";"; - } - } - } - } - fout << std::endl << std::endl; -} - -void write_spherical(std::vector> data, Axis axis, int iteration, std::string base_path) -{ - std::ofstream test_fout; - - int it = 0; - for (std::vector fields : data) - { - it++; - for (int c = static_cast(Component::EX); c <= static_cast(Component::BZ); ++c) - { - std::cout << "Write : " << it << " -- " << c << std::endl; - if (axis == Axis::X) - { - test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv"); - Field field = fields[c]; - for (int k = 0; k < field.get_Nk(); ++k) - { - for (int j = 0; j < field.get_Nj(); ++j) - { - test_fout << field(field.get_Ni() / 2, j, k); - if (j == field.get_Nj() - 1) - { - test_fout << std::endl; - } - else - { - test_fout << ";"; - } - } - } - test_fout.close(); - } - else if (axis == Axis::Y) - { - test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv"); - Field field = fields[c]; - for (int i = 0; i < field.get_Nj(); ++i) - { - for (int k = 0; k < field.get_Nk(); ++k) - { - test_fout << field(i, field.get_Nj() / 2, k); - if (k == field.get_Nk() - 1) - { - test_fout << std::endl; - } - else - { - test_fout << ";"; - } - } - } - test_fout.close(); - } - else - { - test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv"); - Field field = fields[c]; - for (int i = 0; i < field.get_Nj(); ++i) - { - for (int j = 0; j < field.get_Nk(); ++j) - { - test_fout << field(i, j, field.get_Nk() / 2); - if (j == field.get_Nk() - 1) - { - test_fout << std::endl; - } - else - { - test_fout << ";"; - } - } - } - test_fout.close(); - } - } - } -} +using namespace FDTDstruct; -void write_plane(FDTD& test, Axis axis, char* file_path) -{ - std::ofstream test_fout; - test_fout.open(file_path); +void write_x(Field& this_field, std::ofstream& fout); +void write_y(Field& this_field, std::ofstream& fout); +void write_z(Field& this_field, std::ofstream& fout); - if (!test_fout.is_open()) - { - throw std::runtime_error("ERROR: Failed to open " + static_cast(file_path)); - } - for (int i = static_cast(Component::EX); i <= static_cast(Component::BZ); ++i) - { - if (axis == Axis::X) - write_x(test.get_field(static_cast(i)), test_fout); - else if (axis == Axis::Y) - write_y(test.get_field(static_cast(i)), test_fout); - else - write_z(test.get_field(static_cast(i)), test_fout); - } - test_fout.close(); -} +void write_spherical(std::vector fields, Axis axis, std::string base_path, int it); diff --git a/samples/sample.cpp b/samples/sample.cpp index 97c5a93..62c86d9 100644 --- a/samples/sample.cpp +++ b/samples/sample.cpp @@ -1,5 +1,4 @@ #define _USE_MATH_DEFINES -#define __USE_SPHERICAL_WAVE__ #include #include @@ -8,7 +7,6 @@ #include #include "test_FDTD.h" -#include "Writer.h" Axis get_axis(Component field_E, Component field_B) { @@ -36,86 +34,6 @@ Axis get_axis(Component field_E, Component field_B) return selected_axis; } -void plane_wave(std::vector arguments, - std::vector numbers, char* outfile_path) -{ - // Saving data from the console - Component selected_field_1 = static_cast(std::atoi(arguments[1])); - Component selected_field_2 = static_cast(std::atoi(arguments[2])); - Component fld_tested = static_cast(std::atoi(arguments[3])); - bool write_flag = static_cast(std::atoi(arguments[4])); - - Component flds_selected[2] = { selected_field_1, selected_field_2 }; - - // Saving data from the txt file - int grid_sizes[3] = { numbers[0], numbers[1], numbers[2] }; // { Nx, Ny , Nz } - double sizes_x[2] = { numbers[3], numbers[4] }; // { ax, bx } - double sizes_y[2] = { numbers[5], numbers[6] }; // { ay, by } - double sizes_z[2] = { numbers[7], numbers[8] }; // { az, bz } - int iterations_num = static_cast(numbers[9]); - double time = numbers[10]; - double time_step = time / static_cast(iterations_num); // dt - - // Initialization of the initializing function and the true solution function - std::function initial_func = [](double x, double size[2]) - { - return sin(2.0 * M_PI * (x - size[0]) / (size[1] - size[0])); - }; - - // Determination of the wave propagation axis - Axis axis = get_axis(flds_selected[0], flds_selected[1]); - - // Initialization of the structures and method - SelectedFields current_fields - { - flds_selected[0], - flds_selected[1], - }; - - Parameters params - { - grid_sizes[0], - grid_sizes[1], - grid_sizes[2], - sizes_x[0], - sizes_x[1], - sizes_y[0], - sizes_y[1], - sizes_z[0], - sizes_z[1], - (sizes_x[1] - sizes_x[0]) / static_cast(grid_sizes[0]), - (sizes_y[1] - sizes_y[0]) / static_cast(grid_sizes[1]), - (sizes_z[1] - sizes_z[0]) / static_cast(grid_sizes[2]) - }; - FDTD method(params, time_step, 0.0); - - // Meaningful calculations - Test_FDTD test(params); - try - { - test.initial_filling(method, current_fields, iterations_num, initial_func); - } - catch(const std::exception& except) - { - std::cout << except.what() << std::endl; - } - - std::vector> data = method.update_fields(iterations_num); - - if (write_flag == true) - { - // Writing the results to a file - try - { - write_plane(method, axis, outfile_path); - } - catch (const std::exception& except) - { - std::cout << except.what() << std::endl; - } - } -} - void spherical_wave(int n, int it, char* base_path = "") { CurrentParameters cur_param @@ -163,15 +81,7 @@ void spherical_wave(int n, int it, char* base_path = "") // Meaningful calculations Test_FDTD test(params); test.initiialize_current(method, cur_param, it, cur_func); - std::vector> data = method.update_fields(it); - try - { - write_spherical(data, Axis::X, it, base_path); - } - catch (const std::exception& except) - { - std::cout << except.what() << std::endl; - } + method.update_fields(it, true, Axis::X, ""); } int main(int argc, char* argv[]) @@ -182,46 +92,9 @@ int main(int argc, char* argv[]) std::vector arguments(argv, argv + argc); if (argc == 1) { -#ifdef __USE_SPHERICAL_WAVE__ int N = 70; int Iterations = 410; spherical_wave(N, Iterations, "../../PlotScript/"); -#endif - -#ifndef __USE_SPHERICAL_WAVE__ - source_fin.open("../../PlotScript/Source.txt"); - outfile_path = "../../PlotScript/OutFile.csv"; - - struct StringComponent { - char* EX = "0"; - char* EY = "1"; - char* EZ = "2"; - char* BX = "3"; - char* BY = "4"; - char* BZ = "5"; - } field; - const size_t size_tmp = 5; - char* tmp[size_tmp]{ "1", field.EX, field.BZ, field.BZ, "1" }; - - arguments.clear(); - arguments.reserve(size_tmp); - std::copy(std::begin(tmp), std::end(tmp), std::back_inserter(arguments)); - - if (!source_fin.is_open()) - { - std::cout << "ERROR: Failed to open Source.txt" << std::endl; - return 1; - } - std::vector numbers; - double number; - while (source_fin >> number) - { - numbers.push_back(number); - } - source_fin.close(); - - plane_wave(arguments, numbers, outfile_path); -#endif } else if (argc == 4) { @@ -231,21 +104,8 @@ int main(int argc, char* argv[]) } else { - source_fin.open("Source.txt"); - outfile_path = "OutFile.csv"; - if (!source_fin.is_open()) - { - std::cout << "ERROR: Failed to open Source.txt" << std::endl; - return 1; - } - std::vector numbers; - double number; - while (source_fin >> number) - { - numbers.push_back(number); - } - plane_wave(arguments, numbers, outfile_path); - source_fin.close(); + std::cout << "ERROR: Incorrect number of parameters" << std::endl; + exit(1); } return 0; } diff --git a/sln/FDTD/CMakeLists.txt b/sln/FDTD/CMakeLists.txt index 18d6c9f..f1ff02b 100644 --- a/sln/FDTD/CMakeLists.txt +++ b/sln/FDTD/CMakeLists.txt @@ -1,3 +1,3 @@ set(target "FDTD") -add_library(${target} STATIC ../../include/FDTD.h ../../src/FDTD.cpp) +add_library(${target} STATIC ../../include/FDTD.h ../../include/Field.h ../../src/FDTD.cpp ../../src/Field.cpp) diff --git a/src/FDTD.cpp b/src/FDTD.cpp index 1e5125f..247792a 100644 --- a/src/FDTD.cpp +++ b/src/FDTD.cpp @@ -4,47 +4,7 @@ #include #include -//#define __DEBUG_PML__ -//#define __OUTPUT_SIGMA__ - -Field::Field(const int _Ni = 1, const int _Nj = 1, const int _Nk = 1) - : Ni(_Ni), Nj(_Nj), Nk(_Nk) -{ - int size = Ni * Nj * Nk; - field = std::vector(size, 0.0); -} - -Field& Field::operator= (const Field& other) -{ - if (this != &other) - { - field = other.field; - Ni = other.Ni; - Nj = other.Nj; - Nk = other.Nk; - } - return *this; -} - -double& Field::operator() (int i, int j, int k) -{ - int i_isMinusOne = (i == -1); - int j_isMinusOne = (j == -1); - int k_isMinusOne = (k == -1); - int i_isNi = (i == Ni); - int j_isNj = (j == Nj); - int k_isNk = (k == Nk); - - int truly_i = (Ni - 1) * i_isMinusOne + i * - !(i_isMinusOne || i_isNi); - int truly_j = (Nj - 1) * j_isMinusOne + j * - !(j_isMinusOne || j_isNj); - int truly_k = (Nk - 1) * k_isMinusOne + k * - !(k_isMinusOne || k_isNk); - - int index = truly_i + truly_j * Ni + truly_k * Ni * Nj; - return field[index]; -} +#include "Field.h" void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], double SGm, std::function dist) @@ -58,10 +18,10 @@ void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], { EsigmaX(i, j, k) = SGm * std::pow((static_cast(dist(i, j, k))) / - static_cast(pml_size_i), FDTDconst::N); + static_cast(pml_size_i), FDTDconst::N); BsigmaX(i, j, k) = SGm * std::pow((static_cast(dist(i, j, k))) / - static_cast(pml_size_i), FDTDconst::N); + static_cast(pml_size_i), FDTDconst::N); } } } @@ -107,7 +67,7 @@ void FDTD::set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], } } -FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) : +FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) : parameters(_parameters), dt(_dt), pml_percent(_pml_percent) { if (parameters.Ni <= 0 || @@ -117,13 +77,13 @@ FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) : { throw std::invalid_argument("ERROR: invalid parameters"); } - Ex = Ey = Ez = Bx = By = Bz + Ex = Ey = Ez = Bx = By = Bz = Field(parameters.Ni, parameters.Nj, parameters.Nk); - Exy = Exz = Eyx = Eyz = Ezx = Ezy + Exy = Exz = Eyx = Eyz = Ezx = Ezy = Field(parameters.Ni, parameters.Nj, parameters.Nk); - Bxy = Bxz = Byx = Byz = Bzx = Bzy + Bxy = Bxz = Byx = Byz = Bzx = Bzy = Field(parameters.Ni, parameters.Nj, parameters.Nk); - EsigmaX = EsigmaY = EsigmaZ = BsigmaX = BsigmaY = BsigmaZ + EsigmaX = EsigmaY = EsigmaZ = BsigmaX = BsigmaY = BsigmaZ = Field(parameters.Ni, parameters.Nj, parameters.Nk); pml_size_i = static_cast(static_cast(parameters.Ni) * pml_percent); @@ -178,15 +138,15 @@ void FDTD::update_E(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t) { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - Ex(i, j, k) = Ex(i, j, k) - 4.0 * FDTDconst::PI * dt * Jx[t](i, j, k) + - FDTDconst::C * dt * ((Bz(i, j, k) - Bz(i, j - 1, k)) / dy - - (By(i, j, k) - By(i, j, k - 1)) / dz); - Ey(i, j, k) = Ey(i, j, k) - 4.0 * FDTDconst::PI * dt * Jy[t](i, j, k) + - FDTDconst::C * dt * ((Bx(i, j, k) - Bx(i, j, k - 1)) / dz - - (Bz(i, j, k) - Bz(i - 1, j, k)) / dx); - Ez(i, j, k) = Ez(i, j, k) - 4.0 * FDTDconst::PI * dt * Jz[t](i, j, k) + - FDTDconst::C * dt * ((By(i, j, k) - By(i - 1, j, k)) / dx - - (Bx(i, j, k) - Bx(i, j - 1, k)) / dy); + Ex(i, j, k) = Ex(i, j, k) - 4.0 * FDTDconst::PI * dt * Jx[t](i, j, k) + + FDTDconst::C * dt * ((Bz(i, j, k) - Bz(i, j - 1, k)) / dy - + (By(i, j, k) - By(i, j, k - 1)) / dz); + Ey(i, j, k) = Ey(i, j, k) - 4.0 * FDTDconst::PI * dt * Jy[t](i, j, k) + + FDTDconst::C * dt * ((Bx(i, j, k) - Bx(i, j, k - 1)) / dz - + (Bz(i, j, k) - Bz(i - 1, j, k)) / dx); + Ez(i, j, k) = Ez(i, j, k) - 4.0 * FDTDconst::PI * dt * Jz[t](i, j, k) + + FDTDconst::C * dt * ((By(i, j, k) - By(i - 1, j, k)) / dx - + (Bx(i, j, k) - Bx(i, j - 1, k)) / dy); } } } @@ -205,15 +165,15 @@ void FDTD::update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2]) { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - Bx(i, j, k) += FDTDconst::C * dt / 2.0 * - ((Ey(i, j, k + 1) - Ey(i, j, k)) / dz - - (Ez(i, j + 1, k) - Ez(i, j, k)) / dy); - By(i, j, k) += FDTDconst::C * dt / 2.0 * - ((Ez(i + 1, j, k) - Ez(i, j, k)) / dx - - (Ex(i, j, k + 1) - Ex(i, j, k)) / dz); - Bz(i, j, k) += FDTDconst::C * dt / 2.0 * - ((Ex(i, j + 1, k) - Ex(i, j, k)) / dy - - (Ey(i + 1, j, k) - Ey(i, j, k)) / dx); + Bx(i, j, k) += FDTDconst::C * dt / 2.0 * + ((Ey(i, j, k + 1) - Ey(i, j, k)) / dz - + (Ez(i, j + 1, k) - Ez(i, j, k)) / dy); + By(i, j, k) += FDTDconst::C * dt / 2.0 * + ((Ez(i + 1, j, k) - Ez(i, j, k)) / dx - + (Ex(i, j, k + 1) - Ex(i, j, k)) / dz); + Bz(i, j, k) += FDTDconst::C * dt / 2.0 * + ((Ex(i, j + 1, k) - Ex(i, j, k)) / dy - + (Ey(i + 1, j, k) - Ey(i, j, k)) / dx); } } } @@ -234,7 +194,7 @@ void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) - { + { for (int j = bounds_j[0]; j < bounds_j[1]; j++) { for (int k = bounds_k[0]; k < bounds_k[1]; k++) @@ -284,7 +244,7 @@ void FDTD::update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) double dz = parameters.dz; double PMLcoef2_x, PMLcoef2_y, PMLcoef2_z; - + #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) { @@ -330,67 +290,13 @@ void FDTD::update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) } } -double FDTD::calc_reflection(Field E_start[3], Field B_start[3], - Field E_final[3], Field B_final[3]) -{ - double energy_start = 0.0; - double energy_final = 0.0; - -#pragma omp parallel for collapse(2) - for (int i = pml_size_i; i < parameters.Ni - pml_size_i; i++) - { - for (int j = pml_size_j; j < parameters.Nj - pml_size_j; j++) - { - for (int k = pml_size_k; k < parameters.Nk - pml_size_k; k++) - { - energy_start += std::pow(E_start[0](i, j, k), 2.0) + - std::pow(E_start[1](i, j, k), 2.0) + - std::pow(E_start[2](i, j, k), 2.0) + - std::pow(B_start[0](i, j, k), 2.0) + - std::pow(B_start[1](i, j, k), 2.0) + - std::pow(B_start[2](i, j, k), 2.0); - - energy_final += std::pow(E_final[0](i, j, k), 2.0) + - std::pow(E_final[1](i, j, k), 2.0) + - std::pow(E_final[2](i, j, k), 2.0) + - std::pow(B_final[0](i, j, k), 2.0) + - std::pow(B_final[1](i, j, k), 2.0) + - std::pow(B_final[2](i, j, k), 2.0); - } - } - } - - return std::sqrt(energy_final / energy_start); -} - -double FDTD::calc_max_value(Field& field) -{ - double max_val = 0.0; - -#pragma omp parallel for collapse(2) - for (int i = 0; i < parameters.Ni - pml_size_i; i++) - { - for (int j = 0; j < parameters.Nj - pml_size_j; j++) - { - for (int k = 0; k < parameters.Nk - pml_size_k; k++) - { - if (std::abs(field(i, j, k)) > std::abs(max_val)) - { - max_val = field(i, j, k); - } - } - } - } - return max_val; -} - -std::vector> FDTD::update_fields(const int time) +std::vector FDTD::update_fields(const int time, bool write_result, Axis write_axis, std::string base_path) { if (time < 0) { throw std::invalid_argument("ERROR: Invalid update field argument"); } - std::vector> return_data; + std::vector return_data; if (pml_percent == 0.0) { @@ -404,11 +310,13 @@ std::vector> FDTD::update_fields(const int time) update_B(size_i_main, size_j_main, size_k_main); std::vector new_iteration{ Ex, Ey, Ez, Bx, By, Bz }; - return_data.push_back(new_iteration); + return_data = new_iteration; + if (write_result) + write_spherical(return_data, write_axis, base_path, t); } return return_data; } - + // Defining areas of computation int size_i_main[] = { pml_size_i, parameters.Ni - pml_size_i }; int size_j_main[] = { pml_size_j, parameters.Nj - pml_size_j }; @@ -434,29 +342,29 @@ std::vector> FDTD::update_fields(const int time) int size_zx_upper_j_pml[] = { parameters.Nj - pml_size_j, parameters.Nj }; // Definition of functions for calculating the distance to the interface - std::function calc_distant_i_up = - [=](int i, int j, int k) { + std::function calc_distant_i_up = + [=](int i, int j, int k) { return i + 1 + pml_size_i - parameters.Ni; }; - std::function calc_distant_j_up = - [=](int i, int j, int k) { + std::function calc_distant_j_up = + [=](int i, int j, int k) { return j + 1 + pml_size_j - parameters.Nj; }; - std::function calc_distant_k_up = - [=](int i, int j, int k) { + std::function calc_distant_k_up = + [=](int i, int j, int k) { return k + 1 + pml_size_k - parameters.Nk; }; - std::function calc_distant_i_low = - [=](int i, int j, int k) { + std::function calc_distant_i_low = + [=](int i, int j, int k) { return pml_size_i - i; }; - std::function calc_distant_j_low = - [=](int i, int j, int k) { + std::function calc_distant_j_low = + [=](int i, int j, int k) { return pml_size_j - j; }; - std::function calc_distant_k_low = - [=](int i, int j, int k) { + std::function calc_distant_k_low = + [=](int i, int j, int k) { return pml_size_k - k; }; @@ -469,26 +377,19 @@ std::vector> FDTD::update_fields(const int time) / (static_cast(pml_size_k) * parameters.dz); // Calculation of permittivity in the cells - set_sigma_z(size_i_solid, size_j_solid, size_xy_lower_k_pml, - SGm_z, calc_distant_k_low); - set_sigma_y(size_i_solid, size_zx_lower_j_pml, size_k_solid, - SGm_y, calc_distant_j_low); - set_sigma_x(size_yz_lower_i_pml, size_j_solid, size_k_solid, - SGm_x, calc_distant_i_low); - - set_sigma_z(size_i_solid, size_j_solid, size_xy_upper_k_pml, - SGm_z, calc_distant_k_up); - set_sigma_y(size_i_solid, size_zx_upper_j_pml, size_k_solid, - SGm_y, calc_distant_j_up); - set_sigma_x(size_yz_upper_i_pml, size_j_solid, size_k_solid, - SGm_x, calc_distant_i_up); - -#ifdef __DEBUG_PML__ - double max_val_1 = 0.0; - double max_val_2 = 0.0; - int t_start = 50; - int t_final = 400; -#endif // __DEBUG_PML__ + set_sigma_z(size_i_solid, size_j_solid, size_xy_lower_k_pml, + SGm_z, calc_distant_k_low); + set_sigma_y(size_i_solid, size_zx_lower_j_pml, size_k_solid, + SGm_y, calc_distant_j_low); + set_sigma_x(size_yz_lower_i_pml, size_j_solid, size_k_solid, + SGm_x, calc_distant_i_low); + + set_sigma_z(size_i_solid, size_j_solid, size_xy_upper_k_pml, + SGm_z, calc_distant_k_up); + set_sigma_y(size_i_solid, size_zx_upper_j_pml, size_k_solid, + SGm_y, calc_distant_j_up); + set_sigma_x(size_yz_upper_i_pml, size_j_solid, size_k_solid, + SGm_x, calc_distant_i_up); for (int t = 0; t < time; t++) { @@ -516,57 +417,11 @@ std::vector> FDTD::update_fields(const int time) update_B(size_i_main, size_j_main, size_k_main); -#ifdef __DEBUG_PML__ - if (t == t_start) - { - max_val_1 = calc_max_value(Ex); - } - if (t == t_final) - { - max_val_2 = calc_max_value(Ex); - } - -#endif //__DEBUG_PML__ - -#ifndef __OUTPUT_SIGMA__ std::vector new_iteration{ Ex, Ey, Ez, Bx, By, Bz }; -#endif // __OUTPUT_SIGMA__ - -#ifdef __OUTPUT_SIGMA__ - - EsigmaX(parameters.Ni / 2, 1, 1) = 100; - - std::vector new_iteration{ EsigmaX, EsigmaY, EsigmaZ, - BsigmaX, BsigmaY, BsigmaZ }; - return_data.push_back(new_iteration); - if (t == 0) break; -#endif // __OUTPUT_SIGMA__ - - return_data.push_back(new_iteration); + return_data = new_iteration; + if (write_result) + write_spherical(return_data, write_axis, base_path, t); } -#ifdef __DEBUG_PML__ - Field E_start[] = { return_data[t_start][0], - return_data[t_start][1], - return_data[t_start][2] }; - Field B_start[] = { return_data[t_start][3], - return_data[t_start][4], - return_data[t_start][5] }; - - Field E_final[] = { return_data[t_final][0], - return_data[t_final][1], - return_data[t_final][2] }; - Field B_final[] = { return_data[t_final][3], - return_data[t_final][4], - return_data[t_final][5] }; - - std::cout << "reflection: " << - calc_reflection(E_start, B_start, E_final, B_final) - << std::endl; - - std::cout << "before: " << max_val_1 << std::endl; - std::cout << "after: " << max_val_2 << std::endl; -#endif // __DEBUG_PML__ - return return_data; } diff --git a/src/Field.cpp b/src/Field.cpp new file mode 100644 index 0000000..e0b2966 --- /dev/null +++ b/src/Field.cpp @@ -0,0 +1,41 @@ +#include "Field.h" + +Field::Field(const int _Ni, const int _Nj, const int _Nk) + : Ni(_Ni), Nj(_Nj), Nk(_Nk) +{ + int size = Ni * Nj * Nk; + field = std::vector(size, 0.0); +} + +Field& Field::operator= (const Field& other) +{ + if (this != &other) + { + field = other.field; + Ni = other.Ni; + Nj = other.Nj; + Nk = other.Nk; + } + return *this; +} + +double& Field::operator() (int i, int j, int k) +{ + int i_isMinusOne = (i == -1); + int j_isMinusOne = (j == -1); + int k_isMinusOne = (k == -1); + int i_isNi = (i == Ni); + int j_isNj = (j == Nj); + int k_isNk = (k == Nk); + + int truly_i = (Ni - 1) * i_isMinusOne + i * + !(i_isMinusOne || i_isNi); + int truly_j = (Nj - 1) * j_isMinusOne + j * + !(j_isMinusOne || j_isNj); + int truly_k = (Nk - 1) * k_isMinusOne + k * + !(k_isMinusOne || k_isNk); + + int index = truly_i + truly_j * Ni + truly_k * Ni * Nj; + return field[index]; +} + diff --git a/src/Writer.cpp b/src/Writer.cpp new file mode 100644 index 0000000..cef8cc1 --- /dev/null +++ b/src/Writer.cpp @@ -0,0 +1,141 @@ +#include "Writer.h" + +void write_x(Field& this_field, std::ofstream& fout) +{ + for (int k = 0; k < this_field.get_Nk(); ++k) + { + for (int j = 0; j < this_field.get_Nj(); ++j) + { + for (int i = 0; i < this_field.get_Ni(); ++i) + { + fout << this_field(i, j, k); + if (i == this_field.get_Ni() - 1) + { + fout << std::endl; + } + else + { + fout << ";"; + } + } + } + } + fout << std::endl << std::endl; +} +void write_y(Field& this_field, std::ofstream& fout) +{ + for (int k = 0; k < this_field.get_Nk(); ++k) + { + for (int i = 0; i < this_field.get_Ni(); ++i) + { + for (int j = 0; j < this_field.get_Nj(); ++j) + { + fout << this_field(i, j, k); + if (j == this_field.get_Nj() - 1) + { + fout << std::endl; + } + else + { + fout << ";"; + } + } + } + } + fout << std::endl << std::endl; +} +void write_z(Field& this_field, std::ofstream& fout) +{ + for (int i = 0; i < this_field.get_Ni(); ++i) + { + for (int j = 0; j < this_field.get_Nj(); ++j) + { + for (int k = 0; k < this_field.get_Nk(); ++k) + { + fout << this_field(i, j, k); + if (k == this_field.get_Nk() - 1) + { + fout << std::endl; + } + else + { + fout << ";"; + } + } + } + } + fout << std::endl << std::endl; +} + +void write_spherical(std::vector fields, Axis axis, std::string base_path, int it) +{ + std::ofstream test_fout; + + for (int c = static_cast(Component::EX); c <= static_cast(Component::BZ); ++c) + { + std::cout << "Write : " << it << " -- " << c << std::endl; + if (axis == Axis::X) + { + test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv"); + Field field = fields[c]; + for (int k = 0; k < field.get_Nk(); ++k) + { + for (int j = 0; j < field.get_Nj(); ++j) + { + test_fout << field(field.get_Ni() / 2, j, k); + if (j == field.get_Nj() - 1) + { + test_fout << std::endl; + } + else + { + test_fout << ";"; + } + } + } + test_fout.close(); + } + else if (axis == Axis::Y) + { + test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv"); + Field field = fields[c]; + for (int i = 0; i < field.get_Nj(); ++i) + { + for (int k = 0; k < field.get_Nk(); ++k) + { + test_fout << field(i, field.get_Nj() / 2, k); + if (k == field.get_Nk() - 1) + { + test_fout << std::endl; + } + else + { + test_fout << ";"; + } + } + } + test_fout.close(); + } + else + { + test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv"); + Field field = fields[c]; + for (int i = 0; i < field.get_Nj(); ++i) + { + for (int j = 0; j < field.get_Nk(); ++j) + { + test_fout << field(i, j, field.get_Nk() / 2); + if (j == field.get_Nk() - 1) + { + test_fout << std::endl; + } + else + { + test_fout << ";"; + } + } + } + test_fout.close(); + } + } +}