From f5c926e08d94733ef6343b7c9128238be46c2857 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 6 Dec 2021 13:48:27 +0100
Subject: [PATCH 01/23] rewrite some function signatures and scopes (NOT
 RUNNABLE)

---
 src/BTCSDiffusion.cpp | 23 +++++++++--------------
 src/BTCSDiffusion.hpp | 33 ++++++++-------------------------
 2 files changed, 17 insertions(+), 39 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index fcebc92..251eb59 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -16,37 +16,32 @@
 #include <algorithm>
 #include <iomanip>
 #include <iostream>
+#include <tuple>
 
-const BCSide BTCSDiffusion::LEFT = 0;
-const BCSide BTCSDiffusion::RIGHT = 1;
+const int BTCSDiffusion::BC_NEUMANN = 0;
+const int BTCSDiffusion::BC_DIRICHLET = 1;
 
 BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) {
   this->grid_dim = 1;
 
   // per default use Neumann condition with gradient of 0 at the end of the grid
-  this->bc.resize(2, -1);
+  this->bc.resize(2, std::tuple<int, double>(0,0.));
 }
 BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) {
 
-  this->grid_dim = 2;
+  // 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);
+  // 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;
+  // this->grid_dim = 3;
   // TODO: reserve memory for boundary conditions
 }
 
-void BTCSDiffusion::setBoundaryCondition(std::vector<double> input,
-                                         BCSide side) {
-  if (this->grid_dim == 1) {
-    bc[side] = input[0];
-  }
-}
 void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
                              double timestep) {
   // calculate dx
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index cadc1e7..8d0bc3e 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -2,19 +2,17 @@
 #define BTCSDIFFUSION_H_
 
 #include <Eigen/Sparse>
+#include <tuple>
 #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;
 
+typedef std::vector<std::tuple<int,double>> boundary_condition;
+
 /*!
  * Class implementing a solution for a 1/2/3D diffusion equation using backward
  * euler.
@@ -22,15 +20,9 @@ typedef Eigen::Triplet<double> T;
 class BTCSDiffusion {
 
 public:
-  /*!
-   * Set left boundary condition.
-   */
-  static const BCSide LEFT;
 
-  /*!
-   * Set right boundary condition.
-   */
-  static const BCSide RIGHT;
+    static const int BC_NEUMANN;
+    static const int BC_DIRICHLET;
 
   /*!
    * Create 1D-diffusion module.
@@ -56,17 +48,6 @@ public:
    */
   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.
    *
@@ -79,7 +60,9 @@ public:
                 double timestep);
 
 private:
-  std::vector<double> bc;
+
+  boundary_condition bc;
+
   int grid_dim;
   int dim_x;
   int dim_y;

From 8d272741011eacf5924751e50d2085df30f23560 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 6 Dec 2021 20:03:40 +0100
Subject: [PATCH 02/23] use private datatypes to adress solver matrix

---
 src/BTCSDiffusion.cpp | 143 ++++++++++++++++++++++++++++--------------
 src/BTCSDiffusion.hpp |  27 ++++++--
 src/main.cpp          |   7 ++-
 3 files changed, 121 insertions(+), 56 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 251eb59..802c9db 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -6,12 +6,6 @@
 #include <Eigen/src/Core/Matrix.h>
 #include <Eigen/src/Core/util/Constants.h>
 #include <Eigen/src/OrderingMethods/Ordering.h>
-#include <Eigen/src/SparseCholesky/SimplicialCholesky.h>
-#include <Eigen/src/SparseCore/SparseMap.h>
-#include <Eigen/src/SparseCore/SparseMatrix.h>
-#include <Eigen/src/SparseCore/SparseMatrixBase.h>
-#include <Eigen/src/SparseLU/SparseLU.h>
-#include <Eigen/src/SparseQR/SparseQR.h>
 
 #include <algorithm>
 #include <iomanip>
@@ -25,15 +19,15 @@ BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) {
   this->grid_dim = 1;
 
   // per default use Neumann condition with gradient of 0 at the end of the grid
-  this->bc.resize(2, std::tuple<int, double>(0,0.));
+  this->bc.resize(2, std::tuple<int, double>(0, 0.));
 }
 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);
+  // // 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) {
@@ -42,16 +36,20 @@ BTCSDiffusion::BTCSDiffusion(int x, int y, int z)
   // TODO: reserve memory for boundary conditions
 }
 
-void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
-                             double timestep) {
+void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
+                               double bc_right, std::vector<double> &alpha) {
   // calculate dx
   double dx = 1. / (this->dim_x - 1);
 
   // 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);
+  // set sizes of private and yet allocated vectors
+  b_vector.resize(size);
+  x_vector.resize(size);
+
+  // Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
+  // Eigen::VectorXd x_out(size);
 
   /*
    * Initalization of matrix A
@@ -59,42 +57,42 @@ void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
    * https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
    */
 
-  std::vector<T> tripletList;
-  tripletList.reserve(c.size() * 3 + bc.size());
+  // std::vector<T> tripletList;
+  // tripletList.reserve(c.size() * 3 + bc.size());
 
-  int A_line = 0;
+  // 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);
+  // // 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);
 
-    tripletList.push_back(T(A_line, i, (-1. - 2. * sx)));
+  //   tripletList.push_back(T(A_line, i, (-1. - 2. * sx)));
 
-    tripletList.push_back(T(A_line, i - 1, sx));
-    tripletList.push_back(T(A_line, i + 1, sx));
+  //   tripletList.push_back(T(A_line, i - 1, sx));
+  //   tripletList.push_back(T(A_line, i + 1, sx));
 
-    b[A_line] = -c[i - 1];
-    A_line++;
-  }
+  //   b[A_line] = -c[i - 1];
+  //   A_line++;
+  // }
 
-  // append left and right boundary conditions/ghost zones
-  tripletList.push_back(T(A_line, 0, 1));
+  // // 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];
+  // // 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];
 
-  A_line++;
-  tripletList.push_back(T(A_line, size - 1, 1));
-  // b[A_line] = bc[1];
-  if (bc[1] == -1)
-    b[A_line] = c[c.size() - 1];
-  else
-    b[A_line] = this->bc[1];
+  // A_line++;
+  // tripletList.push_back(T(A_line, size - 1, 1));
+  // // b[A_line] = bc[1];
+  // if (bc[1] == -1)
+  //   b[A_line] = c[c.size() - 1];
+  // else
+  //   b[A_line] = this->bc[1];
 
   /*
    * Begin to solve the equation system
@@ -103,22 +101,71 @@ void BTCSDiffusion::simulate(std::vector<double> &c, std::vector<double> &alpha,
    * TODO: remove output
    */
 
-  Eigen::SparseMatrix<double> A(size, size);
-  A.setFromTriplets(tripletList.begin(), tripletList.end());
+  A_matrix.resize(size, size);
+  A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
+
+  A_matrix.insert(0, 0) = bc_left;
+  A_matrix.insert(size - 1, size - 1) = bc_right;
+
+  for (int i = 1; i < this->dim_x + 1; i++) {
+    double sx = (alpha[i - 1] * time_step) / (dx * dx);
+
+    A_matrix.insert(i, i) = -1. - 2. * sx;
+    A_matrix.insert(i, i - 1) = sx;
+    A_matrix.insert(i, i + 1) = sx;
+
+    b_vector[i] = -c[i - 1];
+
+    // tripletList.push_back(T(A_line, i, (-1. - 2. * sx)));
+
+    // tripletList.push_back(T(A_line, i - 1, sx));
+    // tripletList.push_back(T(A_line, i + 1, sx));
+
+    // b[A_line] = -c[i - 1];
+    // A_line++;
+  }
+
+  // Eigen::SparseMatrix<double> A(size, size);
+  // A.setFromTriplets(tripletList.begin(), tripletList.end());
 
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       solver;
-  solver.analyzePattern(A);
+  solver.analyzePattern(A_matrix);
 
-  solver.factorize(A);
+  solver.factorize(A_matrix);
 
   std::cout << solver.lastErrorMessage() << std::endl;
 
-  x_out = solver.solve(b);
+  x_vector = solver.solve(b_vector);
 
