fix 13: add version with shifts

PlotScript/Source.txt

@@ -1,8 +1,8 @@
 32
-17
+32
 0
 1
 0
-0.5
-5e-17
-50e-13
+1
+5e-16
+5e-12

PlotScript/main.py

@@ -5,7 +5,7 @@
 
 components = {1: "Ex", 2: "Ey", 3: "Ez", 4: "Bx", 5: "By", 6: "Bz"}
 nums_com = {"Ex": 0, "Ey": 1, "Ez": 2, "Bx": 3, "By": 4, "Bz": 5}
-
+shift_flag = "1"
 
 def get_plot(num, data, size_n, size_x, time):
     x = 0.
@@ -26,7 +26,7 @@ def get_plot(num, data, size_n, size_x, time):
     plt.plot(X, V)
     plt.xlabel('X')
     plt.ylabel(components[num])
-    plt.title('Plot ' + components[num] + " (Time: " + str(time) + ")")
+    plt.title(f"Plot {components[num]} (Time: {str(time)})")
     plt.grid(True)
     plt.savefig(f"Plots/plt_{components[num]}")
     plt.show()
@@ -40,7 +40,7 @@ def execute_cpp(field_1, field_2, field_to_plot, source_nums):
     print("\n" + field_to_plot + ":\n")
 
     cpp_executable = "src/sample.exe"
-    args = [cpp_executable, str(num_field_1), str(num_field_2), str(num_field_to_plot)]
+    args = [cpp_executable, str(num_field_1), str(num_field_2), str(num_field_to_plot), shift_flag]
     try:
         subprocess.run(args, check=True)
     except subprocess.CalledProcessError:
@@ -91,7 +91,7 @@ def execute_cpp(field_1, field_2, field_to_plot, source_nums):
 
     execute_cpp("Ex", "Bz", "Ex", numbers)
     execute_cpp("Ey", "Bz", "Ey", numbers)
-    execute_cpp("Ez", "Bx", "Ez", numbers)
+    execute_cpp("Ez", "By", "Ez", numbers)
     execute_cpp("Ez", "Bx", "Bx", numbers)
     execute_cpp("Ez", "By", "By", numbers)
-    execute_cpp("Ey", "Bz", "Bz", numbers)
+    execute_cpp("Ex", "Bz", "Bz", numbers)

include/FDTD.h

@@ -12,8 +12,9 @@ class Cell_number
 private:
     int border;
     int current;
+    int start;
 public:
-    explicit Cell_number(int max_num, int start_num = 0) : border(max_num), current(start_num) {}
+    explicit Cell_number(int max_num, int start_num = 0) : border(max_num), current(start_num), start(start_num) {}
     int operator+ (int) const;
     int operator- (int) const;
     int operator* ();
@@ -53,4 +54,5 @@ class FDTD
     Field& get_field(Component);
 
     void update_field(const double&);
+    void shifted_update_field(const double&);
 };

samples/sample.cpp

@@ -1,31 +1,61 @@
 #define _USE_MATH_DEFINES
+//#define __DEBUG__
 
 #include <iostream>
 #include <iomanip>
 #include <fstream>
+
 #include "FDTD.h"
 
