Added comments

src/BTCSDiffusion.cpp

@@ -13,6 +13,7 @@
 #include <Eigen/src/SparseLU/SparseLU.h>
 #include <Eigen/src/SparseQR/SparseQR.h>
 
+#include <algorithm>
 #include <iomanip>
 #include <iostream>
 
@@ -21,18 +22,23 @@ const BCSide BTCSDiffusion::RIGHT = 1;
 
 BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) {
   this->grid_dim = 1;
-  this->bc.reserve(2);
+
+  // per default use Neumann condition with gradient of 0 at the end of the grid
+  this->bc.resize(2, -1);
 }
 BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_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)
     : dim_x(x), dim_y(y), dim_z(z) {
 
   this->grid_dim = 3;
-  //TODO: reserve memory for boundary conditions
+  // TODO: reserve memory for boundary conditions
 }
 
 void BTCSDiffusion::setBoundaryCondition(std::vector<double> input,
@@ -43,17 +49,27 @@ void BTCSDiffusion::setBoundaryCondition(std::vector<double> input,
 }
 void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
                              double timestep) {
+  // calculate dx
   double dx = 1. / this->dim_x;
 
+  // calculate size needed for A matrix and b,x vectors
   int size = this->dim_x + 2;
 
   Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
   Eigen::VectorXd x_out(size);
+
+  /*
+   * Initalization of matrix A
+   * This is done by triplets. See:
+   * https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
+   */
+
   std::vector<T> tripletList;
   tripletList.reserve(c.size() * 3 + bc.size());
 
   int A_line = 0;
 
+  // For all concentrations create one row in matrix A
   for (int i = 1; i < this->dim_x + 1; i++) {
     double sx = (alpha[i - 1] * timestep) / (dx * dx);
 
@@ -66,9 +82,14 @@ void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
     A_line++;
   }
 
+  // append left and right boundary conditions/ghost zones
   tripletList.push_back(T(A_line, 0, 1));
+
+  // if value is -1 apply Neumann condition with given gradient
+  // TODO: set specific gradient
   if (bc[0] == -1)
     b[A_line] = c[0];
+  // else apply given Dirichlet condition
   else
     b[A_line] = this->bc[0];
 
@@ -80,6 +101,13 @@ void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
   else
     b[A_line] = this->bc[1];
 
+  /*
+   * Begin to solve the equation system
+   *
+   * At this point there is some debugging output in the code.
+   * TODO: remove output
+   */
+
   Eigen::SparseMatrix<double> A(size, size);
   A.setFromTriplets(tripletList.begin(), tripletList.end());
 

src/BTCSDiffusion.hpp

@@ -1,23 +1,80 @@
 #ifndef BTCSDIFFUSION_H_
 #define BTCSDIFFUSION_H_
 
-#include <vector>
 #include <Eigen/Sparse>
+#include <vector>
 
+/*!
+ * Type defining the side of given boundary condition.
+ */
 typedef int BCSide;
+
+/*!
+ * Datatype to fill the sparse matrix which is used to solve the equation
+ * system.
+ */
 typedef Eigen::Triplet<double> T;
 
+/*!
+ * Class implementing a solution for a 1/2/3D diffusion equation using backward
+ * euler.
+ */
 class BTCSDiffusion {
 
 public:
+  /*!
+   * Set left boundary condition.
+   */
   static const BCSide LEFT;
+
+  /*!
+   * Set right boundary condition.
+   */
   static const BCSide RIGHT;
 
+  /*!
+   * Create 1D-diffusion module.
+   *
+   * @param x Count of cells in x direction.
+   */
   BTCSDiffusion(int x);
+
+  /*!
+   * Currently not implemented: Create 2D-diffusion module.
+   *
+   * @param x Count of cells in x direction.
+   * @param y Count of cells in y direction.
+   */
   BTCSDiffusion(int x, int y);
+
+  /*!
+   * Currently not implemented: Create 3D-diffusion module.
+   *
+   * @param x Count of cells in x direction.
+   * @param y Count of cells in y direction.
+   * @param z Count of cells in z direction.
+   */
   BTCSDiffusion(int x, int y, int z);
 
+  /*!
+   * Sets internal boundary condition at the end of the grid/ghost zones.
+   * Currently only implemented for 1D diffusion.
+   *
+   * @param input Vector containing all the values to initialize the ghost
+   * zones.
+   * @param side Sets the side of the boundary condition. See BCSide for more
+   * information.
+   */
   void setBoundaryCondition(std::vector<double> input, BCSide side);
+
+  /*!
+   * 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).
+   * @param alpha Vector of diffusioncoefficients for each grid element.
+   * @param timestep Time (in seconds ?) to simulate.
+   */
   void simulate(std::vector<double> &c, std::vector<double> &alpha,
                 double timestep);
 

src/main.cpp

@@ -19,7 +19,9 @@ int main(int argc, char *argv[]) {
   BTCSDiffusion diffu(x);
 
   diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT);
-  diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT);
+  // we don't need this since Neumann condition with gradient of 0 is set per
+  // default
+  // diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT);
 
   for (int i = 0; i < 100; i++) {
     diffu.simulate(input, alpha, 1.);