-  std::cout << std::setprecision(10) << x_out << std::endl << std::endl;
+  std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
 
   for (int i = 0; i < c.size(); i++) {
-    c[i] = x_out[i + 1];
+    c[i] = x_vector[i + 1];
   }
 }
+
+void BTCSDiffusion::setTimestep(double time_step) {
+  this->time_step = time_step;
+}
+
+void BTCSDiffusion::simulate(std::vector<double> &c,
+                             std::vector<double> &alpha) {
+  if (this->grid_dim == 1) {
+    double bc_left = getBCFromTuple(0);
+    double bc_right = getBCFromTuple(1);
+
+    simulate1D(c, bc_left, bc_right, alpha);
+  }
+}
+
+double BTCSDiffusion::getBCFromTuple(int index) {
+  double val = std::get<1>(bc[index]);
+
+  return val;
+}
+
+void BTCSDiffusion::setBoundaryCondition(int index, double val, int type) {
+  std::get<0>(bc[index]) = val;
+  std::get<1>(bc[index]) = type;
+}
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 8d0bc3e..27b2aed 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -2,6 +2,7 @@
 #define BTCSDIFFUSION_H_
 
 #include <Eigen/Sparse>
+#include <Eigen/src/SparseCore/SparseMatrixBase.h>
 #include <tuple>
 #include <vector>
 
@@ -11,7 +12,7 @@
  */
 typedef Eigen::Triplet<double> T;
 
-typedef std::vector<std::tuple<int,double>> boundary_condition;
+typedef std::vector<std::tuple<int, double>> boundary_condition;
 
 /*!
  * Class implementing a solution for a 1/2/3D diffusion equation using backward
@@ -20,9 +21,8 @@ typedef std::vector<std::tuple<int,double>> boundary_condition;
 class BTCSDiffusion {
 
 public:
-
-    static const int BC_NEUMANN;
-    static const int BC_DIRICHLET;
+  static const int BC_NEUMANN;
+  static const int BC_DIRICHLET;
 
   /*!
    * Create 1D-diffusion module.
@@ -56,13 +56,28 @@ public:
    * @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);
+  void simulate(std::vector<double> &c, std::vector<double> &alpha);
+
+  void setTimestep(double time_step);
+
+  void setBoundaryCondition(int index, double val, int type);
 
 private:
+  void simulate1D(std::vector<double> &c, double bc_left, double bc_right,
+                  std::vector<double> &alpha);
+  void simulate2D(std::vector<double> &c);
+  void simulate3D(std::vector<double> &c);
+
+  double getBCFromTuple(int index);
 
   boundary_condition bc;
 
+  Eigen::SparseMatrix<double> A_matrix;
+  Eigen::VectorXd b_vector;
+  Eigen::VectorXd x_vector;
+
+  double time_step;
+
   int grid_dim;
   int dim_x;
   int dim_y;
diff --git a/src/main.cpp b/src/main.cpp
index 2836c7c..9b376f5 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -18,13 +18,16 @@ int main(int argc, char *argv[]) {
 
   BTCSDiffusion diffu(x);
 
-  diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT);
+  diffu.setBoundaryCondition(0, 5. * std::pow(10, -6), BTCSDiffusion::BC_DIRICHLET);
+  diffu.setTimestep(1.);
+
+  // diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT);
   // 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.);
+    diffu.simulate(input, alpha);
   }
 
   return 0;

From 76640da6cb6b8765e45a769b6dac0280c7e8449d Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 6 Dec 2021 20:16:08 +0100
Subject: [PATCH 03/23] set boundary conditions hard

---
 src/BTCSDiffusion.cpp | 14 ++++++++++----
 src/BTCSDiffusion.hpp |  2 +-
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 802c9db..8023970 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -8,9 +8,11 @@
 #include <Eigen/src/OrderingMethods/Ordering.h>
 
 #include <algorithm>
+#include <cmath>
 #include <iomanip>
 #include <iostream>
 #include <tuple>
+#include <vector>
 
 const int BTCSDiffusion::BC_NEUMANN = 0;
 const int BTCSDiffusion::BC_DIRICHLET = 1;
@@ -125,6 +127,8 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
     // A_line++;
   }
 
+  std::cout << A_matrix << std::endl;
+
   // Eigen::SparseMatrix<double> A(size, size);
   // A.setFromTriplets(tripletList.begin(), tripletList.end());
 
@@ -152,15 +156,17 @@ void BTCSDiffusion::setTimestep(double time_step) {
 void BTCSDiffusion::simulate(std::vector<double> &c,
                              std::vector<double> &alpha) {
   if (this->grid_dim == 1) {
-    double bc_left = getBCFromTuple(0);
-    double bc_right = getBCFromTuple(1);
+    // double bc_left = getBCFromTuple(0);
+    // double bc_right = getBCFromTuple(1);
+    double bc_left = 5. * std::pow(10,-6);
+    double bc_right = c[this->dim_x -1];
 
     simulate1D(c, bc_left, bc_right, alpha);
   }
 }
 
-double BTCSDiffusion::getBCFromTuple(int index) {
-  double val = std::get<1>(bc[index]);
+double BTCSDiffusion::getBCFromTuple(int index, std::vector<double> &c) {
+  double val = std::get<0>(bc[index]);
 
   return val;
 }
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 27b2aed..afa724d 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -68,7 +68,7 @@ private:
   void simulate2D(std::vector<double> &c);
   void simulate3D(std::vector<double> &c);
 
-  double getBCFromTuple(int index);
+  double getBCFromTuple(int index, std::vector<double> &c);
 
   boundary_condition bc;
 

From 0a5c271a4b7f132ceb858b8879fe553c5a596549 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 6 Dec 2021 20:18:33 +0100
Subject: [PATCH 04/23] change values of bc @ A and b

---
 src/BTCSDiffusion.cpp | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 8023970..f6977f2 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -106,8 +106,11 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   A_matrix.resize(size, size);
   A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
 
-  A_matrix.insert(0, 0) = bc_left;
-  A_matrix.insert(size - 1, size - 1) = bc_right;
+  A_matrix.insert(0, 0) = 1;
+  A_matrix.insert(size - 1, size - 1) = 1;
+
+  b_vector[0] = bc_left;
+  b_vector[size-1] = bc_right;
 
   for (int i = 1; i < this->dim_x + 1; i++) {
     double sx = (alpha[i - 1] * time_step) / (dx * dx);

From 496e4ee3bcb2f17db6c544b1e380b425c54e5924 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 13:23:42 +0100
Subject: [PATCH 05/23] use internal boundary condition mechanism

---
 src/BTCSDiffusion.cpp | 36 ++++++++++++++++++++++++------------
 src/BTCSDiffusion.hpp |  7 ++++---
 2 files changed, 28 insertions(+), 15 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index f6977f2..a261ee3 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -97,7 +97,9 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   //   b[A_line] = this->bc[1];
 
   /*
-   * Begin to solve the equation system
+   * 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
@@ -110,7 +112,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   A_matrix.insert(size - 1, size - 1) = 1;
 
   b_vector[0] = bc_left;
-  b_vector[size-1] = bc_right;
+  b_vector[size - 1] = bc_right;
 
   for (int i = 1; i < this->dim_x + 1; i++) {
     double sx = (alpha[i - 1] * time_step) / (dx * dx);
@@ -130,7 +132,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
     // A_line++;
   }
 
-  std::cout << A_matrix << std::endl;
+  // std::cout << A_matrix << std::endl;
 
   // Eigen::SparseMatrix<double> A(size, size);
   // A.setFromTriplets(tripletList.begin(), tripletList.end());
@@ -159,22 +161,32 @@ void BTCSDiffusion::setTimestep(double time_step) {
 void BTCSDiffusion::simulate(std::vector<double> &c,
                              std::vector<double> &alpha) {
   if (this->grid_dim == 1) {
-    // double bc_left = getBCFromTuple(0);
-    // double bc_right = getBCFromTuple(1);
-    double bc_left = 5. * std::pow(10,-6);
-    double bc_right = c[this->dim_x -1];
+    double bc_left = getBCFromTuple(0, c[0]);
+    double bc_right = getBCFromTuple(1, c[c.size() - 1]);
+    // double bc_left = 5. * std::pow(10,-6);
+    // double bc_right = c[this->dim_x -1];
 
     simulate1D(c, bc_left, bc_right, alpha);
   }
 }
 
-double BTCSDiffusion::getBCFromTuple(int index, std::vector<double> &c) {
-  double val = std::get<0>(bc[index]);
+double BTCSDiffusion::getBCFromTuple(int index, double nearest_value) {
+  double val;
+  int type = std::get<0>(bc[index]);
+
+  if (type == BTCSDiffusion::BC_NEUMANN) {
+    // TODO implement gradient here
+    val = nearest_value;
+  } else if (type == BTCSDiffusion::BC_DIRICHLET) {
+    val = std::get<1>(bc[index]);
+  } else {
+    // some error handling here. Type was set to wrong value.
+  }
 
   return val;
 }
 
-void BTCSDiffusion::setBoundaryCondition(int index, double val, int type) {
-  std::get<0>(bc[index]) = val;
-  std::get<1>(bc[index]) = type;
+void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) {
+  std::get<0>(bc[index]) = type;
+  std::get<1>(bc[index]) = val;
 }
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index afa724d..11864ec 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -12,7 +12,8 @@
  */
 typedef Eigen::Triplet<double> T;
 
