Cleanup and refactoring of library.

- Update doxygen compatible comments for all public methods of BTCSDiffusion
- Remove commented code lines as we do noit need them anymore.
- Added comments at key points at source file of BTCSDiffusion.
- Refactor main.cpp to new function signature of =setBoundaryCondtion=

src/BTCSDiffusion.cpp

@@ -67,28 +67,6 @@ void BTCSDiffusion::updateInternals() {
 
   bc.resize(cells, {BTCSDiffusion::BC_CLOSED, 0});
 }
-// BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) {
-//   this->grid_dim = 1;
-//   this->dx = 1. / (x - 1);
-
-//   // per default use Neumann condition with gradient of 0 at the end of the
-//   grid this->bc.resize(2, std::tuple<bctype,
-//   double>(BTCSDiffusion::BC_CONSTANT, 0.));
-// }
-// BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) {
-
-//   // this->grid_dim = 2;
-
-//   // this->bc.reserve(x * 2 + y * 2);
-//   // // per default use Neumann condition with gradient of 0 at the end of
-//   the
-//   // grid std::fill(this->bc.begin(), this->bc.end(), -1);
-// }
-// BTCSDiffusion::BTCSDiffusion(int x, int y, int z) : n_x(x), n_y(y), n_z(z) {
-
-//   // this->grid_dim = 3;
-//   // TODO: reserve memory for boundary conditions
-// }
 
 void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
                                boundary_condition right,
@@ -97,11 +75,11 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
 
   bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT);
   bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT);
-  int loop_end = size + !right_is_constant;
-
-  // we need 2 more grid cells for ghost cells
-  // size = size + 2;
 
+  //The sizes for matrix and vectors of the equation system is defined by the
+  //actual size of the input vector and if the system is (partially) closed.
+  //Then we will need ghost nodes. So this variable will give the count of ghost
+  //nodes.
   int bc_offset = !left_is_constant + !right_is_constant;
   ;
 
@@ -113,11 +91,10 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
    * Begin to solve the equation system using LU solver of Eigen.
    *
    * But first fill the A matrix and b vector.
-   *
-   * At this point there is some debugging output in the code.
-   * TODO: remove output
    */
 
+  // Set boundary condition for ghost nodes (for closed or flux system) or outer
+  // inlet nodes (constant boundary condition)
   A_matrix.resize(size + bc_offset, size + bc_offset);
   A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3));
 
