use private datatypes to adress solver matrix

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>
@@ -32,8 +26,8 @@ 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
+  // // per default use Neumann condition with gradient of 0 at the end of the
-  // std::fill(this->bc.begin(), this->bc.end(), -1);
+  // 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,
+void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
-                             double timestep) {
+                               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);
+  // set sizes of private and yet allocated vectors
-  Eigen::VectorXd x_out(size);
+  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;
+  // std::vector<T> tripletList;
-  tripletList.reserve(c.size() * 3 + bc.size());
+  // 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 all concentrations create one row in matrix A
-  for (int i = 1; i < this->dim_x + 1; i++) {
+  // for (int i = 1; i < this->dim_x + 1; i++) {
-    double sx = (alpha[i - 1] * timestep) / (dx * dx);
+  //   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];
+  //   b[A_line] = -c[i - 1];
-    A_line++;
+  //   A_line++;
-  }
+  // }
 
-  // append left and right boundary conditions/ghost zones
+  // // append left and right boundary conditions/ghost zones
-  tripletList.push_back(T(A_line, 0, 1));
+  // tripletList.push_back(T(A_line, 0, 1));
 
-  // if value is -1 apply Neumann condition with given gradient
+  // // if value is -1 apply Neumann condition with given gradient
-  // TODO: set specific gradient
+  // // TODO: set specific gradient
-  if (bc[0] == -1)
+  // if (bc[0] == -1)
-    b[A_line] = c[0];
+  //   b[A_line] = c[0];
-  // else apply given Dirichlet condition
+  // // else apply given Dirichlet condition
-  else
+  // else
-    b[A_line] = this->bc[0];
+  //   b[A_line] = this->bc[0];
 
-  A_line++;
+  // A_line++;
-  tripletList.push_back(T(A_line, size - 1, 1));
+  // tripletList.push_back(T(A_line, size - 1, 1));
-  // b[A_line] = bc[1];
+  // // b[A_line] = bc[1];
-  if (bc[1] == -1)
+  // if (bc[1] == -1)
-    b[A_line] = c[c.size() - 1];
+  //   b[A_line] = c[c.size() - 1];
-  else
+  // else
-    b[A_line] = this->bc[1];
+  //   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_matrix.resize(size, size);
-  A.setFromTriplets(tripletList.begin(), tripletList.end());
+  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;
+}

src/BTCSDiffusion.hpp

@@ -2,6 +2,7 @@
 #define BTCSDIFFUSION_H_
 
 #include <Eigen/Sparse>
+#include <Eigen/src/SparseCore/SparseMatrixBase.h>
 #include <tuple>
 #include <vector>
 
@@ -20,7 +21,6 @@ typedef std::vector<std::tuple<int,double>> boundary_condition;
 class BTCSDiffusion {
 
 public:
-
   static const int BC_NEUMANN;
   static const int BC_DIRICHLET;
 
@@ -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,
+  void simulate(std::vector<double> &c, std::vector<double> &alpha);
-                double timestep);
+
+  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;

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;