-typedef std::vector<std::tuple<int, double>> boundary_condition;
+typedef int bctype;
+typedef std::vector<std::tuple<bctype, double>> boundary_condition;
 
 /*!
  * Class implementing a solution for a 1/2/3D diffusion equation using backward
@@ -60,7 +61,7 @@ public:
 
   void setTimestep(double time_step);
 
-  void setBoundaryCondition(int index, double val, int type);
+  void setBoundaryCondition(int index, double val, bctype type);
 
 private:
   void simulate1D(std::vector<double> &c, double bc_left, double bc_right,
@@ -68,7 +69,7 @@ private:
   void simulate2D(std::vector<double> &c);
   void simulate3D(std::vector<double> &c);
 
-  double getBCFromTuple(int index, std::vector<double> &c);
+  double getBCFromTuple(int index, double nearest_value);
 
   boundary_condition bc;
 

From 74ab002d496d2d529e75b82fdc42a4cfc2791f8a Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 13:26:04 +0100
Subject: [PATCH 06/23] change datatype of tuple in constructor

---
 src/BTCSDiffusion.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index a261ee3..efe7e05 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -21,7 +21,7 @@ BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) {
   this->grid_dim = 1;
 
   // per default use Neumann condition with gradient of 0 at the end of the grid
-  this->bc.resize(2, std::tuple<int, double>(0, 0.));
+  this->bc.resize(2, std::tuple<bctype, double>(0, 0.));
 }
 BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) {
 

From e62545592364bc38227fe85125714ef60bf7a957 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 13:47:32 +0100
Subject: [PATCH 07/23] cleanup

---
 src/BTCSDiffusion.cpp | 58 -------------------------------------------
 1 file changed, 58 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index efe7e05..9f70440 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -50,52 +50,6 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   b_vector.resize(size);
   x_vector.resize(size);
 
-  // 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);
-
-  //   tripletList.push_back(T(A_line, i, (-1. - 2. * sx)));
-
-  //   tripletList.push_back(T(A_line, i - 1, sx));
-  //   tripletList.push_back(T(A_line, i + 1, sx));
-
-  //   b[A_line] = -c[i - 1];
-  //   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];
-
-  // A_line++;
-  // tripletList.push_back(T(A_line, size - 1, 1));
-  // // b[A_line] = bc[1];
-  // if (bc[1] == -1)
-  //   b[A_line] = c[c.size() - 1];
-  // else
-  //   b[A_line] = this->bc[1];
-
   /*
    * Begin to solve the equation system using LU solver of Eigen.
    *
@@ -123,20 +77,8 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
 
     b_vector[i] = -c[i - 1];
 
-    // tripletList.push_back(T(A_line, i, (-1. - 2. * sx)));
-
-    // tripletList.push_back(T(A_line, i - 1, sx));
-    // tripletList.push_back(T(A_line, i + 1, sx));
-
-    // b[A_line] = -c[i - 1];
-    // A_line++;
   }
 
-  // std::cout << A_matrix << std::endl;
-
-  // Eigen::SparseMatrix<double> A(size, size);
-  // A.setFromTriplets(tripletList.begin(), tripletList.end());
-
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       solver;
   solver.analyzePattern(A_matrix);

From 1bde33cd793d87ca3dd0dd445cd3dfe37cdb7e47 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 13:48:21 +0100
Subject: [PATCH 08/23] avoid return of unitialized variable

---
 src/BTCSDiffusion.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 9f70440..ab7d2ea 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -113,7 +113,7 @@ void BTCSDiffusion::simulate(std::vector<double> &c,
 }
 
 double BTCSDiffusion::getBCFromTuple(int index, double nearest_value) {
-  double val;
+  double val = -1;
   int type = std::get<0>(bc[index]);
 
   if (type == BTCSDiffusion::BC_NEUMANN) {
@@ -122,7 +122,7 @@ double BTCSDiffusion::getBCFromTuple(int index, double nearest_value) {
   } else if (type == BTCSDiffusion::BC_DIRICHLET) {
     val = std::get<1>(bc[index]);
   } else {
-    // some error handling here. Type was set to wrong value.
+    // TODO: implement error handling here. Type was set to wrong value.
   }
 
   return val;

From dcecc2dd723d656ca6d82b9a46d5ee117cd5af7f Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 13:59:08 +0100
Subject: [PATCH 09/23] added comments for public methods

---
 src/BTCSDiffusion.hpp | 35 ++++++++++++++++++++++++++++++++++-
 1 file changed, 34 insertions(+), 1 deletion(-)

diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 11864ec..c5ea7bb 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -12,7 +12,20 @@
  */
 typedef Eigen::Triplet<double> T;
 
+/*!
+ * Defines both types of boundary condition as a datatype.
+ */
 typedef int bctype;
+
+/*!
+ * 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;
 
 /*!
@@ -22,7 +35,13 @@ typedef std::vector<std::tuple<bctype, double>> boundary_condition;
 class BTCSDiffusion {
 
 public:
+  /*!
+   * Defines a Neumann boundary condition.
+   */
   static const int BC_NEUMANN;
+  /*!
+   * Defines a Dirichlet boundary condition.
+   */
   static const int BC_DIRICHLET;
 
   /*!
@@ -55,12 +74,26 @@ public:
    * @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);
 
+  /*!
+   * Set the timestep of the simulation
+   *
+   * @param time_step Time step (in seconds ???)
+   */
   void setTimestep(double time_step);
 
+  /*!
+   * Set the boundary condition of the given grid. This is done by defining an
+   * 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 val Value of the boundary condition (gradient for Neumann, exact
+   * value for Dirichlet).
+   * @param Type of the grid cell.
+   */
   void setBoundaryCondition(int index, double val, bctype type);
 
 private:

From 7cb18e62c1f60ba15fbdcb988c62281430ef9829 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 14:13:48 +0100
Subject: [PATCH 10/23] update header files

---
 src/BTCSDiffusion.cpp |  6 ------
 src/BTCSDiffusion.hpp |  1 -
 src/main.cpp          | 12 ++++++------
 3 files changed, 6 insertions(+), 13 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index ab7d2ea..369c525 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -1,14 +1,8 @@
 #include "BTCSDiffusion.hpp"
 
-#include <Eigen/SparseCholesky>
 #include <Eigen/SparseLU>
-#include <Eigen/SparseQR>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
-#include <Eigen/src/OrderingMethods/Ordering.h>
 
 #include <algorithm>
-#include <cmath>
 #include <iomanip>
 #include <iostream>
 #include <tuple>
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index c5ea7bb..c64e9db 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -2,7 +2,6 @@
 #define BTCSDIFFUSION_H_
 
 #include <Eigen/Sparse>
-#include <Eigen/src/SparseCore/SparseMatrixBase.h>
 #include <tuple>
 #include <vector>
 
