From f4253f2e6ab3fe37f14bd76dde206ae75bb3f874 Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Tue, 8 Feb 2022 11:22:18 +0100
Subject: [PATCH 1/2] Replace internal handling of vectors by Eigen library
 members

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

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 92855dc..870479c 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -2,6 +2,8 @@
 
 #include <Eigen/SparseLU>
 
+#include <Eigen/src/Core/Map.h>
+#include <Eigen/src/Core/Matrix.h>
 #include <algorithm>
 #include <cassert>
 #include <iomanip>
@@ -68,7 +70,7 @@ void BTCSDiffusion::updateInternals() {
   bc.resize(cells, {BTCSDiffusion::BC_CLOSED, 0});
 }
 
-void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
+void BTCSDiffusion::simulate1D(Eigen::Map<Eigen::VectorXd> &c, boundary_condition left,
                                boundary_condition right,
                                const std::vector<double> &alpha, double dx,
                                int size) {
@@ -151,7 +153,9 @@ void BTCSDiffusion::setTimestep(double time_step) {
 void BTCSDiffusion::simulate(std::vector<double> &c,
                              const std::vector<double> &alpha) {
   if (this->grid_dim == 1) {
-    simulate1D(c, bc[0], bc[grid_cells[0] + 1], alpha, this->deltas[0],
+    assert(c.size() == grid_cells[0]);
+    Eigen::Map<Eigen::VectorXd> c_in(c.data(), this->grid_cells[0]);
+    simulate1D(c_in, bc[0], bc[grid_cells[0] + 1], alpha, this->deltas[0],
                this->grid_cells[0]);
   }
 }
diff --git a/src/BTCSDiffusion.hpp b/src/BTCSDiffusion.hpp
index ea3473d..d44a407 100644
--- a/src/BTCSDiffusion.hpp
+++ b/src/BTCSDiffusion.hpp
@@ -2,6 +2,8 @@
 #define BTCSDIFFUSION_H_
 
 #include <Eigen/Sparse>
+#include <Eigen/src/Core/Map.h>
+#include <Eigen/src/Core/Matrix.h>
 #include <tuple>
 #include <type_traits>
 #include <vector>
@@ -135,7 +137,7 @@ private:
   } boundary_condition;
   typedef Eigen::Triplet<double> T;
 
-  void simulate1D(std::vector<double> &c, boundary_condition left,
+  void simulate1D(Eigen::Map<Eigen::VectorXd> &c, boundary_condition left,
                   boundary_condition right, const std::vector<double> &alpha,
                   double dx, int size);
   void simulate2D(std::vector<double> &c);

From 8de34ad65b1213e45247ccc756121f271b9c4cad Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Tue, 8 Feb 2022 11:26:54 +0100
Subject: [PATCH 2/2] Replace copying of output vector by oneliner

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

diff --git a/src/BTCSDiffusion.cpp b/src/BTCSDiffusion.cpp
index 870479c..ceef611 100644
--- a/src/BTCSDiffusion.cpp
+++ b/src/BTCSDiffusion.cpp
@@ -70,7 +70,8 @@ void BTCSDiffusion::updateInternals() {
   bc.resize(cells, {BTCSDiffusion::BC_CLOSED, 0});
 }
 
-void BTCSDiffusion::simulate1D(Eigen::Map<Eigen::VectorXd> &c, boundary_condition left,
+void BTCSDiffusion::simulate1D(Eigen::Map<Eigen::VectorXd> &c,
+                               boundary_condition left,
                                boundary_condition right,
                                const std::vector<double> &alpha, double dx,
                                int size) {
@@ -78,10 +79,10 @@ void BTCSDiffusion::simulate1D(Eigen::Map<Eigen::VectorXd> &c, boundary_conditio
   bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT);
   bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT);
 
-  //The sizes for matrix and vectors of the equation system is defined by the
-  //actual size of the input vector and if the system is (partially) closed.
-  //Then we will need ghost nodes. So this variable will give the count of ghost
-  //nodes.
+  // The sizes for matrix and vectors of the equation system is defined by the
+  // actual size of the input vector and if the system is (partially) closed.
+  // Then we will need ghost nodes. So this variable will give the count of
+  // ghost nodes.
   int bc_offset = !left_is_constant + !right_is_constant;
   ;
 
@@ -140,10 +141,8 @@ void BTCSDiffusion::simulate1D(Eigen::Map<Eigen::VectorXd> &c, boundary_conditio
 
   x_vector = solver.solve(b_vector);
 
-  //fill solution back in place into =c= vector
-  for (int i = 0, j = i + !left_is_constant; i < c.size(); i++, j++) {
-    c[i] = x_vector[i + !left_is_constant];
-  }
+  // write back result to input/output vector
+  c = x_vector.segment(!left_is_constant, c.size());
 }
 
 void BTCSDiffusion::setTimestep(double time_step) {