-void initial_filling_x(FDTD& test, Component field_1, Component field_2, int size_N[2], double size_d, double size_wave[2])
+double sign;
+double shift = 0.0;
+
+double get_vector_sign(Component vect_1, Component vect_2)
+{
+    if (vect_1 == Component::EY && vect_2 == Component::BZ || vect_1 == Component::EZ && vect_2 == Component::BX)
+    {
+        return 1.0;
+    }
+    else if (vect_1 == Component::EX && vect_2 == Component::BZ || vect_1 == Component::EZ && vect_2 == Component::BY)
+    {
+        return -1.0;
+    }
+    else return 0.0;
+}
+
+void initial_filling_x(FDTD& test, Component field_1, Component field_2, int size_N[2], double size_d, double size_wave[2], bool shifted = true)
 {
-    double x = 0;
+    double x = 0.0;
+    double x_b = 0.0;
     for (int i = 0; i < size_N[0]; x += size_d, ++i)
     {
         for (int j = 0; j < size_N[1]; ++j)
         {
-            test.get_field(field_1)(i, j) = sin(2 * M_PI * (x - size_wave[0]) / (size_wave[1] - size_wave[0]));
-            test.get_field(field_2)(i, j) = sin(2 * M_PI * (x - size_wave[0]) / (size_wave[1] - size_wave[0]));
+            test.get_field(field_1)(i, j) = sign * sin(2.0 * M_PI * (x - size_wave[0]) / (size_wave[1] - size_wave[0]));
+            if (shifted)
+            {
+                x_b = x + size_d / 2.0;
+            }
+            else x_b = x;
+            test.get_field(field_2)(i, j) = sin(2.0 * M_PI * (x_b - size_wave[0]) / (size_wave[1] - size_wave[0]));
         }
     }
 }
-void initial_filling_y(FDTD& test, Component field_1, Component field_2, int size_N[2], double size_d, double size_wave[2])
+void initial_filling_y(FDTD& test, Component field_1, Component field_2, int size_N[2], double size_d, double size_wave[2], bool shifted = true)
 {
-    double y = 0;
+    double y = 0.0;
+    double y_b = 0.0;
     for (int j = 0; j < size_N[1]; y += size_d, ++j)
     {
         for (int i = 0; i < size_N[0]; ++i)
         {
-            test.get_field(field_1)(i, j) = sin(2 * M_PI * (y - size_wave[0]) / (size_wave[1] - size_wave[0]));
-            test.get_field(field_2)(i, j) = sin(2 * M_PI * (y - size_wave[0]) / (size_wave[1] - size_wave[0]));
+            test.get_field(field_1)(i, j) = sign * sin(2.0 * M_PI * (y - size_wave[0]) / (size_wave[1] - size_wave[0]));
+            if (shifted)
+            {
+                y_b = y + size_d / 2.0;
+            }
+            else y_b = y;
+            test.get_field(field_2)(i, j) = sin(2.0 * M_PI * (y_b - size_wave[0]) / (size_wave[1] - size_wave[0]));
         }
     }
 }
@@ -69,27 +99,27 @@ void write_y(Field& this_field, std::ofstream& fout)
     fout << std::endl << std::endl;
 }
 
-double max_abs_error_x(Field& this_field, int size_N[2], double size_x[2], double size_d[2], double t)
+double max_abs_error_x(Field& this_field, int size_N[2], double size_x[2], double size_d[2], double t, bool shifted = true)
 {
     double extr_n = 0.0;
-    double x = 0;
+    double x = (shifted) ? shift : 0.0;
     int j = 0;
     for (int i = 0; i < this_field.get_Ni(); ++i, x += size_d[0])
     {
-        double this_n = fabs(fabs(this_field(i, j)) - fabs(sin(2 * M_PI * (x - size_x[0] - FDTD_Const::C * t) / (size_x[1] - size_x[0]))));
+        double this_n = fabs(this_field(i, j) - sign * sin(2.0 * M_PI * (x - size_x[0] - FDTD_Const::C * t) / (size_x[1] - size_x[0])));
         if (this_n > extr_n)
             extr_n = this_n;
     }
     return extr_n;
 }
