From 04763406cbc3fe77d1fff10a8584a5f2ff08cbfd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Max=20L=C3=BCbke?= <mluebke@gfz-potsdam.de>
Date: Tue, 26 Oct 2021 11:54:08 +0200
Subject: [PATCH] udpate header + function

---
 src/diffusion.cpp | 25 ++++++++++++++++++++++++-
 src/diffusion.hpp |  4 ++--
 2 files changed, 26 insertions(+), 3 deletions(-)

diff --git a/src/diffusion.cpp b/src/diffusion.cpp
index e668388..11410ff 100644
--- a/src/diffusion.cpp
+++ b/src/diffusion.cpp
@@ -2,6 +2,13 @@
 
 #include <Eigen/SparseLU>
 #include <Eigen/src/Core/Matrix.h>
+#include <Eigen/src/Core/util/Constants.h>
+#include <Eigen/src/OrderingMethods/Ordering.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 <iostream>
 
 void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
             double timestep) {
@@ -11,7 +18,8 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
 
   int size = x * y;
 
-  Eigen::VectorXd d = Eigen::VectorXd::Constant(size, 0);
+  Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
+  Eigen::VectorXd x_out(size);
   std::vector<T> tripletList;
   tripletList.reserve(((x - 1) * (y - 1) * 5));
 
@@ -25,6 +33,21 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
       tripletList.push_back(T((i * x) + j, ((i * x) + j) - x, sy));
       tripletList.push_back(T((i * x) + j, ((i * x) + j) + 1, sx));
       tripletList.push_back(T((i * x) + j, ((i * x) + j) - 1, sx));
+
+      b[(i * x) + j] = -c[(i * x) + j];
     }
   }
+
+  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.factorize(A);
+
+  std::cout << A << std::endl;
+
+  x_out = solver.solve(b);
 }
diff --git a/src/diffusion.hpp b/src/diffusion.hpp
index 5e04332..6645f8c 100644
--- a/src/diffusion.hpp
+++ b/src/diffusion.hpp
@@ -6,6 +6,6 @@
 
 typedef Eigen::Triplet<double> T;
 
-extern void BTCS2D(int x, int y, std::vector<double> &c, double alpha,
-                   double timestep);
+extern void BTCS2D(int x, int y, std::vector<double> &c,
+                   std::vector<double> &alpha, double timestep);
 #endif // DIFFUSION_H_