From eb8b8566c6e6466d5fd4af75b1a68c740230dac9 Mon Sep 17 00:00:00 2001 From: Cory Date: Mon, 1 Jan 2024 15:57:03 +0100 Subject: [PATCH] Make some changes --- lib/WaveSimulation.cpp | 25 +++++++++++++++++++----- lib/utils.cpp | 43 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 5 deletions(-) diff --git a/lib/WaveSimulation.cpp b/lib/WaveSimulation.cpp index a78fd52..83f5888 100644 --- a/lib/WaveSimulation.cpp +++ b/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 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 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 &steps) +void WaveSimulation::simulate(std::string outfile, std::vector &steps) { // Create path and proceed if successful. if (!utils::mkpath(utils::dirname(outfile))) { diff --git a/lib/utils.cpp b/lib/utils.cpp index 89752ae..ae07678 100644 --- a/lib/utils.cpp +++ b/lib/utils.cpp @@ -113,4 +113,47 @@ void print_sp_matrix_structure(const arma::sp_cx_mat &A) cout << endl; } +std::vector split(const std::string &s, char delim) { + std::vector 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