-double max_abs_error_y(Field& this_field, int size_N[2], double size_y[2], double size_d[2], double t)
+double max_abs_error_y(Field& this_field, int size_N[2], double size_y[2], double size_d[2], double t, bool shifted = true)
 {
     double extr_n = 0.0;
-    double y = 0;
+    double y = (shifted) ? shift : 0.0;
     int i = 0;
     for (int j = 0; j < this_field.get_Nj(); ++j, y += size_d[1])
     {
-        double this_n = fabs(fabs(this_field(i, j)) - fabs(sin(2 * M_PI * (y - size_y[0] - FDTD_Const::C * t) / (size_y[1] - size_y[0]))));
+        double this_n = fabs(this_field(i, j) - sign * sin(2.0 * M_PI * (y - size_y[0] - FDTD_Const::C * t) / (size_y[1] - size_y[0])));
         if (this_n > extr_n)
             extr_n = this_n;
     }
@@ -99,7 +129,13 @@ double max_abs_error_y(Field& this_field, int size_N[2], double size_y[2], doubl
 void write_all(FDTD& test, char axis)
 {
     std::ofstream test_fout;
-    test_fout.open("../../../PlotScript/OutFile.csv"); //
+
+#ifndef __DEBUG__
+    test_fout.open("OutFile.csv");
+#else
+    test_fout.open("../../../PlotScript/OutFile.csv");
+#endif
+
     if (!test_fout.is_open())
     {
         std::cout << "ERROR: Failed to open OutFile.csv" << std::endl;
@@ -116,21 +152,34 @@ void write_all(FDTD& test, char axis)
     test_fout.close();
 }
 
-int main() //int argc, char* argv[]
+#ifndef __DEBUG__
+int main(int argc, char* argv[])
+#else
+int main()
+#endif
 {
     std::vector<double> numbers;
     std::ifstream source_fin;
-    source_fin.open("../../../PlotScript/Source.txt"); //
+
+#ifndef __DEBUG__
+    source_fin.open("Source.txt");
+#else
+    source_fin.open("../../../PlotScript/Source.txt");
+#endif
 
     if (!source_fin.is_open())
     {
         std::cout << "ERROR: Failed to open Source.txt" << std::endl;
         return 1;
     }
-    const char* argv[4] = { "0", "2", "3", "2" };
+
+#ifdef __DEBUG__
+    const char* argv[5] = {"0", "0", "5", "0", "0"};
+#endif
 
     double number;
-    while (source_fin >> number) {
+    while (source_fin >> number) 
+    {
         numbers.push_back(number);
     }
 
@@ -145,39 +194,58 @@ int main() //int argc, char* argv[]
     Component fld_2 = static_cast<Component>(std::atoi(argv[2]));
     Component fld_3 = static_cast<Component>(std::atoi(argv[3]));
 
+    bool shifted_version = std::atoi(argv[4]);
+    sign = get_vector_sign(fld_1, fld_2);
+
     char selected_axis;
     if (fld_1 == Component::EY && fld_2 == Component::BZ || fld_1 == Component::EZ && fld_2 == Component::BY)
     {
         selected_axis = 'x';
 
-        initial_filling_x(test_1, fld_1, fld_2, arr_N, arr_d[0], arr_x);
+        initial_filling_x(test_1, fld_1, fld_2, arr_N, arr_d[0], arr_x, shifted_version);
 
-        test_1.update_field(numbers[7]);
+        if (shifted_version)
+        {
+            test_1.shifted_update_field(numbers[7]);
+        }
+        else test_1.update_field(numbers[7]);
 
-        std::cout << "Maximum absolute error: " << max_abs_error_x(test_1.get_field(fld_3), arr_N, arr_x, arr_d, numbers[7]) << std::endl;
+        if (static_cast<int>(fld_3) > static_cast<int>(Component::EZ))
+        {
+            sign = 1.0;
+            shift = arr_d[0] / 2;
+        }
+        std::cout << "Maximum absolute error: " << max_abs_error_x(test_1.get_field(fld_3), arr_N, arr_x, arr_d, numbers[7], shifted_version) << std::endl;
     }
     else if (fld_1 == Component::EX && fld_2 == Component::BZ || fld_1 == Component::EZ && fld_2 == Component::BX)
     {
         selected_axis = 'y';
 
-        initial_filling_y(test_1, fld_1, fld_2, arr_N, arr_d[1], arr_y);
+        initial_filling_y(test_1, fld_1, fld_2, arr_N, arr_d[1], arr_y, shifted_version);
 
-        test_1.update_field(numbers[7]);
+        if (shifted_version)
+        {
+            test_1.shifted_update_field(numbers[7]);
+        }
+        else test_1.update_field(numbers[7]);
 
-        std::cout << "Maximum absolute error: " << max_abs_error_y(test_1.get_field(fld_3), arr_N, arr_y, arr_d, numbers[7]) << std::endl;
+        if (static_cast<int>(fld_3) > static_cast<int>(Component::EZ))
+        {
+            sign = 1.0;
+            shift = arr_d[1] / 2;
+        }
+        std::cout << "Maximum absolute error: " << max_abs_error_y(test_1.get_field(fld_3), arr_N, arr_y, arr_d, numbers[7], shifted_version) << std::endl;
     }
     else
     {
         std::cout << "ERROR" << std::endl;
         exit(1);
     }
-
     write_all(test_1, selected_axis);
 
     source_fin.close();
 
     std::cout << "The work of sample.exe is completed without errors" << std::endl;
 
-
     return 0;
 }

src/FDTD.cpp

@@ -4,7 +4,7 @@ int Cell_number::operator+ (int number) const
 {
     if (current + number >= border)
     {
-        return 0;
+        return -1 + number;
     }
     else
     {
@@ -22,24 +22,28 @@ int Cell_number::operator- (int number) const
         return current - number;
     }
 }
+
 int Cell_number::operator* ()
 {
     return current;
 }
+
 Cell_number& Cell_number::operator++ ()
 {
     ++current;
     if (current > border)
     {
-        current = 0;
+        current = start;
     }
     return *this;
 }
+
 bool Cell_number::operator< (int other)
 {
     return current < other;
 }
 
+
 Field::Field(const int _Ni = 1, const int _Nj = 1) : Ni(_Ni), Nj(_Nj)
 {
     int size = Ni * Nj;
@@ -63,6 +67,7 @@ double& Field::operator() (int _i, int _j)
     return field[index];
 }
 
+
 FDTD::FDTD(int size_grid[2], double size_x[2], double size_y[2], double _dt) : dt(_dt)
 {
     Ni = size_grid[0];
@@ -122,3 +127,31 @@ void FDTD::update_field(const double& time)
         }
     }
 }
+void FDTD::shifted_update_field(const double& time)
+{
+    for (double t = 0; t < time; t += dt)
+    {
+        for (Cell_number j(Nj); j < Nj; ++j)
+        {
+            for (Cell_number i(Ni); i < Ni; ++i)
+            {
+                Ex(*i, *j) += FDTD_Const::C * dt * (Bz(*i, *j) - Bz(*i, j - 1)) / dy;
+
+                Ey(*i, *j) -= FDTD_Const::C * dt * (Bz(*i, *j) - Bz(i - 1, *j)) / dx;
+
+                Ez(*i, *j) += FDTD_Const::C * dt * ((By(*i, *j) - By(i - 1, *j)) / dx - (Bx(*i, *j) - Bx(*i, j - 1)) / dy);
+            }
+        }
+        for (Cell_number j(Nj); j < Nj; ++j)
+        {
+            for (Cell_number i(Ni); i < Ni; ++i)
+            {
+                Bx(*i, *j) -= FDTD_Const::C * dt * (Ez(*i, j + 1) - Ez(*i, *j)) / dy;
+
+                By(*i, *j) += FDTD_Const::C * dt * (Ez(i + 1, *j) - Ez(*i, *j)) / dx;
+
+                Bz(*i, *j) -= FDTD_Const::C * dt * ((Ey(i + 1, *j) - Ey(*i, *j)) / dx - (Ex(*i, j + 1) - Ex(*i, *j)) / dy);
+            }
+        }
+    }
+}