Minimum, messy code

lib/WaveSimulation.cpp

@@ -0,0 +1,180 @@
+/** @file WaveSimulation.cpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 0.1
+ *
+ *  @brief Implementation of the WaveSimulation class.
+ *
+ *  @bug No known bugs
+ * */
+#include "WaveSimulation.hpp"
+
+#include "utils.hpp"
+
+#include <cmath>
+#include <complex>
+#include <cstdint>
+#include <cstdlib>
+
+WaveSimulation::WaveSimulation(double h, double dt, double T)
+{
+    this->dt = dt;
+    this->h = h;
+    this->T = T;
+    this->M = 1. / h;
+    this->N = M - 2;
+    this->V.set_size(this->N, this->N);
+    this->V.fill(0.);
+    this->U.set_size(this->N, this->N);
+    this->U.fill(0.);
+}
+
+void WaveSimulation::solve(std::ofstream &ofile)
+{
+    ofile << this->N << '\n';
+    uint32_t iterations = this->T / this->dt;
+    for (size_t i = 0; i < iterations; i++) {
+        this->write_U(ofile);
+        this->step();
+    }
+}
+
+void WaveSimulation::step()
+{
+    DEBUG("Inside step");
+    arma::cx_vec tmp = this->B * this->U.as_col();
+    arma::spsolve(this->U, this->A, tmp);
+}
+
+void WaveSimulation::build_A()
+{
+    // Create the diagonal
+    arma::cx_vec diagonal(this->N * this->N);
+
+    // 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);
+    }
+
+    // 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);
+
+    // Set the size of A
+    this->A.set_size(this->N * this->N, this->N * this->N);
+
+    // 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;
+    }
+
+    // Fill the last sub/sup-diagonals
+    this->A.diag() = diagonal;
+    this->A.diag(-this->N).fill(-r);
+    this->A.diag(this->N).fill(-r);
+}
+
+void WaveSimulation::build_B()
+{
+    std::complex<double> r = (1._i * this->dt) / (2 * h * h);
+
+    // Create the diagonal
+    arma::cx_vec diagonal(this->N * this->N);
+
+    for (size_t i = 0; i < diagonal.size(); i++) {
+        diagonal(i) = 1. - 4. * r - (1._i * this->dt / 2.) * this->V(i);
+    }
+
+    // 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);
+
+    // Set the size of B
+    this->B.set_size(this->N * this->N, this->N * this->N);
+
+    // Fill in the values in the submatrix diagonal
+    for (size_t i = 0; i < this->B.n_cols; i += this->N) {
+        this->B.submat(i, i, i + this->N - 1, i + this->N - 1) = sub_matrix;
+    }
+
+    // Fill the last sub/sup-diagonals
+    this->B.diag() = diagonal;
+    this->B.diag(-this->N).fill(r);
+    this->B.diag(this->N).fill(r);
+}
+
+void WaveSimulation::initialize_U(double x_c, double y_c, double sigma_x,
+                                  double sigma_y, double p_x, double p_y)
+{
+    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;
+        diff_x = x - x_c;
+        for (size_t i = 0; i < this->U.n_rows; i++) {
+            y = i * h;
+            diff_y = y - y_c;
+            this->U(i, j) =
+                std::exp(-(diff_x * diff_x) / (2. * sigma_x * sigma_x)
+                         - (diff_y * diff_y) / (2. * sigma_y * sigma_y)
+                         + p_x * x * 1._i + p_y * y * 1._i);
+            sum += this->U(i, j) * std::conj(this->U(i, j));
+        }
+    }
+    if (std::abs(sum.imag()) > 1e-7) {
+        abort();
+    }
+    double norm = 1. / std::sqrt(sum.real());
+
+    this->U.for_each([norm](std::complex<double> &el) { el *= norm; });
+}
+
+void WaveSimulation::write_U(std::ofstream &ofile)
+{
+    this->U.for_each(
+        [&ofile](std::complex<double> el) { ofile << el << '\t'; });
+    ofile << '\n';
+}
+
+void WaveSimulation::build_V(double thickness, double pos_x,
+                             double aperture_separation, double aperture,
+                             uint32_t slits)
+{
+    uint32_t mid_y = this->N / 2 - (this->N % 2 == 0);
+    arma::cx_vec res;
+
+    if (slits % 2 == 0) {
+        res = arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10));
+        for (size_t i=0; i < slits; i+=2) {
+            res = arma::join_cols(res, arma::cx_vec(aperture/this->h,arma::fill::zeros));
+            res = arma::join_cols(arma::cx_vec(aperture/this->h,arma::fill::zeros), res);
+            res = arma::join_cols(res, arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)));
+            res = arma::join_cols(arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)), res);
+        }
+    }
+    else {
+        res = arma::cx_vec(aperture/this->h,arma::fill::value(0));
+        for (size_t i=0; i < slits-1; i+=2) {
+            res = arma::join_cols(res, arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)));
+            res = arma::join_cols(arma::cx_vec(aperture_separation/this->h,arma::fill::value(1e10)), res);
+            res = arma::join_cols(res, arma::cx_vec(aperture/this->h,arma::fill::zeros));
+            res = arma::join_cols(arma::cx_vec(aperture/this->h,arma::fill::zeros), res);
+        }
+    }
+    if (res.size() > this->N) {
+        abort();
+    }
+    uint32_t fill = (this->N - res.size()) / 2;
+    res = arma::join_cols(arma::cx_vec(fill, arma::fill::value(1e10)), res);
+    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;
+    for (size_t i=0; i < thickness/this->h; i++) {
+        this->V.col(start+i) = res;
+    }
+}

lib/literals.cpp

@@ -0,0 +1,17 @@
+/** @file literals.cpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 1.0
+ *
+ *  @brief The implementation of the literals.
+ *
+ *  @bug No known bugs
+ * */
+#include "literals.hpp"
+
+std::complex<double> operator""_i(long double magnitude)
+{
+    return std::complex<double>(0.,magnitude);
+}

lib/utils.cpp

@@ -70,4 +70,47 @@ std::string concatpath(const std::string &left, const std::string &right)
     }
 }
 
+void print_sp_matrix_structure(const arma::sp_cx_mat &A)
+{
+    using namespace std;
+    using namespace arma;
+
+    // Declare a C-style 2D array of strings.
+    string S[A.n_rows][A.n_cols];
+
+    // Initialise all the strings to " ".
+    for (int i = 0; i < A.n_rows; i++) {
+        for (int j = 0; j < A.n_cols; j++) {
+            S[i][j] = " ";
+        }
+    }
+
+    // Next, we want to set the string to a dot at each non-zero element.
+    // To do this we use the special loop iterator from the sp_cx_mat class
+    // to help us loop over only the non-zero matrix elements.
+    sp_cx_mat::const_iterator it = A.begin();
+    sp_cx_mat::const_iterator it_end = A.end();
+
+    int nnz = 0;
+    for (; it != it_end; ++it) {
+        S[it.row()][it.col()] = "•";
+        nnz++;
+    }
+
+    // Finally, print the matrix to screen.
+    cout << endl;
+    for (int i = 0; i < A.n_rows; i++) {
+        cout << "| ";
+        for (int j = 0; j < A.n_cols; j++) {
+            cout << S[i][j] << " ";
+        }
+        cout << "|\n";
+    }
+
+    cout << endl;
+    cout << "matrix size: " << A.n_rows << "x" << A.n_cols << endl;
+    cout << "non-zero elements: " << nnz << endl;
+    cout << endl;
+}
+
 } // namespace utils