Merge branch '9-add-easy-setting-of-boundary-conditions' into 'main'

Resolve "Add simplified setting of boundary conditions"

Closes #9

See merge request mluebke/diffusion!20

.gitlab-ci.yml

@@ -3,36 +3,51 @@ image: sobc/gitlab-ci
 stages:
     - build
     - test
+    - run
     - static_analyze
     - dynamic_analyze
 
 before_script:
     - apt-get update && apt-get install -y libeigen3-dev
 
-build:
+build_debug:
+    stage: build
+    artifacts:
+        paths:
+            - debug/app/1D
+            - debug/app/2D
+        expire_in: 100s
+    script:
+        - mkdir debug && cd debug
+        - cmake -DCMAKE_BUILD_TYPE=Debug -DENABLE_TESTING=ON ..
+        - make -j$(nproc)
+
+build_release:
     stage: build
     artifacts:
         paths:
             - build/app/1D
             - build/app/2D
             - build/app/Comp2D
+            - build/test/test
         expire_in: 100s
     script:
         - mkdir build && cd build
-        - cmake -DCMAKE_BUILD_TYPE=Debug ..
+        - cmake -DCMAKE_BUILD_TYPE=GenericOpt -DENABLE_TESTING=ON ..
-        - make
+        - make -j$(nproc)
+
+test:
+    stage: test
+    script:
+        - ./build/test/test
 
 run_1D:
-    stage: test
+    stage: run
-    dependencies: 
-        - build
     script:
         - ./build/app/1D
 
 run_2D:
-    stage: test
+    stage: run
-    dependencies: 
-       - build
     script:
         - ./build/app/2D
 
@@ -46,11 +61,11 @@ lint:
 memcheck_1D:
     stage: dynamic_analyze
     script:
-        - cd build/app
+        - cd debug/app
         - valgrind --tool=memcheck --leak-check=full --show-leak-kinds=all --track-origins=yes ./1D 1>/dev/null
 
 memcheck_2D:
     stage: dynamic_analyze
     script:
-        - cd build/app
+        - cd debug/app
         - valgrind --tool=memcheck --leak-check=full --show-leak-kinds=all --track-origins=yes ./2D 1>/dev/null

CMakeLists.txt

@@ -29,5 +29,13 @@ if(USE_UNSAFE_MATH_OPT)
   add_compile_options(-ffast-math)
 endif()
 
+option(ENABLE_TESTING
+  "Run tests after succesfull compilation"
+  OFF)
+
 add_subdirectory(app)
 add_subdirectory(src)
+
+if(ENABLE_TESTING)
+  add_subdirectory(test)
+endif()

app/main_1D.cpp

@@ -1,5 +1,5 @@
+#include "diffusion/BTCSBoundaryCondition.hpp"
 #include <diffusion/BTCSDiffusion.hpp>
-#include <diffusion/BoundaryCondition.hpp>
 
 #include <iomanip>
 #include <iostream> // for std