diff --git a/src/main.cpp b/src/main.cpp
index 9b376f5..b4ef449 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,8 +1,7 @@
-#include "BTCSDiffusion.hpp"
-#include <cmath>
-#include <iostream>
-#include <vector>
-
+#include "BTCSDiffusion.hpp" // for BTCSDiffusion, BTCSDiffusion::BC_DIRICHLET
+#include <algorithm>         // for copy, max
+#include <iostream>          // for std
+#include <vector>            // for vector
 using namespace std;
 
 int main(int argc, char *argv[]) {
@@ -18,7 +17,8 @@ int main(int argc, char *argv[]) {
 
   BTCSDiffusion diffu(x);
 
-  diffu.setBoundaryCondition(0, 5. * std::pow(10, -6), BTCSDiffusion::BC_DIRICHLET);
+  diffu.setBoundaryCondition(0, 5. * std::pow(10, -6),
+                             BTCSDiffusion::BC_DIRICHLET);
   diffu.setTimestep(1.);
 
   // diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT);

From 30a9dbafb8410d666efe20527369c0fe5f3d8115 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 14:13:58 +0100
Subject: [PATCH 11/23] update gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index a072bc9..1e841cf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ build/
 compile_commands.json
 .cache/
 .ccls-cache/
+/iwyu/

From ca0fe9678b0fa37a0bfffe1e3538971bd212bdd8 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 19:36:06 +0100
Subject: [PATCH 12/23] implement changes as discussed in meeting on 12-13-2021

---
 src/BTCSDiffusion.cpp | 41 ++++++++++++++++++++---------------------
 src/BTCSDiffusion.hpp | 21 ++++++++++++---------
 2 files changed, 32 insertions(+), 30 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 369c525..d5d7ac1 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -11,13 +11,14 @@
 const int BTCSDiffusion::BC_NEUMANN = 0;
 const int BTCSDiffusion::BC_DIRICHLET = 1;
 
-BTCSDiffusion::BTCSDiffusion(int x) : dim_x(x) {
+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>(0, 0.));
+  this->bc.resize(2, std::tuple<bctype, double>(BTCSDiffusion::BC_NEUMANN, 0.));
 }
-BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) {
+BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) {
 
   // this->grid_dim = 2;
 
@@ -25,20 +26,19 @@ BTCSDiffusion::BTCSDiffusion(int x, int y) : dim_x(x), dim_y(y) {
   // // 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) {
+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, double bc_left,
-                               double bc_right, std::vector<double> &alpha) {
-  // calculate dx
-  double dx = 1. / (this->dim_x - 1);
+                               double bc_right,
+                               const std::vector<double> &alpha, double dx,
+                               int size) {
 
-  // calculate size needed for A matrix and b,x vectors
-  int size = this->dim_x + 2;
+  // we need 2 more grid cells for ghost cells
+  size = size + 2;
 
   // set sizes of private and yet allocated vectors
   b_vector.resize(size);
@@ -62,7 +62,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   b_vector[0] = bc_left;
   b_vector[size - 1] = bc_right;
 
-  for (int i = 1; i < this->dim_x + 1; i++) {
+  for (int i = 1; i < this->n_x + 1; i++) {
     double sx = (alpha[i - 1] * time_step) / (dx * dx);
 
     A_matrix.insert(i, i) = -1. - 2. * sx;
@@ -70,7 +70,6 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
     A_matrix.insert(i, i + 1) = sx;
 
     b_vector[i] = -c[i - 1];
-
   }
 
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
@@ -95,24 +94,24 @@ void BTCSDiffusion::setTimestep(double time_step) {
 }
 
 void BTCSDiffusion::simulate(std::vector<double> &c,
-                             std::vector<double> &alpha) {
+                             const std::vector<double> &alpha) {
   if (this->grid_dim == 1) {
-    double bc_left = getBCFromTuple(0, c[0]);
-    double bc_right = getBCFromTuple(1, c[c.size() - 1]);
-    // double bc_left = 5. * std::pow(10,-6);
-    // double bc_right = c[this->dim_x -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_left, bc_right, alpha);
+    simulate1D(c, bc_left, bc_right, alpha, this->dx, this->n_x);
   }
 }
 
-double BTCSDiffusion::getBCFromTuple(int index, double nearest_value) {
+double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
+                                     double neighbor_alpha) {
   double val = -1;
   int type = std::get<0>(bc[index]);
 
   if (type == BTCSDiffusion::BC_NEUMANN) {
-    // TODO implement gradient here
-    val = nearest_value;
+    val = neighbor_c + (this->time_step / (dx * dx)) * neighbor_alpha *
+                           std::get<1>(bc[index]);
   } else if (type == BTCSDiffusion::BC_DIRICHLET) {
     val = std::get<1>(bc[index]);
   } else {
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index c64e9db..355c646 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -48,7 +48,7 @@ public:
    *
    * @param x Count of cells in x direction.
    */
-  BTCSDiffusion(int x);
+  explicit BTCSDiffusion(int x);
 
   /*!
    * Currently not implemented: Create 2D-diffusion module.
@@ -56,7 +56,7 @@ public:
    * @param x Count of cells in x direction.
    * @param y Count of cells in y direction.
    */
-  BTCSDiffusion(int x, int y);
+  explicit BTCSDiffusion(int x, int y);
 
   /*!
    * Currently not implemented: Create 3D-diffusion module.
@@ -65,7 +65,7 @@ public:
    * @param y Count of cells in y direction.
    * @param z Count of cells in z direction.
    */
-  BTCSDiffusion(int x, int y, int z);
+  explicit BTCSDiffusion(int x, int y, int z);
 
   /*!
    * With given ghost zones simulate diffusion. Only 1D allowed at this moment.
@@ -74,7 +74,7 @@ public:
    * continious memory (Row-wise).
    * @param alpha Vector of diffusioncoefficients for each grid element.
    */
-  void simulate(std::vector<double> &c, std::vector<double> &alpha);
+  void simulate(std::vector<double> &c, const std::vector<double> &alpha);
 
   /*!
    * Set the timestep of the simulation
@@ -97,11 +97,11 @@ public:
 
 private:
   void simulate1D(std::vector<double> &c, double bc_left, double bc_right,
-                  std::vector<double> &alpha);
+                  const std::vector<double> &alpha, double dx, int size);
   void simulate2D(std::vector<double> &c);
   void simulate3D(std::vector<double> &c);
 
-  double getBCFromTuple(int index, double nearest_value);
+  double getBCFromTuple(int index, double nearest_value, double neighbor_alpha);
 
   boundary_condition bc;
 
@@ -112,9 +112,12 @@ private:
   double time_step;
 
   int grid_dim;
-  int dim_x;
-  int dim_y;
-  int dim_z;
+  int n_x;
+  double dx;
+  int n_y;
+  double dy;
+  int n_z;
+  double dz;
 };
 
 #endif // BTCSDIFFUSION_H_

From a1bdcf84a78b8529c41f7da595244da3c0962eef Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 13 Dec 2021 19:38:38 +0100
Subject: [PATCH 13/23] cleanup of main file

---
 src/main.cpp | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index b4ef449..9735f11 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,26 +6,25 @@ using namespace std;
 
 int main(int argc, char *argv[]) {
 
+  // count of grid cells
   int x = 20;
 
+  // create input + diffusion coefficients for each grid cell
   std::vector<double> alpha(x, 1 * pow(10, -1));
   std::vector<double> input(x, 1 * std::pow(10, -6));
-  std::vector<double> bc_left, bc_right;
-
-  bc_left.push_back(5. * std::pow(10, -6));
-  bc_right.push_back(-1);
 
+  // create instance of diffusion module
   BTCSDiffusion diffu(x);
 
+  // set the boundary condition for the left ghost cell to dirichlet
   diffu.setBoundaryCondition(0, 5. * std::pow(10, -6),
                              BTCSDiffusion::BC_DIRICHLET);
+
+  // set timestep for simulation to 1 second
   diffu.setTimestep(1.);
 
-  // diffu.setBoundaryCondition(bc_left, BTCSDiffusion::LEFT);
-  // we don't need this since Neumann condition with gradient of 0 is set per
-  // default
-  // diffu.setBoundaryCondition(bc_right, BTCSDiffusion::RIGHT);
-
+  // loop 100 times
+  // output is currently generated by the method itself
   for (int i = 0; i < 100; i++) {
     diffu.simulate(input, alpha);
   }

From e8dae917d5236a6c80e888eea805e92e5eb82e7f Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Fri, 14 Jan 2022 12:36:39 +0100
Subject: [PATCH 14/23] replace names of boundary condition variables

---
 src/BTCSDiffusion.cpp | 11 ++++++-----
 src/BTCSDiffusion.hpp | 14 ++++++++++----
 2 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index d5d7ac1..1db5f60 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -8,15 +8,16 @@
 #include <tuple>
 #include <vector>
 
-const int BTCSDiffusion::BC_NEUMANN = 0;
-const int BTCSDiffusion::BC_DIRICHLET = 1;
+const int BTCSDiffusion::BC_CONSTANT = 0;
+const int BTCSDiffusion::BC_CLOSED = 1;
+const int BTCSDiffusion::BC_FLUX = 2;
 
 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_NEUMANN, 0.));
+  this->bc.resize(2, std::tuple<bctype, double>(BTCSDiffusion::BC_CONSTANT, 0.));
 }
 BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) {
 
@@ -109,10 +110,10 @@ double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
   double val = -1;
   int type = std::get<0>(bc[index]);
 
-  if (type == BTCSDiffusion::BC_NEUMANN) {
+  if (type == BTCSDiffusion::BC_CONSTANT) {
     val = neighbor_c + (this->time_step / (dx * dx)) * neighbor_alpha *
                            std::get<1>(bc[index]);
-  } else if (type == BTCSDiffusion::BC_DIRICHLET) {
+  } else if (type == BTCSDiffusion::BC_CLOSED) {
     val = std::get<1>(bc[index]);
   } else {
     // TODO: implement error handling here. Type was set to wrong value.
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 355c646..9720ec9 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -35,13 +35,19 @@ class BTCSDiffusion {
 
 public:
   /*!
-   * Defines a Neumann boundary condition.
+   * Defines a constant/Dirichlet boundary condition.
    */
-  static const int BC_NEUMANN;
+  static const int BC_CONSTANT;
+
   /*!
-   * Defines a Dirichlet boundary condition.
+   * Defines a closed/Neumann boundary condition.
    */
-  static const int BC_DIRICHLET;
+    static const int BC_CLOSED;
+
+    /*!
+     * Defines a flux/Cauchy boundary condition.
+     */
+    static const int BC_FLUX;
 
   /*!
    * Create 1D-diffusion module.

From 29fc70ce1ab18ae46cc8d0fb841e2af190f755c5 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 17 Jan 2022 14:20:53 +0100
Subject: [PATCH 15/23] rewrite initialization of module

ATTENTION: code will throw errors on compilation!
---
 src/BTCSDiffusion.cpp | 63 ++++++++++++++++++++++++++++++++-----------
 src/BTCSDiffusion.hpp | 58 ++++++++++++++++++++-------------------
 2 files changed, 79 insertions(+), 42 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 1db5f60..1b45749 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -3,6 +3,7 @@
 #include <Eigen/SparseLU>
 
 #include <algorithm>
+#include <cassert>
 #include <iomanip>
 #include <iostream>
 #include <tuple>
@@ -12,26 +13,58 @@ const int BTCSDiffusion::BC_CONSTANT = 0;
 const int BTCSDiffusion::BC_CLOSED = 1;
 const int BTCSDiffusion::BC_FLUX = 2;
 
-BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) {
-  this->grid_dim = 1;
-  this->dx = 1. / (x - 1);
+BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
+  assert(dim <= 3);
 
-  // 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.));
+  grid_cells.resize(dim, 1);
+  spatial_discretization.resize(dim, 1);
+  deltas.resize(dim, 1);
 }
-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);
+std::vector<int> BTCSDiffusion::getNumberOfGridCells() {
+  return this->grid_cells;
 }
-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
+std::vector<int> BTCSDiffusion::getSpatialDiscretization() {
+  return this->spatial_discretization;
 }
+void BTCSDiffusion::setNumberOfGridCells(std::vector<int> &n_grid) {
+  grid_cells = n_grid;
+  assert(grid_cells.size() == grid_dim);
+  updateDeltas();
+}
+void BTCSDiffusion::setSpatialDiscretization(std::vector<int> &s_grid) {
+  spatial_discretization = s_grid;
+  assert(spatial_discretization.size() == grid_dim);
+  updateDeltas();
+}
+
+void BTCSDiffusion::updateDeltas() {
+  for (int i = 0; i < grid_dim; i++) {
+    deltas[i] = (double)spatial_discretization[i] / grid_cells[i];
+  }
+}
+// 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, double bc_left,
                                double bc_right,
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 9720ec9..e849224 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -42,36 +42,41 @@ public:
   /*!
    * Defines a closed/Neumann boundary condition.
    */
-    static const int BC_CLOSED;
+  static const int BC_CLOSED;
 
-    /*!
-     * Defines a flux/Cauchy boundary condition.
-     */
-    static const int BC_FLUX;
+  /*!
+   * Defines a flux/Cauchy boundary condition.
+   */
+  static const int BC_FLUX;
 
   /*!
    * Create 1D-diffusion module.
    *
    * @param x Count of cells in x direction.
    */
-  explicit BTCSDiffusion(int x);
+  BTCSDiffusion(unsigned int dim);
 
-  /*!
-   * 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);
+  std::vector<int> getNumberOfGridCells();
+  std::vector<int> getSpatialDiscretization();
+  void setNumberOfGridCells(std::vector<int> &n_grid);
+  void setSpatialDiscretization(std::vector<int> &s_grid);
 
-  /*!
-   * 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.
-   */
-  explicit BTCSDiffusion(int x, int y, int z);
+  // /*!
+  //  * 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.
+  //  *
+  //  * @param x Count of cells in x direction.
+  //  * @param y Count of cells in y direction.
+  //  * @param z Count of cells in z direction.
+  //  */
+  // explicit BTCSDiffusion(int x, int y, int z);
 
   /*!
    * With given ghost zones simulate diffusion. Only 1D allowed at this moment.
@@ -109,6 +114,8 @@ private:
 
   double getBCFromTuple(int index, double nearest_value, double neighbor_alpha);
 
+  void updateDeltas();
+
   boundary_condition bc;
 
   Eigen::SparseMatrix<double> A_matrix;
@@ -118,12 +125,9 @@ private:
   double time_step;
 
   int grid_dim;
-  int n_x;
-  double dx;
-  int n_y;
-  double dy;
-  int n_z;
-  double dz;
+  std::vector<int> grid_cells;
+  std::vector<int> spatial_discretization;
+  std::vector<double> deltas;
 };
 
 #endif // BTCSDIFFUSION_H_

From 3fa39fdc367c974eee5aeecc3af9525d02f643df Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 17 Jan 2022 14:33:43 +0100
Subject: [PATCH 16/23] update main + library variables

STILL NO RUNNABLE CODE!
---
 src/BTCSDiffusion.cpp | 27 ++++++++++++++-------------
 src/main.cpp          | 21 +++++++++++++++------
 2 files changed, 29 insertions(+), 19 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 1b45749..372caef 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -75,8 +75,8 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   size = size + 2;
 
   // set sizes of private and yet allocated vectors
-  b_vector.resize(size);
-  x_vector.resize(size);
+  b_vector.resize(size + 2);
+  x_vector.resize(size + 2);
 
   /*
    * Begin to solve the equation system using LU solver of Eigen.
@@ -87,23 +87,23 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
    * TODO: remove output
    */
 
-  A_matrix.resize(size, size);
-  A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
+  A_matrix.resize(size + 2, size + 2);
+  A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, 3));
 
   A_matrix.insert(0, 0) = 1;
-  A_matrix.insert(size - 1, size - 1) = 1;
+  A_matrix.insert(1, 1) = 1;
 
   b_vector[0] = bc_left;
-  b_vector[size - 1] = bc_right;
+  b_vector[1] = bc_right;
 
-  for (int i = 1; i < this->n_x + 1; i++) {
-    double sx = (alpha[i - 1] * time_step) / (dx * dx);
+  for (int i = 2; i < size + 2; i++) {
+    double sx = (alpha[(i - 2) - 1] * time_step) / (dx * dx);
 
     A_matrix.insert(i, i) = -1. - 2. * sx;
     A_matrix.insert(i, i - 1) = sx;
     A_matrix.insert(i, i + 1) = sx;
 
-    b_vector[i] = -c[i - 1];
+    b_vector[i] = -c[(i - 2)];
   }
 
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
@@ -119,7 +119,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
 
   for (int i = 0; i < c.size(); i++) {
-    c[i] = x_vector[i + 1];
+    c[i] = x_vector[i + 2];
   }
 }
 
@@ -134,7 +134,8 @@ void BTCSDiffusion::simulate(std::vector<double> &c,
     double bc_right =
         getBCFromTuple(1, c[c.size() - 1], alpha[alpha.size() - 1]);
 
-    simulate1D(c, bc_left, bc_right, alpha, this->dx, this->n_x);
+    simulate1D(c, bc_left, bc_right, alpha, this->deltas[0],
+               this->grid_cells[0]);
   }
 }
 
@@ -144,8 +145,8 @@ double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
   int type = std::get<0>(bc[index]);
 
   if (type == BTCSDiffusion::BC_CONSTANT) {
-    val = neighbor_c + (this->time_step / (dx * dx)) * neighbor_alpha *
-                           std::get<1>(bc[index]);
+    val = neighbor_c + (this->time_step / (this->deltas[0] * this->deltas[0])) *
+                           neighbor_alpha * std::get<1>(bc[index]);
   } else if (type == BTCSDiffusion::BC_CLOSED) {
     val = std::get<1>(bc[index]);
   } else {
diff --git a/src/main.cpp b/src/main.cpp
index 9735f11..d556c64 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,19 +6,28 @@ using namespace std;
 
 int main(int argc, char *argv[]) {
 
-  // count of grid cells
-  int x = 20;
+  // dimension of grid
+  int dim = 1;
+
+  int n = 20;
 
   // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(x, 1 * pow(10, -1));
-  std::vector<double> input(x, 1 * std::pow(10, -6));
+  std::vector<double> alpha(n, 1 * pow(10, -1));
+  std::vector<double> input(n, 1 * std::pow(10, -6));
 
   // create instance of diffusion module
-  BTCSDiffusion diffu(x);
+  BTCSDiffusion diffu(dim);
+
+  std::vector<int> vec_n = diffu.getNumberOfGridCells();
+
+
+  vec_n[0] = n;
+
+  diffu.setNumberOfGridCells(vec_n);
 
   // set the boundary condition for the left ghost cell to dirichlet
   diffu.setBoundaryCondition(0, 5. * std::pow(10, -6),
-                             BTCSDiffusion::BC_DIRICHLET);
+                             BTCSDiffusion::BC_CONSTANT);
 
   // set timestep for simulation to 1 second
   diffu.setTimestep(1.);

From 52c1f472f63a7298920e55f34889fc292b852860 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Wed, 19 Jan 2022 18:07:25 +0100
Subject: [PATCH 17/23] rewrite to runnable code

---
 src/BTCSDiffusion.cpp | 41 +++++++++++++++++++++++++++--------------
 src/BTCSDiffusion.hpp |  2 +-
 src/main.cpp          |  1 -
 3 files changed, 28 insertions(+), 16 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 372caef..2fe91ce 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -30,18 +30,30 @@ std::vector<int> BTCSDiffusion::getSpatialDiscretization() {
 void BTCSDiffusion::setNumberOfGridCells(std::vector<int> &n_grid) {
   grid_cells = n_grid;
   assert(grid_cells.size() == grid_dim);
-  updateDeltas();
+  updateInternals();
 }
 void BTCSDiffusion::setSpatialDiscretization(std::vector<int> &s_grid) {
   spatial_discretization = s_grid;
   assert(spatial_discretization.size() == grid_dim);
-  updateDeltas();
+  updateInternals();
 }
 
-void BTCSDiffusion::updateDeltas() {
+void BTCSDiffusion::updateInternals() {
   for (int i = 0; i < grid_dim; i++) {
     deltas[i] = (double)spatial_discretization[i] / grid_cells[i];
   }
+
+  switch (grid_dim) {
+  case 1:
+    bc.resize(2, std::tuple<bctype, double>(BTCSDiffusion::BC_CLOSED, 0.));
+    break;
+  case 2:
+    bc.resize(2 * grid_cells[0] + 2 * grid_cells[1], std::tuple<bctype, double>(BTCSDiffusion::BC_CLOSED, 0.));
+    break;
+  case 3:
+    // TODO
+    break;
+  }
 }
 // BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) {
 //   this->grid_dim = 1;
@@ -72,7 +84,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
                                int size) {
 
   // we need 2 more grid cells for ghost cells
-  size = size + 2;
+  // size = size + 2;
 
   // set sizes of private and yet allocated vectors
   b_vector.resize(size + 2);
@@ -91,19 +103,19 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, 3));
 
   A_matrix.insert(0, 0) = 1;
-  A_matrix.insert(1, 1) = 1;
+  A_matrix.insert(size + 1, size + 1) = 1;
 
   b_vector[0] = bc_left;
-  b_vector[1] = bc_right;
+  b_vector[size + 1] = bc_right;
 
-  for (int i = 2; i < size + 2; i++) {
-    double sx = (alpha[(i - 2) - 1] * time_step) / (dx * dx);
+  for (int i = 1; i < size + 1; i++) {
+    double sx = (alpha[i - 1] * time_step) / (dx * dx);
 
     A_matrix.insert(i, i) = -1. - 2. * sx;
     A_matrix.insert(i, i - 1) = sx;
     A_matrix.insert(i, i + 1) = sx;
 
-    b_vector[i] = -c[(i - 2)];
+    b_vector[i] = -c[i - 1];
   }
 
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
@@ -119,7 +131,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
 
   for (int i = 0; i < c.size(); i++) {
-    c[i] = x_vector[i + 2];
+    c[i] = x_vector[i + 1];
   }
 }
 
@@ -144,10 +156,11 @@ double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
   double val = -1;
   int type = std::get<0>(bc[index]);
 
-  if (type == BTCSDiffusion::BC_CONSTANT) {
-    val = neighbor_c + (this->time_step / (this->deltas[0] * this->deltas[0])) *
-                           neighbor_alpha * std::get<1>(bc[index]);
-  } else if (type == BTCSDiffusion::BC_CLOSED) {
+  if (type == BTCSDiffusion::BC_CLOSED) {
+    val = neighbor_c;
+    // val = neighbor_c + (this->time_step / (this->deltas[0] * this->deltas[0])) *
+    //                        neighbor_alpha * std::get<1>(bc[index]);
+  } else if (type == BTCSDiffusion::BC_CONSTANT){
     val = std::get<1>(bc[index]);
   } else {
     // TODO: implement error handling here. Type was set to wrong value.
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index e849224..7d0b96f 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -114,7 +114,7 @@ private:
 
   double getBCFromTuple(int index, double nearest_value, double neighbor_alpha);
 
-  void updateDeltas();
+  void updateInternals();
 
   boundary_condition bc;
 
diff --git a/src/main.cpp b/src/main.cpp
index d556c64..356f89d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -20,7 +20,6 @@ int main(int argc, char *argv[]) {
 
   std::vector<int> vec_n = diffu.getNumberOfGridCells();
 
-
   vec_n[0] = n;
 
   diffu.setNumberOfGridCells(vec_n);

From c3d82afed4b27d0ab2a9114ecc6623902b98b362 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Thu, 20 Jan 2022 09:35:51 +0100
Subject: [PATCH 18/23] swap typedefs into class definition

---
 src/BTCSDiffusion.hpp | 33 +++++++++++++++++----------------
 1 file changed, 17 insertions(+), 16 deletions(-)

diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 7d0b96f..2df15c4 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -5,28 +5,12 @@
 #include <tuple>
 #include <vector>
 
-/*!
- * Datatype to fill the sparse matrix which is used to solve the equation
- * system.
- */
-typedef Eigen::Triplet<double> T;
 
 /*!
  * Defines both types of boundary condition as a datatype.
  */
 typedef int bctype;
 
-/*!
- * 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;
-
 /*!
  * Class implementing a solution for a 1/2/3D diffusion equation using backward
  * euler.
@@ -49,6 +33,23 @@ public:
    */
   static const int BC_FLUX;
 
+  /*!
+   * 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.
    *

From e675381683950f4fbd7d5729b9f8cd70a2b38261 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Thu, 20 Jan 2022 09:41:34 +0100
Subject: [PATCH 19/23] Change boundary_condition to struct instead of tuple

---
 src/BTCSDiffusion.cpp | 22 ++++++++++++++--------
 src/BTCSDiffusion.hpp | 18 ++++++++++++++++--
 2 files changed, 30 insertions(+), 10 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 2fe91ce..c9b0837 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -45,10 +45,11 @@ void BTCSDiffusion::updateInternals() {
 
   switch (grid_dim) {
   case 1:
-    bc.resize(2, std::tuple<bctype, double>(BTCSDiffusion::BC_CLOSED, 0.));
+    bc.resize(2, {BTCSDiffusion::BC_CLOSED, 0});
     break;
   case 2:
-    bc.resize(2 * grid_cells[0] + 2 * grid_cells[1], std::tuple<bctype, double>(BTCSDiffusion::BC_CLOSED, 0.));
+    bc.resize(2 * grid_cells[0] + 2 * grid_cells[1],
+              {BTCSDiffusion::BC_CLOSED, 0});
     break;
   case 3:
     // TODO
@@ -154,14 +155,15 @@ void BTCSDiffusion::simulate(std::vector<double> &c,
 double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
                                      double neighbor_alpha) {
   double val = -1;
-  int type = std::get<0>(bc[index]);
+  int type = bc[index].type;
 
   if (type == BTCSDiffusion::BC_CLOSED) {
     val = neighbor_c;
-    // val = neighbor_c + (this->time_step / (this->deltas[0] * this->deltas[0])) *
+    // val = neighbor_c + (this->time_step / (this->deltas[0] *
+    // this->deltas[0])) *
     //                        neighbor_alpha * std::get<1>(bc[index]);
-  } else if (type == BTCSDiffusion::BC_CONSTANT){
-    val = std::get<1>(bc[index]);
+  } else if (type == BTCSDiffusion::BC_CONSTANT) {
+    val = bc[index].value;
   } else {
     // TODO: implement error handling here. Type was set to wrong value.
   }
@@ -170,6 +172,10 @@ double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
 }
 
 void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) {
-  std::get<0>(bc[index]) = type;
-  std::get<1>(bc[index]) = val;
+
+  bc[index].type = type;
+  bc[index].value = val;
+
+  // std::get<0>(bc[index]) = type;
+  // std::get<1>(bc[index]) = val;
 }
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 2df15c4..21a4b37 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -42,7 +42,21 @@ public:
    * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
    * the boundary condition
    */
-  typedef std::vector<std::tuple<bctype, double>> 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
@@ -117,7 +131,7 @@ private:
 
   void updateInternals();
 
-  boundary_condition bc;
+  std::vector<boundary_condition> bc;
 
   Eigen::SparseMatrix<double> A_matrix;
   Eigen::VectorXd b_vector;

From 5606b559c7b2643e53a62174a73122ccc1ce9671 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Thu, 20 Jan 2022 11:01:20 +0100
Subject: [PATCH 20/23] Refactor variable names

---
 src/BTCSDiffusion.cpp | 10 +++++-----
 src/BTCSDiffusion.hpp |  2 +-
 src/main.cpp          |  4 ++--
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index c9b0837..72f7242 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -17,7 +17,7 @@ BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
   assert(dim <= 3);
 
   grid_cells.resize(dim, 1);
-  spatial_discretization.resize(dim, 1);
+  domain_size.resize(dim, 1);
   deltas.resize(dim, 1);
 }
 
@@ -25,7 +25,7 @@ std::vector<int> BTCSDiffusion::getNumberOfGridCells() {
   return this->grid_cells;
 }
 std::vector<int> BTCSDiffusion::getSpatialDiscretization() {
-  return this->spatial_discretization;
+  return this->domain_size;
 }
 void BTCSDiffusion::setNumberOfGridCells(std::vector<int> &n_grid) {
   grid_cells = n_grid;
@@ -33,14 +33,14 @@ void BTCSDiffusion::setNumberOfGridCells(std::vector<int> &n_grid) {
   updateInternals();
 }
 void BTCSDiffusion::setSpatialDiscretization(std::vector<int> &s_grid) {
-  spatial_discretization = s_grid;
-  assert(spatial_discretization.size() == grid_dim);
+  domain_size = s_grid;
+  assert(domain_size.size() == grid_dim);
   updateInternals();
 }
 
 void BTCSDiffusion::updateInternals() {
   for (int i = 0; i < grid_dim; i++) {
-    deltas[i] = (double)spatial_discretization[i] / grid_cells[i];
+    deltas[i] = (double)domain_size[i] / grid_cells[i];
   }
 
   switch (grid_dim) {
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 21a4b37..6a9ebda 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -141,7 +141,7 @@ private:
 
   int grid_dim;
   std::vector<int> grid_cells;
-  std::vector<int> spatial_discretization;
+  std::vector<int> domain_size;
   std::vector<double> deltas;
 };
 
diff --git a/src/main.cpp b/src/main.cpp
index 356f89d..c4549ee 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -13,7 +13,7 @@ int main(int argc, char *argv[]) {
 
   // create input + diffusion coefficients for each grid cell
   std::vector<double> alpha(n, 1 * pow(10, -1));
-  std::vector<double> input(n, 1 * std::pow(10, -6));
+  std::vector<double> field(n, 1 * std::pow(10, -6));
 
   // create instance of diffusion module
   BTCSDiffusion diffu(dim);
@@ -34,7 +34,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(input, alpha);
+    diffu.simulate(field, alpha);
   }
 
   return 0;

From 2cd8a1e5d8a2a625e41304459c642e396204bc36 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Fri, 21 Jan 2022 14:09:30 +0100
Subject: [PATCH 21/23] Update simulation loop to new indexing

---
 src/BTCSDiffusion.cpp | 64 ++++++++++++++++++++++++-------------------
 src/BTCSDiffusion.hpp | 26 +++++++++---------
 2 files changed, 49 insertions(+), 41 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 72f7242..3251a04 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -79,17 +79,24 @@ void BTCSDiffusion::updateInternals() {
 //   // TODO: reserve memory for boundary conditions
 // }
 
-void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
-                               double bc_right,
+void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
+                               boundary_condition right,
                                const std::vector<double> &alpha, double dx,
                                int size) {
 
+  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;
 
+  int bc_offset = !left_is_constant + !right_is_constant;
+  ;
+
   // set sizes of private and yet allocated vectors
-  b_vector.resize(size + 2);
-  x_vector.resize(size + 2);
+  b_vector.resize(size + bc_offset);
+  x_vector.resize(size + bc_offset);
 
   /*
    * Begin to solve the equation system using LU solver of Eigen.
@@ -100,23 +107,27 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
    * TODO: remove output
    */
 
-  A_matrix.resize(size + 2, size + 2);
-  A_matrix.reserve(Eigen::VectorXi::Constant(size + 2, 3));
+  A_matrix.resize(size + bc_offset, size + bc_offset);
+  A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3));
 
   A_matrix.insert(0, 0) = 1;
-  A_matrix.insert(size + 1, size + 1) = 1;
+  b_vector[0] = (left_is_constant ? c[0] : getBCFromFlux(left, c[0], alpha[0]));
 
-  b_vector[0] = bc_left;
-  b_vector[size + 1] = bc_right;
+  A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1;
+  b_vector[size + bc_offset - 1] =
+    (right_is_constant ? c[size - 1] : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
 
-  for (int i = 1; i < size + 1; i++) {
-    double sx = (alpha[i - 1] * time_step) / (dx * dx);
+  // A_matrix.insert(0, 0) = 1;
+  // A_matrix.insert(size + 1, size + 1) = 1;
+
+  for (int i = 1; i < size - right_is_constant; i++) {
+    double sx = (alpha[i - !(left_is_constant)] * time_step) / (dx * dx);
 
     A_matrix.insert(i, i) = -1. - 2. * sx;
     A_matrix.insert(i, i - 1) = sx;
     A_matrix.insert(i, i + 1) = sx;
 
-    b_vector[i] = -c[i - 1];
+    b_vector[i] = -c[i - !(left_is_constant)];
   }
 
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
@@ -132,7 +143,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
 
   for (int i = 0; i < c.size(); i++) {
-    c[i] = x_vector[i + 1];
+    c[i] = x_vector[i + !left_is_constant];
   }
 }
 
@@ -143,27 +154,24 @@ 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]);
+    // 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_left, bc_right, alpha, this->deltas[0],
-               this->grid_cells[0]);
+    simulate1D(c, bc[0], bc[1], alpha, this->deltas[0], this->grid_cells[0]);
   }
 }
 
-double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
-                                     double neighbor_alpha) {
-  double val = -1;
-  int type = bc[index].type;
+inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, double neighbor_c,
+                                    double neighbor_alpha) {
 
-  if (type == BTCSDiffusion::BC_CLOSED) {
+  double val;
+
+  if (bc.type == BTCSDiffusion::BC_CLOSED) {
     val = neighbor_c;
-    // val = neighbor_c + (this->time_step / (this->deltas[0] *
-    // this->deltas[0])) *
-    //                        neighbor_alpha * std::get<1>(bc[index]);
-  } else if (type == BTCSDiffusion::BC_CONSTANT) {
-    val = bc[index].value;
+  } else if (bc.type == BTCSDiffusion::BC_FLUX) {
+    //TODO
+    // val = bc[index].value;
   } else {
     // TODO: implement error handling here. Type was set to wrong value.
   }
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index 6a9ebda..af08e0a 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -5,7 +5,6 @@
 #include <tuple>
 #include <vector>
 
-
 /*!
  * Defines both types of boundary condition as a datatype.
  */
@@ -42,10 +41,10 @@ public:
    * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
    * the boundary condition
    */
-    typedef struct boundary_condition {
-      bctype type;
-      double value;
-    } boundary_condition;
+  typedef struct boundary_condition {
+    bctype type;
+    double value;
+  } boundary_condition;
 
   /*!
    * A boundary condition consists of two features. A type and the according
@@ -58,11 +57,11 @@ public:
    */
   // 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;
+  /*!
+   * Datatype to fill the sparse matrix which is used to solve the equation
+   * system.
+   */
+  typedef Eigen::Triplet<double> T;
 
   /*!
    * Create 1D-diffusion module.
@@ -122,12 +121,13 @@ public:
   void setBoundaryCondition(int index, double val, bctype type);
 
 private:
-  void simulate1D(std::vector<double> &c, double bc_left, double bc_right,
-                  const std::vector<double> &alpha, double dx, int size);
+  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);
 
-  double getBCFromTuple(int index, double nearest_value, double neighbor_alpha);
+  inline double getBCFromFlux(boundary_condition bc, double nearest_value, double neighbor_alpha);
 
   void updateInternals();
 

From f6dbc3fb164ae0fba91cee6e46af7f41728ffec0 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Mon, 24 Jan 2022 12:28:29 +0100
Subject: [PATCH 22/23] Fix loop indexiation

---
 src/BTCSDiffusion.cpp | 13 ++++++++-----
 src/main.cpp          |  2 +-
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 3251a04..6c9d896 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -6,6 +6,7 @@
 #include <cassert>
 #include <iomanip>
 #include <iostream>
+#include <ostream>
 #include <tuple>
 #include <vector>
 
@@ -111,25 +112,27 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
   A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3));
 
   A_matrix.insert(0, 0) = 1;
-  b_vector[0] = (left_is_constant ? c[0] : getBCFromFlux(left, c[0], alpha[0]));
+  b_vector[0] = (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
 
   A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1;
   b_vector[size + bc_offset - 1] =
-    (right_is_constant ? c[size - 1] : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
+    (right_is_constant ? right.value : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
 
   // A_matrix.insert(0, 0) = 1;
   // A_matrix.insert(size + 1, size + 1) = 1;
 
   for (int i = 1; i < size - right_is_constant; i++) {
-    double sx = (alpha[i - !(left_is_constant)] * time_step) / (dx * dx);
+    double sx = (alpha[i + !(left_is_constant)] * time_step) / (dx * dx);
 
     A_matrix.insert(i, i) = -1. - 2. * sx;
     A_matrix.insert(i, i - 1) = sx;
     A_matrix.insert(i, i + 1) = sx;
 
-    b_vector[i] = -c[i - !(left_is_constant)];
+    b_vector[i] = -c[i + !(left_is_constant)];
   }
 
+std::cout << b_vector << "\n" << A_matrix << std::endl;
+
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       solver;
   solver.analyzePattern(A_matrix);
@@ -138,7 +141,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
 
   std::cout << solver.lastErrorMessage() << std::endl;
 
-  x_vector = solver.solve(b_vector);
+  // x_vector = solver.solve(b_vector);
 
   std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
 
diff --git a/src/main.cpp b/src/main.cpp
index c4549ee..31445a2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -33,7 +33,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++) {
+  for (int i = 0; i < 1; i++) {
     diffu.simulate(field, alpha);
   }
 

From 6308cd52af05c9b755cae44bfa166630a3230ae1 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Tue, 25 Jan 2022 09:55:54 +0100
Subject: [PATCH 23/23] Fix uncommenting solving function.

Now x is actually solved with the help of the equation system.
---
 src/BTCSDiffusion.cpp | 19 +++++++++++--------
 src/main.cpp          |  2 +-
 2 files changed, 12 insertions(+), 9 deletions(-)

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 6c9d896..88c816f 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -112,11 +112,13 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
   A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3));
 
   A_matrix.insert(0, 0) = 1;
-  b_vector[0] = (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
+  b_vector[0] =
+      (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
 
   A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1;
   b_vector[size + bc_offset - 1] =
-    (right_is_constant ? right.value : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
+      (right_is_constant ? right.value
+                         : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
 
   // A_matrix.insert(0, 0) = 1;
   // A_matrix.insert(size + 1, size + 1) = 1;
@@ -131,7 +133,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
     b_vector[i] = -c[i + !(left_is_constant)];
   }
 
-std::cout << b_vector << "\n" << A_matrix << std::endl;
+  std::cout << b_vector << "\n" << A_matrix << std::endl;
 
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       solver;
@@ -141,7 +143,7 @@ std::cout << b_vector << "\n" << A_matrix << std::endl;
 
   std::cout << solver.lastErrorMessage() << std::endl;
 
-  // x_vector = solver.solve(b_vector);
+  x_vector = solver.solve(b_vector);
 
   std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
 
@@ -165,16 +167,17 @@ void BTCSDiffusion::simulate(std::vector<double> &c,
   }
 }
 
-inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc, double neighbor_c,
-                                    double neighbor_alpha) {
+inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc,
+                                           double neighbor_c,
+                                           double neighbor_alpha) {
 
   double val;
 
   if (bc.type == BTCSDiffusion::BC_CLOSED) {
     val = neighbor_c;
   } else if (bc.type == BTCSDiffusion::BC_FLUX) {
-    //TODO
-    // val = bc[index].value;
+    // TODO
+    //  val = bc[index].value;
   } else {
     // TODO: implement error handling here. Type was set to wrong value.
   }
diff --git a/src/main.cpp b/src/main.cpp
index 31445a2..c4549ee 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -33,7 +33,7 @@ int main(int argc, char *argv[]) {
 
   // loop 100 times
   // output is currently generated by the method itself
-  for (int i = 0; i < 1; i++) {
+  for (int i = 0; i < 100; i++) {
     diffu.simulate(field, alpha);
   }