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] 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;