@@ -18,7 +18,8 @@ int main(int argc, char *argv[]) {
   // create input + diffusion coefficients for each grid cell
   std::vector<double> alpha(n, 1e-1);
   std::vector<double> field(n, 1e-6);
-  std::vector<boundary_condition> bc(n + 2, {0, 0});
+
+  BTCSBoundaryCondition bc;
 
   // create instance of diffusion module
   BTCSDiffusion diffu(dim);
@@ -26,7 +27,8 @@ int main(int argc, char *argv[]) {
   diffu.setXDimensions(1, n);
 
   // set the boundary condition for the left ghost cell to dirichlet
-  bc[0] = {Diffusion::BC_CONSTANT, 5e-6};
+  bc(BC_SIDE_LEFT) = {BC_TYPE_CONSTANT, 5e-6};
+  // bc[0] = {Diffusion::BC_CONSTANT, 5e-6};
   // diffu.setBoundaryCondition(1, 0, BTCSDiffusion::BC_CONSTANT,
   //                            5. * std::pow(10, -6));
 
@@ -38,7 +40,7 @@ int main(int argc, char *argv[]) {
   // loop 100 times
   // output is currently generated by the method itself
   for (int i = 0; i < 100; i++) {
-    diffu.simulate(field.data(), alpha.data(), bc.data());
+    diffu.simulate(field.data(), alpha.data(), bc);
 
     cout << "Iteration: " << i << "\n\n";
 

app/main_2D.cpp

@@ -1,5 +1,5 @@
+#include "diffusion/BTCSBoundaryCondition.hpp"
 #include <diffusion/BTCSDiffusion.hpp>
-#include <diffusion/BoundaryCondition.hpp>
 #include <iomanip>
 #include <iostream> // for std
 #include <vector>   // for vector
@@ -16,39 +16,42 @@ int main(int argc, char *argv[]) {
   int m = 5;
 
   // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(n * m, 1e-1);
+  std::vector<double> alpha(n * m, 1e-6);
-  std::vector<double> field(n * m, 1e-6);
+  std::vector<double> field(n * m, 0);
-  std::vector<boundary_condition> bc((n + 2) * (m + 2), {0, 0});
+  BTCSBoundaryCondition bc(n, m);
 
   // create instance of diffusion module
   BTCSDiffusion diffu(dim);
 
-  diffu.setXDimensions(1, n);
+  diffu.setXDimensions(n, n);
-  diffu.setYDimensions(1, m);
+  diffu.setYDimensions(m, m);
 
   // set inlet to higher concentration without setting bc
-  field[12] = 5e-3;
+  field[12] = 1;
 
   // set timestep for simulation to 1 second
-  diffu.setTimestep(1.);
+  diffu.setTimestep(1);
 
   cout << setprecision(12);
 
-  for (int t = 0; t < 10; t++) {
+  for (int t = 0; t < 1000; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc.data());
+    diffu.simulate(field.data(), alpha.data(), bc);
 
     cout << "Iteration: " << t << "\n\n";
 
+    double sum = 0;
+
     // iterate through rows
     for (int i = 0; i < m; i++) {
       // iterate through columns
       for (int j = 0; j < n; j++) {
         cout << left << std::setw(20) << field[i * n + j];
+        sum += field[i * n + j];
       }
       cout << "\n";
     }
 
-    cout << "\n" << endl;
+    cout << "sum: " << sum << "\n" << endl;
   }
 
   return 0;

app/main_2D_const.cpp

@@ -1,5 +1,5 @@
+#include "diffusion/BTCSBoundaryCondition.hpp"
 #include <diffusion/BTCSDiffusion.hpp>
-#include <diffusion/BoundaryCondition.hpp>
 #include <iomanip>
 #include <iostream> // for std
 #include <vector>   // for vector
@@ -18,7 +18,7 @@ int main(int argc, char *argv[]) {
   // create input + diffusion coefficients for each grid cell
   std::vector<double> alpha(n * m, 1e-1);
   std::vector<double> field(n * m, 1e-6);
-  std::vector<boundary_condition> bc((n + 2) * (m + 2), {0, 0});
+  BTCSBoundaryCondition bc(n, m);
 
   // create instance of diffusion module
   BTCSDiffusion diffu(dim);
@@ -26,16 +26,20 @@ int main(int argc, char *argv[]) {
   diffu.setXDimensions(1, n);
   diffu.setYDimensions(1, m);
 
-  for (int i = 1; i <= n; i++) {
+  boundary_condition input = {Diffusion::BC_TYPE_CONSTANT, 5e-6};
-    bc[(n + 2) * i] = {Diffusion::BC_CONSTANT, 5e-6};
+
-  }
+  bc.setSide(BC_SIDE_LEFT, input);
+
+  // for (int i = 1; i <= n; i++) {
+  //   bc[(n + 2) * i] = {Diffusion::BC_CONSTANT, 5e-6};
+  // }
   // set timestep for simulation to 1 second
   diffu.setTimestep(1.);
 
   cout << setprecision(12);
 
   for (int t = 0; t < 10; t++) {
-    diffu.simulate(field.data(), alpha.data(), bc.data());
+    diffu.simulate(field.data(), alpha.data(), bc);
 
     cout << "Iteration: " << t << "\n\n";
 

app/main_2D_mdl.cpp

@@ -1,5 +1,5 @@
+#include "diffusion/BTCSBoundaryCondition.hpp"
 #include <diffusion/BTCSDiffusion.hpp>
-#include <diffusion/BoundaryCondition.hpp>
 #include <iomanip>
 #include <iostream> // for std
 #include <vector>   // for vector
@@ -16,9 +16,9 @@ int main(int argc, char *argv[]) {
   int m = 501;
 
   // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(n * m, 1e-1);
+  std::vector<double> alpha(n * m, 1e-3);
   std::vector<double> field(n * m, 0.);
-  std::vector<boundary_condition> bc((n + 2) * (m + 2), {0, 0});
+  BTCSBoundaryCondition bc(n, m);
 
   field[125500] = 1;
 
@@ -51,7 +51,7 @@ int main(int argc, char *argv[]) {
 
   // Now we simulate and output 8 steps à 1 sec
   for (int t = 1; t < 6; t++) {
-    double time = diffu.simulate(field.data(), alpha.data(), bc.data());
+    double time = diffu.simulate(field.data(), alpha.data(), bc);
 
     cerr << "time elapsed: " << time << " seconds" << endl;
 

include/diffusion/BTCSBoundaryCondition.hpp

@@ -0,0 +1,186 @@
+#ifndef BOUNDARYCONDITION_H_
+#define BOUNDARYCONDITION_H_
+
+#include <array>
+#include <bits/stdint-uintn.h>
+#include <stdexcept>
+#include <stdint.h>
+#include <vector>
+
+typedef uint8_t bctype;
+
+namespace Diffusion {
+
+/**
+ * Defines a closed/Neumann boundary condition.
+ */
+static const bctype BC_TYPE_CLOSED = 0;
+
+/**
+ * Defines a flux/Cauchy boundary condition.
+ */
+static const bctype BC_TYPE_FLUX = 1;
+
+/**
+ * Defines a constant/Dirichlet boundary condition.
+ */
+static const bctype BC_TYPE_CONSTANT = 2;
+
+/**
+ * Defines boundary conditions for the left side of the grid.
+ */
+static const uint8_t BC_SIDE_LEFT = 0;
+
+/**
+ * Defines boundary conditions for the right side of the grid.
+ */
+static const uint8_t BC_SIDE_RIGHT = 1;
+
+/**
+ * Defines boundary conditions for the top of the grid.
+ */
+static const uint8_t BC_SIDE_TOP = 2;
+
+/**
+ * Defines boundary conditions for the bottom of the grid.
+ */
+static const uint8_t BC_SIDE_BOTTOM = 3;
+
+/**
+ * Defines the boundary condition type and according value.
+ */
+typedef struct boundary_condition_s {
+  bctype type;  /**< Type of the boundary condition */
+  double value; /**< Value of the boundary condition. Either a concrete value of
+                   concentration for BC_TYPE_CONSTANT or gradient when type is
+                   BC_TYPE_FLUX. Unused if BC_TYPE_CLOSED.*/
+} boundary_condition;
+
+/**
+ * Internal use only.
+ */
+typedef std::array<boundary_condition, 2> bc_tuple;
+
+/**
+ * Class to define the boundary condition of a grid.
+ */
+class BTCSBoundaryCondition {
+public:
+  /**
+   *  Creates a new instance with two elements. Used when defining boundary
+   *  conditions of 1D grids.
+   */
+  BTCSBoundaryCondition();
+
+  /**
+   * Creates a new instance with 4 * max(x,y) elements. Used to describe the
+   * boundary conditions for 2D grids. NOTE: On non-squared grids there are more
+   * elements than needed and no exception is thrown if some index exceeding
+   * grid limits.
+   *
+   * \param x Number of grid cells in x-direction
+   * \param y Number of grid cells in y-direction
+   *
+   */
+  BTCSBoundaryCondition(int x, int y);
+
+  /**
+   * Sets the boundary condition for a specific side of the grid.
+   *
+   * \param side Side for which the given boundary condition should be set.
+   * \param input_bc Instance of struct boundary_condition with desired boundary
+   * condition.
+   *
+   * \throws std::invalid_argument Indicates wrong dimensions of the grid.
+   * \throws std::out_of_range Indicates a out of range value for side.
+   */
+  void setSide(uint8_t side, boundary_condition &input_bc);
+
+  /**
+   * Get both boundary conditions of a given row (left and right).
+   *
+   * \param i Index of row
+   *
+   * \returns Left and right boundary values as an array (defined as data
+   * type bc_tuple).
+   */
+  auto row(uint32_t i) const -> bc_tuple;
+
+  /**
+   * Get both boundary conditions of a given column (top and bottom).
+   *
+   * \param i Index of column
+   *
+   * \returns Top and bottom boundary values as an array (defined as data
+   * type bc_tuple).
+   */
+  auto col(uint32_t i) const -> bc_tuple;
+
+  /**
+   * Create an instance of boundary_condition data type. Can be seen as a helper
+   * function.
+   *
+   * \param type Type of the boundary condition.
+   * \param value According value of condition.
+   *
+   * \returns Instance of boundary_condition
+   */
+  static boundary_condition returnBoundaryCondition(bctype type, double value) {
+    return {type, value};
+  }
+
+private:
+  std::vector<boundary_condition> bc_internal;
+
+  uint8_t dim;
+  uint32_t maxsize;
+
+public:
+  boundary_condition operator()(uint8_t side) const {
+    if (dim != 1) {
+      throw std::invalid_argument(
+          "Only 1D grid support 1 parameter in operator");
+    }
+    if (side > 1) {
+      throw std::out_of_range("1D index out of range");
+    }
+    return bc_internal[side];
+  }
+
+  boundary_condition operator()(uint8_t side, uint32_t i) const {
+    if (dim != 2) {
+      throw std::invalid_argument(
+          "Only 2D grids support 2 parameters in operator");
+    }
+    if (side > 3) {
+      throw std::out_of_range("2D index out of range");
+    }
+    return bc_internal[side * maxsize + i];
+  }
+
+  boundary_condition &operator()(uint8_t side) {
+    if (dim != 1) {
+      throw std::invalid_argument("Explicit setting of bc value with 2 "
+                                  "parameters is only supported for 2D grids");
+    }
+    if (side > 1) {
+      throw std::out_of_range("2D index out of range");
+    }
+    return bc_internal[side];
+  }
+
+  boundary_condition &operator()(uint8_t side, uint32_t i) {
+    if (dim != 2) {
+      throw std::invalid_argument("Explicit setting of bc value with 2 "
+                                  "parameters is only supported for 2D grids");
+    }
+    if (side > 3) {
+      throw std::out_of_range("2D index out of range");
+    }
+    return bc_internal[side * maxsize + i];
+  }
+};
+
+} // namespace Diffusion
+
+#endif // BOUNDARYCONDITION_H_

include/diffusion/BTCSDiffusion.hpp

@@ -1,7 +1,7 @@
 #ifndef BTCSDIFFUSION_H_
 #define BTCSDIFFUSION_H_
 
-#include "BoundaryCondition.hpp"
+#include "diffusion/BTCSBoundaryCondition.hpp"
 
 #include <Eigen/Sparse>
 #include <Eigen/src/Core/Map.h>
@@ -93,11 +93,11 @@ public:
    * state of each grid cell.
    * \param alpha Pointer to memory area of diffusion coefficients for each grid
    * element.
-   * \param bc Pointer to memory setting boundary conidition of each grid cell.
+   * \param bc Instance of boundary condition class with.
    *
    * \return Time in seconds [s] used to simulate one iteration.
    */
-  auto simulate(double *c, double *alpha, Diffusion::boundary_condition *bc)
+  auto simulate(double *c, double *alpha, const BTCSBoundaryCondition &bc)
       -> double;
 
   /*!
@@ -112,38 +112,31 @@ private:
       DMatrixRowMajor;
   typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>
       DVectorRowMajor;
-  typedef Eigen::Matrix<Diffusion::boundary_condition, Eigen::Dynamic,
-                        Eigen::Dynamic, Eigen::RowMajor>
-      BCMatrixRowMajor;
-  typedef Eigen::Matrix<Diffusion::boundary_condition, 1, Eigen::Dynamic,
-                        Eigen::RowMajor>
-      BCVectorRowMajor;
 
-  void simulate_base(DVectorRowMajor &c, const BCVectorRowMajor &bc,
+  static void simulate_base(DVectorRowMajor &c, const bc_tuple &bc,
                      const DVectorRowMajor &alpha, double dx, double time_step,
                      int size, const DVectorRowMajor &d_ortho);
 
   void simulate1D(Eigen::Map<DVectorRowMajor> &c,
                   Eigen::Map<const DVectorRowMajor> &alpha,
-                  Eigen::Map<const BCVectorRowMajor> &bc);
+                  const BTCSBoundaryCondition &bc);
 
   void simulate2D(Eigen::Map<DMatrixRowMajor> &c,
                   Eigen::Map<const DMatrixRowMajor> &alpha,
-                  Eigen::Map<const BCMatrixRowMajor> &bc);
+                  const BTCSBoundaryCondition &bc);
 
-  auto calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
+  static auto calc_d_ortho(const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
-                 const BCMatrixRowMajor &bc, double time_step, double dx)
+                    const BTCSBoundaryCondition &bc, bool transposed,
-      -> DMatrixRowMajor;
+                    double time_step, double dx) -> DMatrixRowMajor;
 
   static void fillMatrixFromRow(Eigen::SparseMatrix<double> &A_matrix,
-                                const DVectorRowMajor &alpha,
+                                const DVectorRowMajor &alpha, int size,
-                                const BCVectorRowMajor &bc, int size, double dx,
+                                double dx, double time_step);
-                                double time_step);
 
   static void fillVectorFromRow(Eigen::VectorXd &b_vector,
                                 const DVectorRowMajor &c,
                                 const DVectorRowMajor &alpha,
-                                const BCVectorRowMajor &bc,
+                                const bc_tuple &bc,
                                 const DVectorRowMajor &d_ortho, int size,
                                 double dx, double time_step);
   void simulate3D(std::vector<double> &c);

include/diffusion/BoundaryCondition.hpp

@@ -1,33 +0,0 @@
-#ifndef BOUNDARYCONDITION_H_
-#define BOUNDARYCONDITION_H_
-
-namespace Diffusion {
-
-/*!
- * Defines both types of boundary condition as a datatype.
- */
-typedef int bctype;
-
-/*!
- * Defines a closed/Neumann boundary condition.
- */
-static const bctype BC_CLOSED = 0;
-
-/*!
- * Defines a flux/Cauchy boundary condition.
- */
-static const bctype BC_FLUX = 1;
-
-/*!
- * Defines a constant/Dirichlet boundary condition.
- */
-static const bctype BC_CONSTANT = 2;
-
-typedef struct boundary_condition {
-  bctype type;
-  double value;
-} boundary_condition;
-
-} // namespace Diffusion
-
-#endif // BOUNDARYCONDITION_H_

src/BTCSBoundaryCondition.cpp

@@ -0,0 +1,56 @@
+#include <bits/stdint-uintn.h>
+#include <diffusion/BTCSBoundaryCondition.hpp>
+#include <stdexcept>
+
+constexpr uint8_t DIM_1D = 2;
+constexpr uint8_t DIM_2D = 4;
+
+Diffusion::BTCSBoundaryCondition::BTCSBoundaryCondition() {
+  this->bc_internal.resize(DIM_1D, {0, 0});
+  this->dim = 1;
+  // this value is actually unused
+  this->maxsize = 1;
+}
+
+Diffusion::BTCSBoundaryCondition::BTCSBoundaryCondition(int x, int y) {
+  this->maxsize = (x >= y ? x : y);
+  this->bc_internal.resize(DIM_2D * maxsize, {0, 0});
+  this->dim = 2;
+}
+
+void Diffusion::BTCSBoundaryCondition::setSide(
+    uint8_t side, Diffusion::boundary_condition &input_bc) {
+  if (this->dim == 1) {
+    throw std::invalid_argument("setSide requires at least a 2D grid");
+  }
+  if (side > 3) {
+    throw std::out_of_range("Invalid range for 2D grid");
+  }
+
+  for (int i = 0; i < maxsize; i++) {
+    this->bc_internal[side * maxsize + i] = input_bc;
+  }
+}
+
+auto Diffusion::BTCSBoundaryCondition::col(uint32_t i) const
+    -> Diffusion::bc_tuple {
+  if (this->dim == 1) {
+    throw std::invalid_argument("Access of column requires at least 2D grid");
+  }
+  if (i >= this->maxsize) {
+    throw std::out_of_range("Index out of range");
+  }
+
+  return {this->bc_internal[BC_SIDE_TOP * this->maxsize + i],
+          this->bc_internal[BC_SIDE_BOTTOM * this->maxsize + i]};
+}
+
+auto Diffusion::BTCSBoundaryCondition::row(uint32_t i) const
+    -> Diffusion::bc_tuple {
+  if (i >= this->maxsize) {
+    throw std::out_of_range("Index out of range");
+  }
+
+  return {this->bc_internal[BC_SIDE_LEFT * this->maxsize + i],
+          this->bc_internal[BC_SIDE_RIGHT * this->maxsize + i]};
+}

src/BTCSDiffusion.cpp

@@ -1,5 +1,5 @@
 #include "diffusion/BTCSDiffusion.hpp"
-#include "diffusion/BoundaryCondition.hpp"
+#include "diffusion/BTCSBoundaryCondition.hpp"
 
 #include <Eigen/SparseLU>
 
@@ -9,6 +9,7 @@
 #include <Eigen/src/SparseCore/SparseMatrixBase.h>
 #include <algorithm>
 #include <array>
+#include <bits/stdint-uintn.h>
 #include <cassert>
 #include <chrono>
 #include <cstddef>
@@ -25,7 +26,7 @@
 #define omp_get_thread_num() 0
 #endif
 
-#define DOUBLE_MACHINE_EPSILON 1.93e-34
+constexpr double DOUBLE_MACHINE_EPSILON = 1.93e-34;
 
 constexpr int BTCS_MAX_DEP_PER_CELL = 3;
 constexpr int BTCS_2D_DT_SIZE = 2;
@@ -87,7 +88,7 @@ void Diffusion::BTCSDiffusion::updateInternals() {
   }
 }
 void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
-                                             const BCVectorRowMajor &bc,
+                                             const Diffusion::bc_tuple &bc,
                                              const DVectorRowMajor &alpha,
                                              double dx, double time_step,
                                              int size,
@@ -103,7 +104,7 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
   b_vector.resize(size + 2);
   x_vector.resize(size + 2);
 
-  fillMatrixFromRow(A_matrix, alpha, bc, size, dx, time_step);
+  fillMatrixFromRow(A_matrix, alpha, size, dx, time_step);
   fillVectorFromRow(b_vector, c, alpha, bc, d_ortho, size, dx, time_step);
 
   // start to solve
@@ -120,7 +121,7 @@ void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
 
 void Diffusion::BTCSDiffusion::simulate1D(
     Eigen::Map<DVectorRowMajor> &c, Eigen::Map<const DVectorRowMajor> &alpha,
-    Eigen::Map<const BCVectorRowMajor> &bc) {
+    const Diffusion::BTCSBoundaryCondition &bc) {
 
   int size = this->grid_cells[0];
   double dx = this->deltas[0];
@@ -128,7 +129,7 @@ void Diffusion::BTCSDiffusion::simulate1D(
 
   DVectorRowMajor input_field = c.row(0);
 
-  simulate_base(input_field, bc, alpha, dx, time_step, size,
+  simulate_base(input_field, bc.row(0), alpha, dx, time_step, size,
                 Eigen::VectorXd::Constant(size, 0));
 
   c.row(0) << input_field;
@@ -136,7 +137,7 @@ void Diffusion::BTCSDiffusion::simulate1D(
 
 void Diffusion::BTCSDiffusion::simulate2D(
     Eigen::Map<DMatrixRowMajor> &c, Eigen::Map<const DMatrixRowMajor> &alpha,
-    Eigen::Map<const BCMatrixRowMajor> &bc) {
+    const Diffusion::BTCSBoundaryCondition &bc) {
 
   int n_rows = this->grid_cells[1];
   int n_cols = this->grid_cells[0];
@@ -144,35 +145,37 @@ void Diffusion::BTCSDiffusion::simulate2D(
 
   double local_dt = this->time_step / BTCS_2D_DT_SIZE;
 
-  DMatrixRowMajor d_ortho = calc_d_ortho(c, alpha, bc, local_dt, dx);
+  DMatrixRowMajor d_ortho = calc_d_ortho(c, alpha, bc, false, local_dt, dx);
 
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_rows; i++) {
     DVectorRowMajor input_field = c.row(i);
-    simulate_base(input_field, bc.row(i + 1), alpha.row(i), dx, local_dt,
+    simulate_base(input_field, bc.row(i), alpha.row(i), dx, local_dt, n_cols,
-                  n_cols, d_ortho.row(i));
+                  d_ortho.row(i));
     c.row(i) << input_field;
   }
 
   dx = this->deltas[1];
 
-  d_ortho = calc_d_ortho(c.transpose(), alpha.transpose(), bc.transpose(),
+  d_ortho =
-                         local_dt, dx);
+      calc_d_ortho(c.transpose(), alpha.transpose(), bc, true, local_dt, dx);
 
 #pragma omp parallel for schedule(dynamic)
   for (int i = 0; i < n_cols; i++) {
     DVectorRowMajor input_field = c.col(i);
-    simulate_base(input_field, bc.col(i + 1), alpha.col(i), dx, local_dt,
+    simulate_base(input_field, bc.col(i), alpha.col(i), dx, local_dt, n_rows,
-                  n_rows, d_ortho.row(i));
+                  d_ortho.row(i));
     c.col(i) << input_field.transpose();
   }
 }
 
-auto Diffusion::BTCSDiffusion::calc_d_ortho(const DMatrixRowMajor &c,
+auto Diffusion::BTCSDiffusion::calc_d_ortho(
-                                            const DMatrixRowMajor &alpha,
+    const DMatrixRowMajor &c, const DMatrixRowMajor &alpha,
-                                            const BCMatrixRowMajor &bc,
+    const Diffusion::BTCSBoundaryCondition &bc, bool transposed,
-                                            double time_step, double dx)
+    double time_step, double dx) -> DMatrixRowMajor {
-    -> DMatrixRowMajor {
+
+  uint8_t upper = (transposed ? Diffusion::BC_SIDE_LEFT : Diffusion::BC_SIDE_TOP);
+  uint8_t lower = (transposed ? Diffusion::BC_SIDE_RIGHT : Diffusion::BC_SIDE_BOTTOM);
 
   int n_rows = c.rows();
   int n_cols = c.cols();
@@ -183,10 +186,10 @@ auto Diffusion::BTCSDiffusion::calc_d_ortho(const DMatrixRowMajor &c,
 
   // first, iterate over first row
   for (int j = 0; j < n_cols; j++) {
-    boundary_condition tmp_bc = bc(0, j + 1);
+    boundary_condition tmp_bc = bc(upper, j);
     double sy = (time_step * alpha(0, j)) / (dx * dx);
 
-    y_values[0] = (tmp_bc.type == Diffusion::BC_CONSTANT
+    y_values[0] = (tmp_bc.type == Diffusion::BC_TYPE_CONSTANT
                        ? tmp_bc.value
                        : getBCFromFlux(tmp_bc, c(0, j), alpha(0, j)));
     y_values[1] = c(0, j);
@@ -213,12 +216,12 @@ auto Diffusion::BTCSDiffusion::calc_d_ortho(const DMatrixRowMajor &c,
 
   // and finally over last row
   for (int j = 0; j < n_cols; j++) {
-    boundary_condition tmp_bc = bc(end + 1, j + 1);
+    boundary_condition tmp_bc = bc(lower, j);
     double sy = (time_step * alpha(end, j)) / (dx * dx);
 
     y_values[0] = c(end - 1, j);
     y_values[1] = c(end, j);
-    y_values[2] = (tmp_bc.type == Diffusion::BC_CONSTANT
+    y_values[2] = (tmp_bc.type == Diffusion::BC_TYPE_CONSTANT
                        ? tmp_bc.value
                        : getBCFromFlux(tmp_bc, c(end, j), alpha(end, j)));
 
@@ -230,73 +233,49 @@ auto Diffusion::BTCSDiffusion::calc_d_ortho(const DMatrixRowMajor &c,
 
 void Diffusion::BTCSDiffusion::fillMatrixFromRow(
     Eigen::SparseMatrix<double> &A_matrix, const DVectorRowMajor &alpha,
-    const BCVectorRowMajor &bc, int size, double dx, double time_step) {
+    int size, double dx, double time_step) {
 
-  Diffusion::boundary_condition left = bc[1];
+  double sx = 0;
-  Diffusion::boundary_condition right = bc[size];
-
-  bool left_constant = (left.type == Diffusion::BC_CONSTANT);
-  bool right_constant = (right.type == Diffusion::BC_CONSTANT);
 
   int A_size = A_matrix.cols();
 
   A_matrix.insert(0, 0) = 1;
 
-  if (left_constant) {
+  sx = (alpha[0] * time_step) / (dx * dx);
-    A_matrix.insert(1, 1) = 1;
+  A_matrix.insert(1, 1) = -1. - 3. * sx;
-  } else {
+  A_matrix.insert(1, 0) = 2. * sx;
-    double sx = (alpha[0] * time_step) / (dx * dx);
+  A_matrix.insert(1, 2) = sx;
-    A_matrix.insert(1, 1) = -1. - 3. * sx;
-    A_matrix.insert(1, 0) = 2. * sx;
-    A_matrix.insert(1, 2) = sx;
-  }
 
   for (int j = 2, k = j - 1; k < size - 1; j++, k++) {
-    double sx = (alpha[k] * time_step) / (dx * dx);
+    sx = (alpha[k] * time_step) / (dx * dx);
-
-    if (bc[k + 1].type == Diffusion::BC_CONSTANT) {
-      A_matrix.insert(j, j) = 1;
-      continue;
-    }
 
     A_matrix.insert(j, j) = -1. - 2. * sx;
     A_matrix.insert(j, (j - 1)) = sx;
     A_matrix.insert(j, (j + 1)) = sx;
   }
 
-  if (right_constant) {
+  sx = (alpha[size - 1] * time_step) / (dx * dx);
-    A_matrix.insert(A_size - 2, A_size - 2) = 1;
+  A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
-  } else {
+  A_matrix.insert(A_size - 2, A_size - 3) = sx;
-    double sx = (alpha[size - 1] * time_step) / (dx * dx);
+  A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
-    A_matrix.insert(A_size - 2, A_size - 2) = -1. - 3. * sx;
-    A_matrix.insert(A_size - 2, A_size - 3) = sx;
-    A_matrix.insert(A_size - 2, A_size - 1) = 2. * sx;
-  }
 
   A_matrix.insert(A_size - 1, A_size - 1) = 1;
 }
 
 void Diffusion::BTCSDiffusion::fillVectorFromRow(
     Eigen::VectorXd &b_vector, const DVectorRowMajor &c,
-    const DVectorRowMajor &alpha, const BCVectorRowMajor &bc,
+    const DVectorRowMajor &alpha, const bc_tuple &bc,
     const DVectorRowMajor &d_ortho, int size, double dx, double time_step) {
 
   Diffusion::boundary_condition left = bc[0];
-  Diffusion::boundary_condition right = bc[size + 1];
+  Diffusion::boundary_condition right = bc[1];
 
-  bool left_constant = (left.type == Diffusion::BC_CONSTANT);
+  bool left_constant = (left.type == Diffusion::BC_TYPE_CONSTANT);
-  bool right_constant = (right.type == Diffusion::BC_CONSTANT);
+  bool right_constant = (right.type == Diffusion::BC_TYPE_CONSTANT);
 
   int b_size = b_vector.size();
 
   for (int j = 0; j < size; j++) {
-    boundary_condition tmp_bc = bc[j + 1];
-
-    if (tmp_bc.type == Diffusion::BC_CONSTANT) {
-      b_vector[j + 1] = tmp_bc.value;
-      continue;
-    }
-
     b_vector[j + 1] = -c[j] + d_ortho[j];
   }
 
@@ -314,8 +293,8 @@ void Diffusion::BTCSDiffusion::setTimestep(double time_step) {
   this->time_step = time_step;
 }
 
-auto Diffusion::BTCSDiffusion::simulate(double *c, double *alpha,
+auto Diffusion::BTCSDiffusion::simulate(
-                                        Diffusion::boundary_condition *bc)
+    double *c, double *alpha, const Diffusion::BTCSBoundaryCondition &bc)
     -> double {
 
   std::chrono::high_resolution_clock::time_point start =
@@ -324,9 +303,8 @@ auto Diffusion::BTCSDiffusion::simulate(double *c, double *alpha,
   if (this->grid_dim == 1) {
     Eigen::Map<DVectorRowMajor> c_in(c, this->grid_cells[0]);
     Eigen::Map<const DVectorRowMajor> alpha_in(alpha, this->grid_cells[0]);
-    Eigen::Map<const BCVectorRowMajor> bc_in(bc, this->grid_cells[0] + 2);
 
-    simulate1D(c_in, alpha_in, bc_in);
+    simulate1D(c_in, alpha_in, bc);
   }
   if (this->grid_dim == 2) {
     Eigen::Map<DMatrixRowMajor> c_in(c, this->grid_cells[1],
@@ -335,10 +313,7 @@ auto Diffusion::BTCSDiffusion::simulate(double *c, double *alpha,
     Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, this->grid_cells[1],
                                                this->grid_cells[0]);
 
-    Eigen::Map<const BCMatrixRowMajor> bc_in(bc, this->grid_cells[1] + 2,
+    simulate2D(c_in, alpha_in, bc);
-                                             this->grid_cells[0] + 2);
-
-    simulate2D(c_in, alpha_in, bc_in);
   }
 
   std::chrono::high_resolution_clock::time_point end =
@@ -357,9 +332,9 @@ inline auto Diffusion::BTCSDiffusion::getBCFromFlux(boundary_condition bc,
 
   double val = 0;
 
-  if (bc.type == Diffusion::BC_CLOSED) {
+  if (bc.type == Diffusion::BC_TYPE_CLOSED) {
     val = neighbor_c;
-  } else if (bc.type == Diffusion::BC_FLUX) {
+  } else if (bc.type == Diffusion::BC_TYPE_FLUX) {
     // TODO
     //  val = bc[index].value;
   } else {

src/CMakeLists.txt

@@ -1,7 +1,4 @@
-set(HEADER_LIST "${BTCSDiffusion_SOURCE_DIR}/include/diffusion/BTCSDiffusion.hpp"
+add_library(BTCSDiffusion STATIC BTCSDiffusion.cpp BTCSBoundaryCondition.cpp)
-  "${BTCSDiffusion_SOURCE_DIR}/include/diffusion/BoundaryCondition.hpp")
-
-add_library(BTCSDiffusion STATIC BTCSDiffusion.cpp ${HEADER_LIST})
 
 target_link_libraries(BTCSDiffusion Eigen3::Eigen)
 

test/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_library(doctest INTERFACE)
+target_include_directories(doctest INTERFACE doctest)
+
+add_executable(test setup.cpp testBoundaryCondition.cpp testDiffusion.cpp)
+target_link_libraries(test doctest BTCSDiffusion)

test/setup.cpp

@@ -0,0 +1,2 @@
+#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
+#include "doctest.h"

test/testBoundaryCondition.cpp

@@ -0,0 +1,83 @@
+#include <diffusion/BTCSBoundaryCondition.hpp>
+#include <doctest.h>
+
+using namespace Diffusion;
+
+#define BC_CONST_VALUE 1e-5
+
+TEST_CASE("1D Boundary Condition") {
+
+  BTCSBoundaryCondition bc;
+  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
+
+  SUBCASE("valid get") { CHECK_EQ(bc(BC_SIDE_LEFT).value, 0); }
+
+  SUBCASE("invalid get") {
+    CHECK_THROWS(bc(BC_SIDE_TOP));
+    CHECK_THROWS(bc(BC_SIDE_LEFT, 1));
+  }
+
+  SUBCASE("valid set") {
+    CHECK_NOTHROW(bc(BC_SIDE_LEFT) = bc_set);
+    CHECK_EQ(bc(BC_SIDE_LEFT).value, bc_set.value);
+    CHECK_EQ(bc(BC_SIDE_LEFT).type, bc_set.type);
+  }
+
+  SUBCASE("invalid set") {
+    CHECK_THROWS(bc(BC_SIDE_TOP) = bc_set);
+    CHECK_THROWS(bc(BC_SIDE_LEFT, 0) = bc_set);
+  }
+
+  SUBCASE("valid row getter") {
+    bc(BC_SIDE_LEFT) = bc_set;
+    bc_tuple tup = bc.row(0);
+
+    CHECK_EQ(tup[0].value, BC_CONST_VALUE);
+    CHECK_EQ(tup[1].value, 0);
+  }
+
+  SUBCASE("invalid row and col getter") {
+    CHECK_THROWS(bc.row(1));
+    CHECK_THROWS(bc.col(0));
+  }
+}
+
+TEST_CASE("2D Boundary Condition") {
+
+  BTCSBoundaryCondition bc(5, 5);
+  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
+
+  SUBCASE("valid get") { CHECK_EQ(bc(BC_SIDE_LEFT, 0).value, 0); }
+
+  SUBCASE("invalid get") {
+    CHECK_THROWS(bc(4, 0));
+    CHECK_THROWS(bc(BC_SIDE_LEFT));
+  }
+
+  SUBCASE("valid set") {
+    CHECK_NOTHROW(bc(BC_SIDE_LEFT, 0) = bc_set);
+    CHECK_EQ(bc(BC_SIDE_LEFT, 0).value, bc_set.value);
+    CHECK_EQ(bc(BC_SIDE_LEFT, 0).type, bc_set.type);
+  }
+
+  SUBCASE("invalid set") {
+    CHECK_THROWS(bc(BC_SIDE_LEFT) = bc_set);
+    CHECK_THROWS(bc(4, 0) = bc_set);
+  }
+
+  SUBCASE("call of setSide") {
+      CHECK_NOTHROW(bc.setSide(BC_SIDE_BOTTOM, bc_set));
+      CHECK_EQ(bc(BC_SIDE_BOTTOM, 1).value, bc_set.value);
+      CHECK_EQ(bc(BC_SIDE_BOTTOM, 1).type, bc_set.type);
+  }
+}
+
+TEST_CASE("Boundary Condition helpers") {
+  boundary_condition bc_set = {BC_TYPE_CONSTANT, BC_CONST_VALUE};
+
+  SUBCASE("return boundary condition skeleton") {
+    boundary_condition bc_test = BTCSBoundaryCondition::returnBoundaryCondition(bc_set.type, bc_set.value);
+    CHECK_EQ(bc_test.value, bc_set.value);
+    CHECK_EQ(bc_test.type, bc_set.type);
+  }
+}

test/testDiffusion.cpp

@@ -0,0 +1,125 @@
+#include <bits/stdint-uintn.h>
+#include <diffusion/BTCSBoundaryCondition.hpp>
+#include <diffusion/BTCSDiffusion.hpp>
+#include <doctest.h>
+#include <vector>
+
+using namespace Diffusion;
+
+#define DIMENSION 2
+#define N 51
+#define M 51
+#define MID 1300
+
+static std::vector<double> alpha(N *M, 1e-3);
+
+static BTCSDiffusion setupDiffu(uint32_t n, uint32_t m) {
+  BTCSDiffusion diffu(DIMENSION);
+
+  diffu.setXDimensions(n, n);
+  diffu.setYDimensions(m, m);
+
+  diffu.setTimestep(1.);
+
+  return diffu;
+}
+
+TEST_CASE("closed boundaries - 1 concentration to 1 - rest 0") {
+  std::vector<double> field(N * M, 0);
+
+  field[MID] = 1;
+
+  BTCSDiffusion diffu = setupDiffu(N, M);
+  BTCSBoundaryCondition bc(N, M);
+
+  uint32_t iterations = 1000;
+  double sum = 0;
+
+  for (int t = 0; t < iterations; t++) {
+    diffu.simulate(field.data(), alpha.data(), bc);
+
+    if (t == iterations - 1) {
+      // iterate through rows
+      for (int i = 0; i < M; i++) {
+        // iterate through columns
+        for (int j = 0; j < N; j++) {
+          sum += field[i * N + j];
+        }
+      }
+    }
+  }
+  CAPTURE(sum);
+  // epsilon of 1e-8
+  CHECK(sum == doctest::Approx(1).epsilon(1e-6));
+}
+
+TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") {
+  std::vector<double> field(N * M, 0);
+
+  field[MID] = 1;
+
+  BTCSDiffusion diffu = setupDiffu(N, M);
+  BTCSBoundaryCondition bc(N, M);
+
+  boundary_condition input = {BC_TYPE_CONSTANT, 0};
+
+  bc.setSide(BC_SIDE_LEFT, input);
+  bc.setSide(BC_SIDE_RIGHT, input);
+  bc.setSide(BC_SIDE_TOP, input);
+  bc.setSide(BC_SIDE_BOTTOM, input);
+
+  uint32_t max_iterations = 20000;
+  bool reached = false;
+
+  int t = 0;
+
+  for (t = 0; t < max_iterations; t++) {
+    diffu.simulate(field.data(), alpha.data(), bc);
+
+    if (field[N * M - 1] > 1e-15) {
+      reached = true;
+      break;
+    }
+  }
+
+  if (!reached) {
+    CAPTURE(field[N * M - 1]);
+    FAIL_CHECK(
+        "Concentration did not reach boundaries after count of iterations: ",
+        t);
+  }
+}
+
+TEST_CASE(
+    "constant top and bottom (1 and 0) - left and right closed - 0 inlet") {
+  std::vector<double> field(N * M, 0);
+
+  BTCSDiffusion diffu = setupDiffu(N, M);
+  BTCSBoundaryCondition bc(N, M);
+
+  boundary_condition top =
+      BTCSBoundaryCondition::returnBoundaryCondition(BC_TYPE_CONSTANT, 1);
+  boundary_condition bottom =
+      BTCSBoundaryCondition::returnBoundaryCondition(BC_TYPE_CONSTANT, 0);
+
+  bc.setSide(BC_SIDE_TOP, top);
+  bc.setSide(BC_SIDE_BOTTOM, bottom);
+
+  uint32_t max_iterations = 100;
+
+  for (int t = 0; t < max_iterations; t++) {
+    diffu.simulate(field.data(), alpha.data(), bc);
+  }
+
+  for (int i = 0; i < N; i++) {
+    double above = field[i];
+    for (int j = 1; j < M; j++) {
+      double curr = field[j * N + i];
+      if (curr > above) {
+        CAPTURE(curr);
+        CAPTURE(above);
+        FAIL("Concentration below is greater than above @ cell ", j * N + i);
+      }
+    }
+  }
+}