@@ -130,11 +107,13 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
       (right_is_constant ? right.value
                          : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
 
-  // A_matrix.insert(0, 0) = 1;
+  // Start filling the A matrix
-  // A_matrix.insert(size + 1, size + 1) = 1;
+  // =i= is used for equation system matrix and vector indexing
-
+  // and =j= for indexing of c,alpha and bc
   for (int i = 1, j = i + !(left_is_constant); i < size - right_is_constant;
        i++, j++) {
+
+    // if current grid cell is considered as constant boundary conditon
     if (bc[j].type == BTCSDiffusion::BC_CONSTANT) {
       A_matrix.insert(i, i) = 1;
       b_vector[i] = bc[j].value;
@@ -150,6 +129,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
     b_vector[i] = -c[j];
   }
 
+  // start to solve
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       solver;
   solver.analyzePattern(A_matrix);
@@ -158,7 +138,8 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
 
   x_vector = solver.solve(b_vector);
 
-  for (int i = 0; i < c.size(); i++) {
+  //fill solution back in place into =c= vector
+  for (int i = 0, j = i + !left_is_constant; i < c.size(); i++, j++) {
     c[i] = x_vector[i + !left_is_constant];
   }
 }
@@ -170,10 +151,6 @@ void BTCSDiffusion::setTimestep(double time_step) {
 void BTCSDiffusion::simulate(std::vector<double> &c,
                              const std::vector<double> &alpha) {
   if (this->grid_dim == 1) {
-    // double bc_left = getBCFromTuple(0, c[0], alpha[0]);
-    // double bc_right =
-    //     getBCFromTuple(1, c[c.size() - 1], alpha[alpha.size() - 1]);
-
     simulate1D(c, bc[0], bc[grid_cells[0] + 1], alpha, this->deltas[0],
                this->grid_cells[0]);
   }
@@ -197,11 +174,8 @@ inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc,
   return val;
 }
 
-void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) {
+void BTCSDiffusion::setBoundaryCondition(int index, bctype type, double value) {
 
   bc[index].type = type;
-  bc[index].value = val;
+  bc[index].value = value;
-
-  // std::get<0>(bc[index]) = type;
-  // std::get<1>(bc[index]) = val;
 }

src/BTCSDiffusion.hpp

@@ -34,76 +34,76 @@ public:
   static const int BC_FLUX;
 
   /*!
-   * A boundary condition consists of two features. A type and the according
+   * Creates a diffusion module.
-   * value. Here we can differentiate between:
    *
-   * - Neumann boundary conditon: type BC_NEUMANN with the value defining the
+   * @param dim Number of dimensions. Should not be greater than 3 and not less
-   * gradient
+   * than 1.
-   * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
-   * the boundary condition
-   */
-  typedef struct boundary_condition {
-    bctype type;
-    double value;
-  } boundary_condition;
-
-  /*!
-   * A boundary condition consists of two features. A type and the according
-   * value. Here we can differentiate between:
-   *
-   * - Neumann boundary conditon: type BC_NEUMANN with the value defining the
-   * gradient
-   * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
-   * the boundary condition
-   */
-  // typedef std::vector<std::tuple<bctype, double>> boundary_condition;
-
-  /*!
-   * Datatype to fill the sparse matrix which is used to solve the equation
-   * system.
-   */
-  typedef Eigen::Triplet<double> T;
-
-  /*!
-   * Create 1D-diffusion module.
-   *
-   * @param x Count of cells in x direction.
    */
   BTCSDiffusion(unsigned int dim);
 
+  /*!
+   * Define the grid in x direction.
+   *
+   * @param domain_size Size of the domain in x direction.
+   * @param n_grid_cells Number of grid cells in x direction the domain is
+   * divided to.
+   */
   void setXDimensions(unsigned int domain_size, unsigned int n_grid_cells);
+
+  /*!
+   * Define the grid in y direction.
+   *
+   * Throws an error if the module wasn't initialized at least as a 2D model.
+   *
+   * @param domain_size Size of the domain in y direction.
+   * @param n_grid_cells Number of grid cells in y direction the domain is
+   * divided to.
+   */
   void setYDimensions(unsigned int domain_size, unsigned int n_grid_cells);
+
+  /*!
+   * Define the grid in z direction.
+   *
+   * Throws an error if the module wasn't initialized at least as a 3D model.
+   *
+   * @param domain_size Size of the domain in z direction.
+   * @param n_grid_cells Number of grid cells in z direction the domain is
+   * divided to.
+   */
   void setZDimensions(unsigned int domain_size, unsigned int n_grid_cells);
 
+  /*!
+   * Returns the number of grid cells in x direction.
+   */
   unsigned int getXGridCellsN();
+  /*!
+   * Returns the number of grid cells in y direction.
+   */
   unsigned int getYGridCellsN();
+  /*!
+   * Returns the number of grid cells in z direction.
+   */
   unsigned int getZGridCellsN();
-  unsigned int getXDomainSize();
-  unsigned int getYDomainSize();
-  unsigned int getZDomainSize();
-  // /*!
-  //  * Currently not implemented: Create 2D-diffusion module.
-  //  *
-  //  * @param x Count of cells in x direction.
-  //  * @param y Count of cells in y direction.
-  //  */
-  // explicit BTCSDiffusion(int x, int y);
 
-  // /*!
+  /*!
-  //  * Currently not implemented: Create 3D-diffusion module.
+   * Returns the domain size in x direction.
-  //  *
+   */
-  //  * @param x Count of cells in x direction.
+  unsigned int getXDomainSize();
-  //  * @param y Count of cells in y direction.
+  /*!
-  //  * @param z Count of cells in z direction.
+   * Returns the domain size in y direction.
-  //  */
+   */
-  // explicit BTCSDiffusion(int x, int y, int z);
+  unsigned int getYDomainSize();
+  /*!
+   * Returns the domain size in z direction.
+   */
+  unsigned int getZDomainSize();
 
   /*!
    * With given ghost zones simulate diffusion. Only 1D allowed at this moment.
    *
    * @param c Vector describing the concentration of one solution of the grid as
-   * continious memory (Row-wise).
+   * continious memory (row major).
-   * @param alpha Vector of diffusioncoefficients for each grid element.
+   * @param alpha Vector of diffusion coefficients for each grid element.
    */
   void simulate(std::vector<double> &c, const std::vector<double> &alpha);
 
@@ -119,23 +119,29 @@ public:
    * index (exact order still to be determined), the type of the boundary
    * condition and the according value.
    *
-   * @param index Index of the boundary condition vector.
+   * @param index Index of the grid cell the boundary condition is applied to.
-   * @param val Value of the boundary condition (gradient for Neumann, exact
+   * @param type Type of the boundary condition. Must be constant, closed or
-   * value for Dirichlet).
+   * flux.
-   * @param Type of the grid cell.
+   * @param value For constant boundary conditions this value is set
+   * during solving. For flux value refers to a gradient of change for this grid
+   * cell. For closed this value has no effect since a gradient of 0 is used.
    */
-  void setBoundaryCondition(int index, double val, bctype type);
+  void setBoundaryCondition(int index, bctype type, double value);
 
 private:
+  typedef struct boundary_condition {
+    bctype type;
+    double value;
+  } boundary_condition;
+  typedef Eigen::Triplet<double> T;
+
   void simulate1D(std::vector<double> &c, boundary_condition left,
                   boundary_condition right, const std::vector<double> &alpha,
                   double dx, int size);
   void simulate2D(std::vector<double> &c);
   void simulate3D(std::vector<double> &c);
-
   inline double getBCFromFlux(boundary_condition bc, double nearest_value,
                               double neighbor_alpha);
-
   void updateInternals();
 
   std::vector<boundary_condition> bc;
@@ -145,8 +151,8 @@ private:
   Eigen::VectorXd x_vector;
 
   double time_step;
-
   int grid_dim;
+
   std::vector<unsigned int> grid_cells;
   std::vector<unsigned int> domain_size;
   std::vector<double> deltas;

src/main.cpp

@@ -1,5 +1,4 @@
 #include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
-#include <Eigen/src/Core/arch/SSE/PacketMath.h>
 #include <algorithm> // for copy, max
 #include <iomanip>
 #include <iostream> // for std
@@ -23,8 +22,8 @@ int main(int argc, char *argv[]) {
   diffu.setXDimensions(1, n);
 
   // set the boundary condition for the left ghost cell to dirichlet
-  diffu.setBoundaryCondition(0, 5. * std::pow(10, -6),
+  diffu.setBoundaryCondition(0, BTCSDiffusion::BC_CONSTANT,
-                             BTCSDiffusion::BC_CONSTANT);
+                             5. * std::pow(10, -6));
 
   // set timestep for simulation to 1 second
   diffu.setTimestep(1.);