lib/WaveSimulation.cpp

@@ -27,10 +27,10 @@ void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x,
     double x, y, diff_x, diff_y;
     std::complex<double> sum = 0.;
     for (size_t j = 0; j < this->U.n_cols; j++) {
-        x = j * h;
+        x = (j+1) * h;
         diff_x = x - x_c;
         for (size_t i = 0; i < this->U.n_rows; i++) {
-            y = i * h;
+            y = (i+1) * h;
             diff_y = y - y_c;
             this->U(i, j) =
                 std::exp(-(diff_x * diff_x) / (2. * sigma_x * sigma_x)
@@ -117,7 +117,8 @@ void WaveSimulation::initialize_V(double thickness, double pos_x, double ap_sep,
     res = arma::join_cols(res, arma::cx_vec(fill + ((this->N - res.size()) % 2),
                                             arma::fill::value(1e10)));
 
-    uint32_t start = pos_x / this->h - thickness / this->h / 2;
+    // Subtract one to compensate for the border being gone
+    uint32_t start = (pos_x / this->h - 1) - thickness / this->h / 2;
     for (size_t i = 0; i < thickness / this->h; i++) {
         this->V.col(start + i) = res;
     }
@@ -125,32 +126,40 @@ void WaveSimulation::initialize_V(double thickness, double pos_x, double ap_sep,
 
 void WaveSimulation::initialize_A()
 {
+    DEBUG("Inside initialize_A");
     // Create the diagonal
     arma::cx_vec diagonal(this->N * this->N);
 
+    DEBUG("Before setting diagonal values");
     // Set diagonal values
     std::complex<double> r = (1._i * this->dt) / (2 * h * h);
     for (size_t i = 0; i < diagonal.size(); i++) {
         diagonal(i) = 1. + 4. * r + (1._i * this->dt / 2.) * this->V(i);
     }
 
+    DEBUG("Create submatrix");
     // Create the submatrix
     arma::cx_mat sub_matrix(this->N, this->N, arma::fill::zeros);
     sub_matrix.diag(-1).fill(-r);
     sub_matrix.diag(1).fill(-r);
 
+    DEBUG("Set the size of A");
     // Set the size of A
     this->A.set_size(this->N * this->N, this->N * this->N);
 
+    DEBUG("Fill in values");
     // Fill in the values in the submatrix diagonal
     for (size_t i = 0; i < this->A.n_cols; i += this->N) {
         this->A.submat(i, i, i + this->N - 1, i + this->N - 1) = sub_matrix;
     }
 
+    DEBUG("Fill in diagonals");
     // Fill the last sub/sup-diagonals
     this->A.diag() = diagonal;
     this->A.diag(-this->N).fill(-r);
+    DEBUG("After diagonal");
     this->A.diag(this->N).fill(-r);
+    DEBUG("After");
 }
 
 void WaveSimulation::initialize_B()
@@ -205,15 +214,21 @@ WaveSimulation::WaveSimulation(double h, double dt, double T, double x_c,
                                double y_c, double sigma_x, double sigma_y,
                                double p_x, double p_y)
 {
+    DEBUG("Inside constructor");
     this->dt = dt;
     this->h = h;
     this->T = T;
     this->M = 1. / h;
     this->N = M - 2;
+    DEBUG("Before V");
     this->initialize_V();
+    DEBUG("Before U");
     this->initialize_U(x_c, y_c, sigma_x, sigma_y, p_x, p_y);
+    DEBUG("Before A");
     this->initialize_A();
+    DEBUG("Before B");
     this->initialize_B();
+    DEBUG("After B");
 }
 
 // Public methods
@@ -226,7 +241,7 @@ void WaveSimulation::step()
     arma::spsolve(this->U, this->A, tmp);
 }
 
-void WaveSimulation::solve(std::string outfile, bool write_each_step)
+void WaveSimulation::simulate(std::string outfile, bool write_each_step)
 {
     // Create path and proceed if successful.
     if (!utils::mkpath(utils::dirname(outfile))) {
@@ -256,7 +271,7 @@ void WaveSimulation::solve(std::string outfile, bool write_each_step)
     this->write_U(ofile);
     ofile.close();
 }
-void WaveSimulation::solve(std::string outfile, std::vector<double> &steps)
+void WaveSimulation::simulate(std::string outfile, std::vector<double> &steps)
 {
     // Create path and proceed if successful.
     if (!utils::mkpath(utils::dirname(outfile))) {

lib/utils.cpp

@@ -113,4 +113,47 @@ void print_sp_matrix_structure(const arma::sp_cx_mat &A)
     cout << endl;
 }
 
+std::vector<std::string> split(const std::string &s, char delim) {
+    std::vector<std::string> result;
+    std::stringstream ss(s);
+    std::string item;
+
+    while (getline(ss, item, delim)) {
+        result.push_back(trim(item));
+    }
+
+    return result;
+}
+
+inline std::string& ltrim(std::string& s, const char* t)
+{
+    s.erase(0, s.find_first_not_of(t));
+    return s;
+}
+
+inline std::string& rtrim(std::string& s, const char* t)
+{
+    s.erase(s.find_last_not_of(t) + 1);
+    return s;
+}
+
+inline std::string& trim(std::string& s, const char* t)
+{
+    return ltrim(rtrim(s, t), t);
+}
+
+inline std::string ltrim_copy(std::string s, const char* t)
+{
+    return ltrim(s, t);
+}
+
+inline std::string rtrim_copy(std::string s, const char* t)
+{
+    return rtrim(s, t);
+}
+
+inline std::string trim_copy(std::string s, const char* t)
+{
+    return trim(s, t);
+}
 } // namespace utils