feat: rewrite library as template library

examples/BTCS_1D_proto_example.cpp

@@ -10,7 +10,7 @@ int main(int argc, char *argv[]) {
 
   // create a linear grid with 20 cells
   int cells = 20;
-  Grid grid = Grid(cells);
+  Grid64 grid(cells);
 
   MatrixXd concentrations = MatrixXd::Constant(1, 20, 0);
   concentrations(0, 0) = 2000;

examples/BTCS_2D_proto_example.cpp

@@ -13,7 +13,7 @@ int main(int argc, char *argv[]) {
   // create a grid with a 20 x 20 field
   int row = 40;
   int col = 50;
-  Grid grid = Grid(row, col);
+  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   // grid.setDomain(20, 20);

examples/CRNI_2D_proto_example.cpp

@@ -6,7 +6,7 @@ using namespace Eigen;
 int main(int argc, char *argv[]) {
   int row = 20;
   int col = 20;
-  Grid grid(row, col);
+  Grid64 grid(row, col);
 
   MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
   concentrations(10, 10) = 2000;

examples/FTCS_1D_proto_example.cpp

@@ -11,7 +11,7 @@ int main(int argc, char *argv[]) {
 
   // create a linear grid with 20 cells
   int cells = 20;
-  Grid grid = Grid(cells);
+  Grid64 grid(cells);
 
   MatrixXd concentrations = MatrixXd::Constant(1, 20, 20);
   grid.setConcentrations(concentrations);

examples/FTCS_2D_proto_closed_mdl.cpp

@@ -30,7 +30,7 @@ int main(int argc, char *argv[]) {
 
   // create a grid with a 20 x 20 field
   int n2 = row / 2 - 1;
-  Grid grid = Grid(row, col);
+  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   // grid.setDomain(20, 20);

examples/FTCS_2D_proto_example.cpp

@@ -23,7 +23,7 @@ int main(int argc, char *argv[]) {
   // create a grid with a 20 x 20 field
   int row = 20;
   int col = 20;
-  Grid grid = Grid(row, col);
+  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   // grid.setDomain(20, 20);

examples/FTCS_2D_proto_example_mdl.cpp

@@ -21,7 +21,7 @@ int main(int argc, char *argv[]) {
   int row = 64;
   int col = 64;
   int n2 = row / 2 - 1;
-  Grid grid = Grid(row, col);
+  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   // grid.setDomain(20, 20);

examples/profiling_openmp.cpp

@@ -29,7 +29,7 @@ int main(int argc, char *argv[]) {
         // myfile << "Iterations: " << iterations[j] << endl;
         for (int k = 0; k < repetition; k++) {
           cout << "Wiederholung: " << k << endl;
-          Grid grid = Grid(n[i], n[i]);
+          Grid64 grid(n[i], n[i]);
           grid.setDomain(1, 1);
 
           MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);

examples/profiling_speedup.cpp

@@ -29,7 +29,7 @@ int main(int argc, char *argv[]) {
         // myfile << "Iterations: " << iterations[j] << endl;
         for (int k = 0; k < repetition; k++) {
           cout << "Wiederholung: " << k << endl;
-          Grid grid = Grid(n[i], n[i]);
+          Grid64 grid(n[i], n[i]);
           grid.setDomain(1, 1);
 
           MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);

examples/reference-FTCS_2D_closed.cpp

@@ -12,7 +12,7 @@ int main(int argc, char *argv[]) {
   int domain_col = 10;
 
   // Grid
-  Grid grid = Grid(row, col);
+  Grid64 grid(row, col);
   grid.setDomain(domain_row, domain_col);
 
   MatrixXd concentrations = MatrixXd::Constant(row, col, 0);

include/tug/Boundary.hpp

@@ -9,6 +9,8 @@
 
 #include "Grid.hpp"
 #include <cstddef>
+#include <cstdint>
+#include <exception>
 
 /**
  * @brief Enum defining the two implemented boundary conditions.
@@ -27,7 +29,7 @@ enum BC_SIDE { BC_SIDE_LEFT, BC_SIDE_RIGHT, BC_SIDE_TOP, BC_SIDE_BOTTOM };
  * These can be flexibly used and combined later in other classes.
  * The class serves as an auxiliary class for structuring the Boundary class.
  */
-class BoundaryElement {
+template <class T> class BoundaryElement {
 public:
   /**
    * @brief Construct a new Boundary Element object for the closed case.
@@ -35,7 +37,7 @@ public:
    *        BC_TYPE_CLOSED, where the value takes -1 and does not hold any
    *        physical meaning.
    */
-  BoundaryElement();
+  BoundaryElement(){};
 
   /**
    * @brief Construct a new Boundary Element object for the constant case.
@@ -45,7 +47,7 @@ public:
    * @param value Value of the constant concentration to be assumed at the
    *              corresponding boundary element.
    */
-  BoundaryElement(double value);
+  BoundaryElement(T _value) : value(_value), type(BC_TYPE_CONSTANT) {}
 
   /**
    * @brief Allows changing the boundary type of a corresponding
@@ -54,7 +56,7 @@ public:
    * @param type Type of boundary condition. Either BC_TYPE_CONSTANT or
                  BC_TYPE_CLOSED.
    */
-  void setType(BC_TYPE type);
+  void setType(BC_TYPE type) { this->type = type; };
 
   /**
    * @brief Sets the value of a boundary condition for the constant case.
@@ -62,7 +64,17 @@ public:
    * @param value Concentration to be considered constant for the
    *              corresponding boundary element.
    */
-  void setValue(double value);
+  void setValue(double value) {
+    if (value < 0) {
+      throw std::invalid_argument("No negative concentration allowed.");
+    }
+    if (type == BC_TYPE_CLOSED) {
+      throw std::invalid_argument(
+          "No constant boundary concentrations can be set for closed "
+          "boundaries. Please change type first.");
+    }
+    this->value = value;
+  }
 
   /**
    * @brief Return the type of the boundary condition, i.e. whether the
@@ -71,25 +83,25 @@ public:
    * @return BC_TYPE Type of boundary condition, either BC_TYPE_CLOSED or
              BC_TYPE_CONSTANT.
    */
-  BC_TYPE getType();
+  BC_TYPE getType() const { return this->type; }
 
   /**
    * @brief Return the concentration value for the constant boundary condition.
    *
    * @return double Value of the concentration.
    */
-  double getValue();
+  T getValue() const { return this->value; }
 
 private:
   BC_TYPE type{BC_TYPE_CLOSED};
-  double value{-1};
+  T value{-1};
 };
 
 /**
  * This class implements the functionality and management of the boundary
  * conditions in the grid to be simulated.
  */
-class Boundary {
+template <class T> class Boundary {
 public:
   /**
    * @brief Creates a boundary object based on the passed grid object and
@@ -98,7 +110,27 @@ public:
    * @param grid Grid object on the basis of which the simulation takes place
    *             and from which the dimensions (in 2D case) are taken.
    */
-  Boundary(Grid grid);
+  Boundary(const Grid<T> &grid)
+      : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) {
+    if (this->dim == 1) {
+      this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(
+          2); // in 1D only left and right boundary
+
+      this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement<T>());
+      this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement<T>());
+    } else if (this->dim == 2) {
+      this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(4);
+
+      this->boundaries[BC_SIDE_LEFT] =
+          std::vector<BoundaryElement<T>>(this->rows, BoundaryElement<T>());
+      this->boundaries[BC_SIDE_RIGHT] =
+          std::vector<BoundaryElement<T>>(this->rows, BoundaryElement<T>());
+      this->boundaries[BC_SIDE_TOP] =
+          std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
+      this->boundaries[BC_SIDE_BOTTOM] =
+          std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
+    }
+  }
 
   /**
    * @brief Sets all elements of the specified boundary side to the boundary
@@ -106,7 +138,21 @@ public:
    *
    * @param side Side to be set to closed, e.g. BC_SIDE_LEFT.
    */
-  void setBoundarySideClosed(BC_SIDE side);
+  void setBoundarySideClosed(BC_SIDE side) {
+    if (this->dim == 1) {
+      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
+        throw std::invalid_argument(
+            "For the one-dimensional case, only the BC_SIDE_LEFT and "
+            "BC_SIDE_RIGHT borders exist.");
+      }
+    }
+
+    const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
+    const int n = is_vertical ? this->rows : this->cols;
+
+    this->boundaries[side] =
+        std::vector<BoundaryElement<T>>(n, BoundaryElement<T>());
+  }
 
   /**
    * @brief Sets all elements of the specified boundary side to the boundary
@@ -117,7 +163,21 @@ public:
    * @param value Concentration to be set for all elements of the specified
    * page.
    */
-  void setBoundarySideConstant(BC_SIDE side, double value);
+  void setBoundarySideConstant(BC_SIDE side, double value) {
+    if (this->dim == 1) {
+      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
+        throw std::invalid_argument(
+            "For the one-dimensional case, only the BC_SIDE_LEFT and "
+            "BC_SIDE_RIGHT borders exist.");
+      }
+    }
+
+    const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
+    const int n = is_vertical ? this->rows : this->cols;
+
+    this->boundaries[side] =
+        std::vector<BoundaryElement<T>>(n, BoundaryElement<T>(value));
+  }
 
   /**
    * @brief Specifically sets the boundary element of the specified side
@@ -128,7 +188,13 @@ public:
    *              boundary side. Must index an element of the corresponding
    * side.
    */
-  void setBoundaryElementClosed(BC_SIDE side, int index);
+  void setBoundaryElemenClosed(BC_SIDE side, int index) {
+    // tests whether the index really points to an element of the boundary side.
+    if ((boundaries[side].size() < index) || index < 0) {
+      throw std::invalid_argument("Index is selected either too large or too small.");
+    }
+    this->boundaries[side][index].setType(BC_TYPE_CLOSED);
+  }
 
   /**
    * @brief Specifically sets the boundary element of the specified side
@@ -142,7 +208,14 @@ public:
    * @param value Concentration value to which the boundary element should be
    set.
    */
-  void setBoundaryElementConstant(BC_SIDE side, int index, double value);
+  void setBoundaryElementConstant(BC_SIDE side, int index, double value) {
+    // tests whether the index really points to an element of the boundary side.
+    if ((boundaries[side].size() < index) || index < 0) {
+      throw std::invalid_argument("Index is selected either too large or too small.");
+    }
+    this->boundaries[side][index].setType(BC_TYPE_CONSTANT);
+    this->boundaries[side][index].setValue(value);
+  }
 
   /**
    * @brief Returns the boundary condition of a specified side as a vector
@@ -150,19 +223,41 @@ public:
    *
    * @param side Boundary side from which the boundary conditions are to be
    * returned.
-   * @return vector<BoundaryElement> Contains the boundary conditions as
+   * @return vector<BoundaryElement<T>> Contains the boundary conditions as
-   * BoundaryElement objects.
+   * BoundaryElement<T> objects.
    */
-  const std::vector<BoundaryElement> getBoundarySide(BC_SIDE side);
+  const std::vector<BoundaryElement<T>> &getBoundarySide(BC_SIDE side) const {
+    if (this->dim == 1) {
+      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
+        throw std::invalid_argument(
+            "For the one-dimensional trap, only the BC_SIDE_LEFT and "
+            "BC_SIDE_RIGHT borders exist.");
+      }
+    }
+    return this->boundaries[side];
+  }
 
   /**
    * @brief Get thes Boundary Side Values as a vector. Value is -1 in case some
    specific boundary is closed.
    *
    * @param side Boundary side for which the values are to be returned.
-   * @return VectorXd Vector with values as doubles.
+   * @return VectorX<T> Vector with values as T.
    */
-  Eigen::VectorXd getBoundarySideValues(BC_SIDE side);
+  Eigen::VectorX<T> getBoundarySideValues(BC_SIDE side) const {
+    const std::size_t length = boundaries[side].size();
+    Eigen::VectorX<T> values(length);
+
+    for (int i = 0; i < length; i++) {
+      if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) {
+        values(i) = -1;
+        continue;
+      }
+      values(i) = getBoundaryElementValue(side, i);
+    }
+
+    return values;
+  }
 
   /**
    * @brief Returns the boundary condition of a specified element on a given
@@ -172,9 +267,15 @@ public:
    * @param index Index of the boundary element on the corresponding
    *              boundary side. Must index an element of the corresponding
    * side.
-   * @return BoundaryElement Boundary condition as a BoundaryElement object.
+   * @return BoundaryElement<T> Boundary condition as a BoundaryElement<T>
+   * object.
    */
-  BoundaryElement getBoundaryElement(BC_SIDE side, int index);
+  BoundaryElement<T> getBoundaryElement(BC_SIDE side, int index) const {
+    if ((boundaries[side].size() < index) || index < 0) {
+      throw std::invalid_argument("Index is selected either too large or too small.");
+    }
+    return this->boundaries[side][index];
+  }
 
   /**
    * @brief Returns the type of a boundary condition, i.e. either BC_TYPE_CLOSED
@@ -186,25 +287,42 @@ public:
    side.
    * @return BC_TYPE Boundary Type of the corresponding boundary condition.
    */
-  BC_TYPE getBoundaryElementType(BC_SIDE side, int index);
+  BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const {
+    if ((boundaries[side].size() < index) || index < 0) {
+      throw std::invalid_argument("Index is selected either too large or too small.");
+    }
+    return this->boundaries[side][index].getType();
+  }
 
   /**
    * @brief Returns the concentration value of a corresponding
-   *        BoundaryElement object if it is a constant boundary condition.
+   *        BoundaryElement<T> object if it is a constant boundary condition.
    *
    * @param side Boundary side in which the boundary condition value is
    *             located.
    * @param index Index of the boundary element on the corresponding
    *              boundary side. Must index an element of the corresponding
    *              side.
-   * @return double Concentration of the corresponding BoundaryElement object.
+   * @return double Concentration of the corresponding BoundaryElement<T>
+   * object.
    */
-  double getBoundaryElementValue(BC_SIDE side, int index);
+  T getBoundaryElementValue(BC_SIDE side, int index) const {
+    if ((boundaries[side].size() < index) || index < 0) {
+      throw std::invalid_argument("Index is selected either too large or too small.");
+    }
+    if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) {
+      throw std::invalid_argument(
+          "A value can only be output if it is a constant boundary condition.");
+    }
+    return this->boundaries[side][index].getValue();
+  }
 
 private:
-  Grid grid; // Boundary is directly dependent on the dimensions of a predefined
+  const std::uint8_t dim;
+  const std::uint32_t cols;
+  const std::uint32_t rows;
 
-  std::vector<std::vector<BoundaryElement>>
+  std::vector<std::vector<BoundaryElement<T>>>
       boundaries; // Vector with Boundary Element information
 };
 

include/tug/Grid.hpp

@@ -10,8 +10,9 @@
 
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <stdexcept>
 
-class Grid {
+template <class T> class Grid {
 public:
   /**
    * @brief Constructs a new 1D-Grid object of a given length, which holds a
@@ -21,7 +22,18 @@ public:
    *
    * @param length Length of the 1D-Grid. Must be greater than 3.
    */
-  Grid(int length);
+  Grid(int length) : col(length), domainCol(length) {
+    if (length <= 3) {
+      throw std::invalid_argument(
+          "Given grid length too small. Must be greater than 3.");
+    }
+
+    this->dim = 1;
+    this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
+
+    this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL);
+    this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL);
+  }
 
   /**
    * @brief Constructs a new 2D-Grid object of given dimensions, which holds a
@@ -34,146 +46,280 @@ public:
    * @param row Length of the 2D-Grid in y-direction. Must be greater than 3.
    * @param col Length of the 2D-Grid in x-direction. Must be greater than 3.
    */
-  Grid(int row, int col);
+  Grid(int _row, int _col)
+      : row(_row), col(_col), domainRow(_row), domainCol(_col) {
+    if (row <= 3 || col <= 3) {
+      throw std::invalid_argument(
+          "Given grid dimensions too small. Must each be greater than 3.");
+    }
+
+    this->dim = 2;
+    this->deltaRow = double(this->domainRow) / double(this->row); // -> 1
+    this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
+
+    this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
+    this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
+    this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
+  }
 
   /**
    * @brief Sets the concentrations matrix for a 1D or 2D-Grid.
    *
-   * @param concentrations An Eigen3 MatrixXd holding the concentrations. Matrix
+   * @param concentrations An Eigen3 MatrixX<T> holding the concentrations.
-   * must have correct dimensions as defined in row and col. (Or length, in 1D
+   * Matrix must have correct dimensions as defined in row and col. (Or length,
-   * case).
+   * in 1D case).
    */
-  void setConcentrations(const Eigen::MatrixXd &concentrations);
+  void setConcentrations(const Eigen::MatrixX<T> &concentrations) {
+    if (concentrations.rows() != this->row ||
+        concentrations.cols() != this->col) {
+      throw std::invalid_argument(
+          "Given matrix of concentrations mismatch with Grid dimensions!");
+    }
+
+    this->concentrations = concentrations;
+  }
 
   /**
    * @brief Gets the concentrations matrix for a Grid.
    *
-   * @return MatrixXd An Eigen3 matrix holding the concentrations and having the
+   * @return MatrixX<T> An Eigen3 matrix holding the concentrations and having
-   *                  same dimensions as the grid.
+   * the same dimensions as the grid.
    */
-  const Eigen::MatrixXd &getConcentrations() { return this->concentrations; }
+  const Eigen::MatrixX<T> &getConcentrations() { return this->concentrations; }
 
   /**
    * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
    * dimensional.
    *
-   * @param alpha An Eigen3 MatrixXd with 1 row holding the alpha coefficients.
+   * @param alpha An Eigen3 MatrixX<T> with 1 row holding the alpha
-   * Matrix columns must have same size as length of grid.
+   * coefficients. Matrix columns must have same size as length of grid.
    */
-  void setAlpha(const Eigen::MatrixXd &alpha);
+  void setAlpha(const Eigen::MatrixX<T> &alpha) {
+    if (dim != 1) {
+      throw std::invalid_argument(
+          "Grid is not one dimensional, you should probably "
+          "use 2D setter function!");
+    }
+    if (alpha.rows() != 1 || alpha.cols() != this->col) {
+      throw std::invalid_argument(
+          "Given matrix of alpha coefficients mismatch with Grid dimensions!");
+    }
+
+    this->alphaX = alpha;
+  }
 
   /**
    * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
    * dimensional.
    *
-   * @param alphaX An Eigen3 MatrixXd holding the alpha coefficients in
+   * @param alphaX An Eigen3 MatrixX<T> holding the alpha coefficients in
    * x-direction. Matrix must be of same size as the grid.
-   * @param alphaY An Eigen3 MatrixXd holding the alpha coefficients in
+   * @param alphaY An Eigen3 MatrixX<T> holding the alpha coefficients in
    * y-direction. Matrix must be of same size as the grid.
    */
-  void setAlpha(const Eigen::MatrixXd &alphaX, const Eigen::MatrixXd &alphaY);
+  void setAlpha(const Eigen::MatrixX<T> &alphaX,
+                const Eigen::MatrixX<T> &alphaY) {
+    if (dim != 2) {
+      throw std::invalid_argument(
+          "Grid is not two dimensional, you should probably "
+          "use 1D setter function!");
+    }
+    if (alphaX.rows() != this->row || alphaX.cols() != this->col) {
+      throw std::invalid_argument(
+          "Given matrix of alpha coefficients in x-direction "
+          "mismatch with GRid dimensions!");
+    }
+    if (alphaY.rows() != this->row || alphaY.cols() != this->col) {
+      throw std::invalid_argument(
+          "Given matrix of alpha coefficients in y-direction "
+          "mismatch with GRid dimensions!");
+    }
+
+    this->alphaX = alphaX;
+    this->alphaY = alphaY;
+  }
 
   /**
    * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
    * dimensional.
    *
-   * @return MatrixXd A matrix with 1 row holding the alpha coefficients.
+   * @return MatrixX<T> A matrix with 1 row holding the alpha coefficients.
    */
-  const Eigen::MatrixXd &getAlpha();
+  const Eigen::MatrixX<T> &getAlpha() const {
+    if (dim != 1) {
+      throw std::invalid_argument(
+          "Grid is not one dimensional, you should probably "
+          "use either getAlphaX() or getAlphaY()!");
+    }
+
+    return this->alphaX;
+  }
 
   /**
    * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid.
    * Grid must be two dimensional.
    *
-   * @return MatrixXd A matrix holding the alpha coefficients in x-direction.
+   * @return MatrixX<T> A matrix holding the alpha coefficients in x-direction.
    */
-  const Eigen::MatrixXd &getAlphaX();
+  const Eigen::MatrixX<T> &getAlphaX() const {
+
+    if (dim != 2) {
+      throw std::invalid_argument(
+          "Grid is not two dimensional, you should probably use getAlpha()!");
+    }
+
+    return this->alphaX;
+  }
 
   /**
    * @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
    * Grid must be two dimensional.
    *
-   * @return MatrixXd A matrix holding the alpha coefficients in y-direction.
+   * @return MatrixX<T> A matrix holding the alpha coefficients in y-direction.
    */
-  const Eigen::MatrixXd &getAlphaY();
+  const Eigen::MatrixX<T> &getAlphaY() const {
+
+    if (dim != 2) {
+      throw std::invalid_argument(
+          "Grid is not two dimensional, you should probably use getAlpha()!");
+    }
+
+    return this->alphaY;
+  }
 
   /**
    * @brief Gets the dimensions of the grid.
    *
    * @return int Dimensions, either 1 or 2.
    */
-  int getDim();
+  int getDim() const { return this->dim; }
 
   /**
    * @brief Gets length of 1D grid. Must be one dimensional grid.
    *
    * @return int Length of 1D grid.
    */
-  int getLength();
+  int getLength() const {
+    if (dim != 1) {
+      throw std::invalid_argument(
+          "Grid is not one dimensional, you should probably "
+          "use getRow() or getCol()!");
+    }
+
+    return col;
+  }
 
   /**
    * @brief Gets the number of rows of the grid.
    *
    * @return int Number of rows.
    */
-  int getRow();
+  int getRow() const { return this->row; }
 
   /**
    * @brief Gets the number of columns of the grid.
    *
    * @return int Number of columns.
    */
-  int getCol();
+  int getCol() const { return this->col; }
 
   /**
    * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional.
    *
    * @param domainLength A double value of the domain length. Must be positive.
    */
-  void setDomain(double domainLength);
+  void setDomain(double domainLength) {
+    if (dim != 1) {
+      throw std::invalid_argument(
+          "Grid is not one dimensional, you should probaly "
+          "use the 2D domain setter!");
+    }
+    if (domainLength <= 0) {
+      throw std::invalid_argument("Given domain length is not positive!");
+    }
+
+    this->domainCol = domainLength;
+    this->deltaCol = double(this->domainCol) / double(this->col);
+  }
 
   /**
    * @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional.
    *
-   * @param domainRow A double value of the domain size in y-direction. Must be
+   * @param domainRow A double value of the domain size in y-direction. Must
-   * positive.
+   * be positive.
-   * @param domainCol A double value of the domain size in x-direction. Must be
+   * @param domainCol A double value of the domain size in x-direction. Must
-   * positive.
+   * be positive.
    */
-  void setDomain(double domainRow, double domainCol);
+  void setDomain(double domainRow, double domainCol) {
+    if (dim != 2) {
+      throw std::invalid_argument(
+          "Grid is not two dimensional, you should probably "
+          "use the 1D domain setter!");
+    }
+    if (domainRow <= 0 || domainCol <= 0) {
+      throw std::invalid_argument("Given domain size is not positive!");
+    }
+
+    this->domainRow = domainRow;
+    this->domainCol = domainCol;
+    this->deltaRow = double(this->domainRow) / double(this->row);
+    this->deltaCol = double(this->domainCol) / double(this->col);
+  }
 
   /**
    * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional.
    *
    * @return double Delta value.
    */
-  double getDelta();
+  T getDelta() const {
+
+    if (dim != 1) {
+      throw std::invalid_argument(
+          "Grid is not one dimensional, you should probably "
+          "use the 2D delta getters");
+    }
+
+    return this->deltaCol;
+  }
 
   /**
    * @brief Gets the delta value in x-direction.
    *
    * @return double Delta value in x-direction.
    */
-  double getDeltaCol();
+  T getDeltaCol() const { return this->deltaCol; }
 
   /**
    * @brief Gets the delta value in y-direction. Must be two dimensional grid.
    *
    * @return double Delta value in y-direction.
    */
-  double getDeltaRow();
+  T getDeltaRow() const {
+    if (dim != 2) {
+      throw std::invalid_argument(
+          "Grid is not two dimensional, meaning there is no "
+          "delta in y-direction!");
+    }
+
+    return this->deltaRow;
+  }
 
 private:
-  int col;                        // number of grid columns
+  int col;                          // number of grid columns
-  int row{1};                     // number of grid rows
+  int row{1};                       // number of grid rows
-  int dim;                        // 1D or 2D
+  int dim;                          // 1D or 2D
-  double domainCol;               // number of domain columns
+  T domainCol;                      // number of domain columns
-  double domainRow{0};            // number of domain rows
+  T domainRow{0};                   // number of domain rows
-  double deltaCol;                // delta in x-direction (between columns)
+  T deltaCol;                       // delta in x-direction (between columns)
-  double deltaRow{0};             // delta in y-direction (between rows)
+  T deltaRow{0};                    // delta in y-direction (between rows)
-  Eigen::MatrixXd concentrations; // Matrix holding grid concentrations
+  Eigen::MatrixX<T> concentrations; // Matrix holding grid concentrations
-  Eigen::MatrixXd alphaX; // Matrix holding alpha coefficients in x-direction
+  Eigen::MatrixX<T> alphaX; // Matrix holding alpha coefficients in x-direction
-  Eigen::MatrixXd alphaY; // Matrix holding alpha coefficients in y-direction
+  Eigen::MatrixX<T> alphaY; // Matrix holding alpha coefficients in y-direction
+
+  static constexpr double MAT_INIT_VAL = 0;
 };
 
+using Grid64 = Grid<double>;
+using Grid32 = Grid<float>;
+
 #endif // GRID_H_

include/tug/Simulation.hpp

@@ -11,6 +11,11 @@
 
 #include "Boundary.hpp"
 #include "Grid.hpp"
+#include <exception>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
@@ -77,7 +82,8 @@ enum TIME_MEASURE {
  * time step, number of simulations, etc.
  *
  */
-class Simulation {
+
+template <class T> class Simulation {
 public:
   /**
    * @brief Set up a simulation environment. The timestep and number of
@@ -91,7 +97,8 @@ public:
    * @param bc Valid boundary condition object
    * @param approach Approach to solving the problem. Either FTCS or BTCS.
    */
-  Simulation(Grid &grid, Boundary &bc, APPROACH approach);
+  Simulation(Grid<T> &_grid, Boundary<T> &_bc, APPROACH _approach)
+      : grid(_grid), bc(_bc), approach(_approach){};
 
   /**
    * @brief Set the option to output the results to a CSV file. Off by default.
@@ -107,7 +114,13 @@ public:
    *                     - CSV_OUTPUT_XTREME: produce csv output with all
    *                       concentration matrices and simulation environment
    */
-  void setOutputCSV(CSV_OUTPUT csv_output);
+  void setOutputCSV(CSV_OUTPUT csv_output) {
+    if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) {
+      throw std::invalid_argument("Invalid CSV output option given!");
+    }
+
+    this->csv_output = csv_output;
+  }
 
   /**
    * @brief Set the options for outputting information to the console. Off by
@@ -122,7 +135,14 @@ public:
    *                        - CONSOLE_OUTPUT_VERBOSE: print all concentration
    * matrices to console
    */
-  void setOutputConsole(CONSOLE_OUTPUT console_output);
+  void setOutputConsole(CONSOLE_OUTPUT console_output) {
+    if (console_output < CONSOLE_OUTPUT_OFF &&
+        console_output > CONSOLE_OUTPUT_VERBOSE) {
+      throw std::invalid_argument("Invalid console output option given!");
+    }
+
+    this->console_output = console_output;
+  }
 
   /**
    * @brief Set the Time Measure option. Off by default.
@@ -133,7 +153,13 @@ public:
    *                     - TIME_MEASURE_ON: Time of simulation run is printed to
    * console
    */
-  void setTimeMeasure(TIME_MEASURE time_measure);
+  void setTimeMeasure(TIME_MEASURE time_measure) {
+    if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) {
+      throw std::invalid_argument("Invalid time measure option given!");
+    }
+
+    this->time_measure = time_measure;
+  }
 
   /**
    * @brief Setting the time step for each iteration step. Time step must be
@@ -141,14 +167,14 @@ public:
    *
    * @param timestep Valid timestep greater than zero.
    */
-  void setTimestep(double timestep);
+  void setTimestep(T timestep);
 
   /**
    * @brief Currently set time step is returned.
    *
    * @return double timestep
    */
-  double getTimestep();
+  T getTimestep() const { return this->timestep; }
 
   /**
    * @brief Set the desired iterations to be calculated. A value greater
@@ -156,7 +182,13 @@ public:
    *
    * @param iterations Number of iterations to be simulated.
    */
-  void setIterations(int iterations);
+  void setIterations(int iterations) {
+    if (iterations <= 0) {
+      throw std::invalid_argument(
+          "Number of iterations must be greater than zero.");
+    }
+    this->iterations = iterations;
+  }
 
   /**
    * @brief Set the desired linear equation solver to be used for BTCS approach.
@@ -165,7 +197,17 @@ public:
    * @param solver Solver to be used. Default is Thomas Algorithm as it is more
    * efficient for tridiagonal Matrices.
    */
-  void setSolver(SOLVER solver);
+  void setSolver(SOLVER solver) {
+    if (this->approach == FTCS_APPROACH) {
+      std::cerr
+          << "Warning: Solver was set, but FTCS approach initialized. Setting "
+             "the solver "
+             "is thus without effect."
+          << std::endl;
+    }
+
+    this->solver = solver;
+  }
 
   /**
    * @brief Set the number of desired openMP Threads.
@@ -175,20 +217,32 @@ public:
    * maximum number of processors is set as the default case during Simulation
    * construction.
    */
-  void setNumberThreads(int num_threads);
+  void setNumberThreads(int num_threads) {
+    if (numThreads > 0 && numThreads <= omp_get_num_procs()) {
+      this->numThreads = numThreads;
+    } else {
+      int maxThreadNumber = omp_get_num_procs();
+      throw std::invalid_argument(
+          "Number of threads exceeds the number of processor cores (" +
+          std::to_string(maxThreadNumber) + ") or is less than 1.");
+    }
+  }
 
   /**
    * @brief Return the currently set iterations to be calculated.
    *
    * @return int Number of iterations.
    */
-  int getIterations();
+  int getIterations() const { return this->iterations; }
 
   /**
    * @brief Outputs the current concentrations of the grid on the console.
    *
    */
-  void printConcentrationsConsole();
+  inline void printConcentrationsConsole() const {
+    std::cout << grid.getConcentrations() << std::endl;
+    std::cout << std::endl;
+  }
 
   /**
    * @brief Creates a CSV file with a name containing the current simulation
@@ -199,7 +253,52 @@ public:
    *
    * @return string Filename with configured simulation parameters.
    */
-  std::string createCSVfile();
+  std::string createCSVfile() const {
+    std::ofstream file;
+    int appendIdent = 0;
+    std::string appendIdentString;
+
+    // string approachString = (approach == 0) ? "FTCS" : "BTCS";
+    const std::string &approachString = this->approach_names[approach];
+    std::string row = std::to_string(grid.getRow());
+    std::string col = std::to_string(grid.getCol());
+    std::string numIterations = std::to_string(iterations);
+
+    std::string filename =
+        approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
+
+    while (std::filesystem::exists(filename)) {
+      appendIdent += 1;
+      appendIdentString = std::to_string(appendIdent);
+      filename = approachString + "_" + row + "_" + col + "_" + numIterations +
+                 "-" + appendIdentString + ".csv";
+    }
+
+    file.open(filename);
+    if (!file) {
+      exit(1);
+    }
+
+    // adds lines at the beginning of verbose output csv that represent the
+    // boundary conditions and their values -1 in case of closed boundary
+    if (csv_output == CSV_OUTPUT_XTREME) {
+      Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "",
+                              " ");
+      file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row)
+           << std::endl; // boundary left
+      file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row)
+           << std::endl; // boundary right
+      file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row)
+           << std::endl; // boundary top
+      file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row)
+           << std::endl; // boundary bottom
+      file << std::endl << std::endl;
+    }
+
+    file.close();
+
+    return filename;
+  }
 
   /**
    * @brief Writes the currently calculated concentration values of the grid
@@ -208,7 +307,19 @@ public:
    * @param filename Name of the file to which the concentration values are
    *                 to be written.
    */
-  void printConcentrationsCSV(const std::string &filename);
+  void printConcentrationsCSV(const std::string &filename) const {
+    std::ofstream file;
+
+    file.open(filename, std::ios_base::app);
+    if (!file) {
+      exit(1);
+    }
+
+    Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols);
+    file << grid.getConcentrations().format(do_not_align) << std::endl;
+    file << std::endl << std::endl;
+    file.close();
+  }
 
   /**
    * @brief Method starts the simulation process with the previously set
@@ -217,7 +328,7 @@ public:
   void run();
 
 private:
-  double timestep{-1};
+  T timestep{-1};
   int iterations{-1};
   int innerIterations{1};
   int numThreads{omp_get_num_procs()};
@@ -225,8 +336,8 @@ private:
   CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF};
   TIME_MEASURE time_measure{TIME_MEASURE_OFF};
 
-  Grid &grid;
+  Grid<T> &grid;
-  Boundary &bc;
+  Boundary<T> &bc;
   APPROACH approach;
   SOLVER solver{THOMAS_ALGORITHM_SOLVER};
 

naaice/BTCS_2D_NAAICE.cpp

@@ -103,23 +103,23 @@ int main(int argc, char *argv[]) {
   // create a grid with a 5 x 10 field
   constexpr int row = 5;
   constexpr int col = 10;
-  Grid grid = Grid(row, col);
+  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   grid.setDomain(0.005, 0.01);
 
   const auto init_values_vec = CSVToVector<double>(INPUT_CONC_FILE);
-  MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col);
+  Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col);
   grid.setConcentrations(concentrations);
 
   const double sum_init = concentrations.sum();
 
   // // (optional) set alphas of the grid, e.g.:
   const auto alphax_vec = CSVToVector<double>(INPUT_ALPHAX_FILE);
-  MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col);
+  Eigen::MatrixXd alphax = rmVecTocmMatrix(alphax_vec, row, col);
 
   constexpr double alphay_val = 5e-10;
-  MatrixXd alphay = MatrixXd::Constant(row, col, alphay_val); // row,col,value
+  Eigen::MatrixXd alphay = Eigen::MatrixXd::Constant(row, col, alphay_val); // row,col,value
   grid.setAlpha(alphax, alphay);
 
   // // ******************

src/BTCS.cpp

@@ -45,13 +45,15 @@ constexpr std::pair<T, T> calcBoundaryCoeffClosed(T alpha_center, T alpha_side,
 }
 
 // creates coefficient matrix for next time step from alphas in x-direction
-static Eigen::SparseMatrix<double>
+template <class T>
-createCoeffMatrix(Eigen::MatrixXd &alpha, std::vector<BoundaryElement> &bcLeft,
+static Eigen::SparseMatrix<T>
-                  std::vector<BoundaryElement> &bcRight, int numCols,
+createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
-                  int rowIndex, double sx) {
+                  const std::vector<BoundaryElement<T>> &bcLeft,
+                  const std::vector<BoundaryElement<T>> &bcRight, int numCols,
+                  int rowIndex, T sx) {
 
   // square matrix of column^2 dimension for the coefficients
-  Eigen::SparseMatrix<double> cm(numCols, numCols);
+  Eigen::SparseMatrix<T> cm(numCols, numCols);
   cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
 
   // left column
@@ -142,14 +144,18 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc,
 
 // creates a solution vector for next time step from the current state of
 // concentrations
-static Eigen::VectorXd createSolutionVector(
+template <class T>
-    Eigen::MatrixXd &concentrations, Eigen::MatrixXd &alphaX,
+static Eigen::VectorX<T>
-    Eigen::MatrixXd &alphaY, std::vector<BoundaryElement> &bcLeft,
+createSolutionVector(const Eigen::MatrixX<T> &concentrations,
-    std::vector<BoundaryElement> &bcRight, std::vector<BoundaryElement> &bcTop,
+                     const Eigen::MatrixX<T> &alphaX,
-    std::vector<BoundaryElement> &bcBottom, int length, int rowIndex, double sx,
+                     const Eigen::MatrixX<T> &alphaY,
-    double sy) {
+                     const std::vector<BoundaryElement<T>> &bcLeft,
+                     const std::vector<BoundaryElement<T>> &bcRight,
+                     const std::vector<BoundaryElement<T>> &bcTop,
+                     const std::vector<BoundaryElement<T>> &bcBottom,
+                     int length, int rowIndex, T sx, T sy) {
 
-  Eigen::VectorXd sv(length);
+  Eigen::VectorX<T> sv(length);
   const std::size_t numRows = concentrations.rows();
 
   // inner rows
@@ -235,10 +241,11 @@ static Eigen::VectorXd createSolutionVector(
 // solver for linear equation system; A corresponds to coefficient matrix,
 // b to the solution vector
 // use of EigenLU solver
-static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix<double> &A,
+template <class T>
-                                        Eigen::VectorXd &b) {
+static Eigen::VectorX<T> EigenLUAlgorithm(Eigen::SparseMatrix<T> &A,
+                                          Eigen::VectorX<T> &b) {
 
-  Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
+  Eigen::SparseLU<Eigen::SparseMatrix<T>> solver;
   solver.analyzePattern(A);
   solver.factorize(A);
 
@@ -248,14 +255,15 @@ static Eigen::VectorXd EigenLUAlgorithm(Eigen::SparseMatrix<double> &A,
 // solver for linear equation system; A corresponds to coefficient matrix,
 // b to the solution vector
 // implementation of Thomas Algorithm
-static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix<double> &A,
+template <class T>
-                                       Eigen::VectorXd &b) {
+static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
+                                         Eigen::VectorX<T> &b) {
   Eigen::Index n = b.size();
 
-  Eigen::VectorXd a_diag(n);
+  Eigen::VectorX<T> a_diag(n);
-  Eigen::VectorXd b_diag(n);
+  Eigen::VectorX<T> b_diag(n);
-  Eigen::VectorXd c_diag(n);
+  Eigen::VectorX<T> c_diag(n);
-  Eigen::VectorXd x_vec = b;
+  Eigen::VectorX<T> x_vec = b;
 
   // Fill diagonals vectors
   b_diag[0] = A.coeff(0, 0);
@@ -291,23 +299,24 @@ static Eigen::VectorXd ThomasAlgorithm(Eigen::SparseMatrix<double> &A,
 }
 
 // BTCS solution for 1D grid
-static void
+template <class T>
-BTCS_1D(Grid &grid, Boundary &bc, double timestep,
+static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
-        Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
+                    Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
-                                      Eigen::VectorXd &b)) {
+                                                    Eigen::VectorX<T> &b)) {
   int length = grid.getLength();
-  double sx = timestep / (grid.getDelta() * grid.getDelta());
+  T sx = timestep / (grid.getDelta() * grid.getDelta());
 
-  Eigen::VectorXd concentrations_t1(length);
+  Eigen::VectorX<T> concentrations_t1(length);
 
-  Eigen::SparseMatrix<double> A;
+  Eigen::SparseMatrix<T> A;
-  Eigen::VectorXd b(length);
+  Eigen::VectorX<T> b(length);
 
-  Eigen::MatrixXd alpha = grid.getAlpha();
+  const auto &alpha = grid.getAlpha();
-  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-  std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
-  Eigen::MatrixXd concentrations = grid.getConcentrations();
+  const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
+  const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
+
+  Eigen::MatrixX<T> concentrations = grid.getConcentrations();
   int rowIndex = 0;
   A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex,
                         sx); // this is exactly same as in 2D
@@ -331,31 +340,32 @@ BTCS_1D(Grid &grid, Boundary &bc, double timestep,
 }
 
 // BTCS solution for 2D grid
-static void
+template <class T>
-BTCS_2D(Grid &grid, Boundary &bc, double timestep,
+static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
-        Eigen::VectorXd (*solverFunc)(Eigen::SparseMatrix<double> &A,
+                    Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
-                                      Eigen::VectorXd &b),
+                                                    Eigen::VectorX<T> &b),
-        int numThreads) {
+                    int numThreads) {
   int rowMax = grid.getRow();
   int colMax = grid.getCol();
-  double sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
+  T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
-  double sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
+  T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
 
-  Eigen::MatrixXd concentrations_t1 =
+  Eigen::MatrixX<T> concentrations_t1 =
-      Eigen::MatrixXd::Constant(rowMax, colMax, 0);
+      Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
-  Eigen::VectorXd row_t1(colMax);
+  Eigen::VectorX<T> row_t1(colMax);
 
-  Eigen::SparseMatrix<double> A;
+  Eigen::SparseMatrix<T> A;
-  Eigen::VectorXd b;
+  Eigen::VectorX<T> b;
 
-  Eigen::MatrixXd alphaX = grid.getAlphaX();
+  auto alphaX = grid.getAlphaX();
-  Eigen::MatrixXd alphaY = grid.getAlphaY();
+  auto alphaY = grid.getAlphaY();
-  std::vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-  std::vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
-  std::vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
-  std::vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
 
-  Eigen::MatrixXd concentrations = grid.getConcentrations();
+  const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
+  const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
+  const auto &bcTop = bc.getBoundarySide(BC_SIDE_TOP);
+  const auto &bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
+
+  Eigen::MatrixX<T> concentrations = grid.getConcentrations();
 
 #pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
   for (int i = 0; i < rowMax; i++) {
@@ -364,8 +374,6 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep,
     b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
                              bcTop, bcBottom, colMax, i, sx, sy);
 
-    Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
-
     row_t1 = solverFunc(A, b);
 
     concentrations_t1.row(i) = row_t1;
@@ -394,7 +402,8 @@ BTCS_2D(Grid &grid, Boundary &bc, double timestep,
 }
 
 // entry point for EigenLU solver; differentiate between 1D and 2D grid
-void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) {
+template <class T>
+void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
   if (grid.getDim() == 1) {
     BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
   } else if (grid.getDim() == 2) {
@@ -406,7 +415,8 @@ void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads) {
 }
 
 // entry point for Thomas algorithm solver; differentiate 1D and 2D grid
-void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) {
+template <class T>
+void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
   if (grid.getDim() == 1) {
     BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
   } else if (grid.getDim() == 2) {
@@ -416,3 +426,13 @@ void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep, int numThreads) {
         "Error: Only 1- and 2-dimensional grids are defined!");
   }
 }
+
+template void BTCS_Thomas(Grid<double> &grid, Boundary<double> &bc,
+                          double timestep, int numThreads);
+template void BTCS_Thomas(Grid<float> &grid, Boundary<float> &bc,
+                          float timestep, int numThreads);
+
+template void BTCS_LU(Grid<double> &grid, Boundary<double> &bc, double timestep,
+                      int numThreads);
+template void BTCS_LU(Grid<float> &grid, Boundary<float> &bc, float timestep,
+                      int numThreads);

src/Boundary.cpp

@@ -1,149 +0,0 @@
-#include "TugUtils.hpp"
-
-#include <cstddef>
-#include <stdexcept>
-#include <tug/Boundary.hpp>
-
-BoundaryElement::BoundaryElement() {}
-
-BoundaryElement::BoundaryElement(double _value)
-    : value(_value), type(BC_TYPE_CONSTANT) {}
-
-void BoundaryElement::setType(BC_TYPE type) { this->type = type; }
-
-void BoundaryElement::setValue(double value) {
-  if (value < 0) {
-    throw_invalid_argument("No negative concentration allowed.");
-  }
-  if (type == BC_TYPE_CLOSED) {
-    throw_invalid_argument(
-        "No constant boundary concentrations can be set for closed "
-        "boundaries. Please change type first.");
-  }
-  this->value = value;
-}
-
-BC_TYPE BoundaryElement::getType() { return this->type; }
-
-double BoundaryElement::getValue() { return this->value; }
-
-Boundary::Boundary(Grid grid) : grid(grid) {
-  if (grid.getDim() == 1) {
-    this->boundaries = std::vector<std::vector<BoundaryElement>>(
-        2); // in 1D only left and right boundary
-
-    this->boundaries[BC_SIDE_LEFT].push_back(BoundaryElement());
-    this->boundaries[BC_SIDE_RIGHT].push_back(BoundaryElement());
-  } else if (grid.getDim() == 2) {
-    this->boundaries = std::vector<std::vector<BoundaryElement>>(4);
-
-    this->boundaries[BC_SIDE_LEFT] =
-        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
-    this->boundaries[BC_SIDE_RIGHT] =
-        std::vector<BoundaryElement>(grid.getRow(), BoundaryElement());
-    this->boundaries[BC_SIDE_TOP] =
-        std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
-    this->boundaries[BC_SIDE_BOTTOM] =
-        std::vector<BoundaryElement>(grid.getCol(), BoundaryElement());
-  }
-}
-
-void Boundary::setBoundarySideClosed(BC_SIDE side) {
-  if (grid.getDim() == 1) {
-    if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
-      throw_invalid_argument(
-          "For the one-dimensional case, only the BC_SIDE_LEFT and "
-          "BC_SIDE_RIGHT borders exist.");
-    }
-  }
-
-  const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
-  const int n = is_vertical ? grid.getRow() : grid.getCol();
-
-  this->boundaries[side] = std::vector<BoundaryElement>(n, BoundaryElement());
-}
-
-void Boundary::setBoundarySideConstant(BC_SIDE side, double value) {
-  if (grid.getDim() == 1) {
-    if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
-      throw_invalid_argument(
-          "For the one-dimensional case, only the BC_SIDE_LEFT and "
-          "BC_SIDE_RIGHT borders exist.");
-    }
-  }
-
-  const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
-  const int n = is_vertical ? grid.getRow() : grid.getCol();
-
-  this->boundaries[side] =
-      std::vector<BoundaryElement>(n, BoundaryElement(value));
-}
-
-void Boundary::setBoundaryElementClosed(BC_SIDE side, int index) {
-  // tests whether the index really points to an element of the boundary side.
-  if ((boundaries[side].size() < index) || index < 0) {
-    throw_invalid_argument("Index is selected either too large or too small.");
-  }
-  this->boundaries[side][index].setType(BC_TYPE_CLOSED);
-}
-
-void Boundary::setBoundaryElementConstant(BC_SIDE side, int index,
-                                          double value) {
-  // tests whether the index really points to an element of the boundary side.
-  if ((boundaries[side].size() < index) || index < 0) {
-    throw_invalid_argument("Index is selected either too large or too small.");
-  }
-  this->boundaries[side][index].setType(BC_TYPE_CONSTANT);
-  this->boundaries[side][index].setValue(value);
-}
-
-const std::vector<BoundaryElement> Boundary::getBoundarySide(BC_SIDE side) {
-  if (grid.getDim() == 1) {
-    if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
-      throw_invalid_argument(
-          "For the one-dimensional trap, only the BC_SIDE_LEFT and "
-          "BC_SIDE_RIGHT borders exist.");
-    }
-  }
-  return this->boundaries[side];
-}
-
-Eigen::VectorXd Boundary::getBoundarySideValues(BC_SIDE side) {
-  const std::size_t length = boundaries[side].size();
-  Eigen::VectorXd values(length);
-
-  for (int i = 0; i < length; i++) {
-    if (getBoundaryElementType(side, i) == BC_TYPE_CLOSED) {
-      values(i) = -1;
-      continue;
-    }
-    values(i) = getBoundaryElementValue(side, i);
-  }
-
-  return values;
-}
-
-BoundaryElement Boundary::getBoundaryElement(BC_SIDE side, int index) {
-  if ((boundaries[side].size() < index) || index < 0) {
-    throw_invalid_argument("Index is selected either too large or too small.");
-  }
-  return this->boundaries[side][index];
-}
-
-BC_TYPE Boundary::getBoundaryElementType(BC_SIDE side, int index) {
-  if ((boundaries[side].size() < index) || index < 0) {
-    throw_invalid_argument("Index is selected either too large or too small.");
-  }
-  return this->boundaries[side][index].getType();
-}
-
-double Boundary::getBoundaryElementValue(BC_SIDE side, int index) {
-  if ((boundaries[side].size() < index) || index < 0) {
-    throw_invalid_argument("Index is selected either too large or too small.");
-  }
-  if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) {
-    throw_invalid_argument(
-        "A value can only be output if it is a constant boundary condition.");
-  }
-  return this->boundaries[side][index].getValue();
-}

src/CMakeLists.txt

@@ -1,4 +1,4 @@
-add_library(tug Boundary.cpp Grid.cpp Simulation.cpp FTCS.cpp BTCS.cpp)
+add_library(tug Simulation.cpp FTCS.cpp BTCS.cpp)
 
 IF(TUG_NAAICE_EXAMPLE)
   target_compile_definitions(tug PRIVATE WRITE_THOMAS_CSV)

src/FTCS.cpp

@@ -19,7 +19,8 @@
 #endif
 
 // calculates horizontal change on one cell independent of boundary type
-static inline double calcHorizontalChange(Grid &grid, int &row, int &col) {
+template <class T>
+static inline T calcHorizontalChange(Grid<T> &grid, int &row, int &col) {
 
   return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
                             grid.getAlphaX()(row, col)) *
@@ -35,7 +36,8 @@ static inline double calcHorizontalChange(Grid &grid, int &row, int &col) {
 }
 
 // calculates vertical change on one cell independent of boundary type
-static inline double calcVerticalChange(Grid &grid, int &row, int &col) {
+template <class T>
+static inline T calcVerticalChange(Grid<T> &grid, int &row, int &col) {
 
   return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
                             grid.getAlphaY()(row, col)) *
@@ -51,10 +53,10 @@ static inline double calcVerticalChange(Grid &grid, int &row, int &col) {
 }
 
 // calculates horizontal change on one cell with a constant left boundary
-static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid,
+template <class T>
-                                                              Boundary &bc,
+static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid,
-                                                              int &row,
+                                                         Boundary<T> &bc,
-                                                              int &col) {
+                                                         int &row, int &col) {
 
   return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
                             grid.getAlphaX()(row, col)) *
@@ -68,8 +70,9 @@ static inline double calcHorizontalChangeLeftBoundaryConstant(Grid &grid,
 }
 
 // calculates horizontal change on one cell with a closed left boundary
-static inline double
+template <class T>
-calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) {
+static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &grid, int &row,
+                                                       int &col) {
 
   return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
                             grid.getAlphaX()(row, col)) *
@@ -78,8 +81,9 @@ calcHorizontalChangeLeftBoundaryClosed(Grid &grid, int &row, int &col) {
 }
 
 // checks boundary condition type for a cell on the left edge of grid
-static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc,
+template <class T>
-                                                      int &row, int &col) {
+static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &bc,
+                                                 int &row, int &col) {
   if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CONSTANT) {
     return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
   } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, col) == BC_TYPE_CLOSED) {
@@ -90,10 +94,10 @@ static inline double calcHorizontalChangeLeftBoundary(Grid &grid, Boundary &bc,
 }
 
 // calculates horizontal change on one cell with a constant right boundary
-static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid,
+template <class T>
-                                                               Boundary &bc,
+static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid,
-                                                               int &row,
+                                                          Boundary<T> &bc,
-                                                               int &col) {
+                                                          int &row, int &col) {
 
   return 2 * grid.getAlphaX()(row, col) *
              bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) -
@@ -107,8 +111,9 @@ static inline double calcHorizontalChangeRightBoundaryConstant(Grid &grid,
 }
 
 // calculates horizontal change on one cell with a closed right boundary
-static inline double
+template <class T>
-calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) {
+static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &grid, int &row,
+                                                        int &col) {
 
   return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
                               grid.getAlphaX()(row, col)) *
@@ -117,8 +122,10 @@ calcHorizontalChangeRightBoundaryClosed(Grid &grid, int &row, int &col) {
 }
 
 // checks boundary condition type for a cell on the right edge of grid
-static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc,
+template <class T>
-                                                       int &row, int &col) {
+static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid,
+                                                  Boundary<T> &bc, int &row,
+                                                  int &col) {
   if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CONSTANT) {
     return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
   } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, col) == BC_TYPE_CLOSED) {
@@ -129,9 +136,10 @@ static inline double calcHorizontalChangeRightBoundary(Grid &grid, Boundary &bc,
 }
 
 // calculates vertical change on one cell with a constant top boundary
-static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid,
+template <class T>
-                                                           Boundary &bc,
+static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid,
-                                                           int &row, int &col) {
+                                                      Boundary<T> &bc, int &row,
+                                                      int &col) {
 
   return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
                             grid.getAlphaY()(row, col)) *
@@ -145,8 +153,9 @@ static inline double calcVerticalChangeTopBoundaryConstant(Grid &grid,
 }
 
 // calculates vertical change on one cell with a closed top boundary
-static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row,
+template <class T>
-                                                         int &col) {
+static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &grid, int &row,
+                                                    int &col) {
 
   return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
                             grid.getConcentrations()(row, col)) *
@@ -155,8 +164,9 @@ static inline double calcVerticalChangeTopBoundaryClosed(Grid &grid, int &row,
 }
 
 // checks boundary condition type for a cell on the top edge of grid
-static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc,
+template <class T>
-                                                   int &row, int &col) {
+static inline T calcVerticalChangeTopBoundary(Grid<T> &grid, Boundary<T> &bc,
+                                              int &row, int &col) {
   if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
     return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
   } else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
@@ -167,10 +177,10 @@ static inline double calcVerticalChangeTopBoundary(Grid &grid, Boundary &bc,
 }
 
 // calculates vertical change on one cell with a constant bottom boundary
-static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid,
+template <class T>
-                                                              Boundary &bc,
+static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid,
-                                                              int &row,
+                                                         Boundary<T> &bc,
-                                                              int &col) {
+                                                         int &row, int &col) {
 
   return 2 * grid.getAlphaY()(row, col) *
              bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) -
@@ -184,8 +194,9 @@ static inline double calcVerticalChangeBottomBoundaryConstant(Grid &grid,
 }
 
 // calculates vertical change on one cell with a closed bottom boundary
-static inline double
+template <class T>
-calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) {
+static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &grid, int &row,
+                                                       int &col) {
 
   return -(calcAlphaIntercell(grid.getAlphaY()(row, col),
                               grid.getAlphaY()(row - 1, col)) *
@@ -194,8 +205,9 @@ calcVerticalChangeBottomBoundaryClosed(Grid &grid, int &row, int &col) {
 }
 
 // checks boundary condition type for a cell on the bottom edge of grid
-static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc,
+template <class T>
-                                                      int &row, int &col) {
+static inline T calcVerticalChangeBottomBoundary(Grid<T> &grid, Boundary<T> &bc,
+                                                 int &row, int &col) {
   if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
     return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
   } else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
@@ -206,12 +218,14 @@ static inline double calcVerticalChangeBottomBoundary(Grid &grid, Boundary &bc,
 }
 
 // FTCS solution for 1D grid
-static void FTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
+template <class T>
+static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
   int colMax = grid.getCol();
-  double deltaCol = grid.getDeltaCol();
+  T deltaCol = grid.getDeltaCol();
 
   // matrix for concentrations at time t+1
-  Eigen::MatrixXd concentrations_t1 = Eigen::MatrixXd::Constant(1, colMax, 0);
+  Eigen::MatrixX<T> concentrations_t1 =
+      Eigen::MatrixX<T>::Constant(1, colMax, 0);
 
   // only one row in 1D case -> row constant at index 0
   int row = 0;
@@ -243,16 +257,17 @@ static void FTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
 }
 
 // FTCS solution for 2D grid
-static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
+template <class T>
+static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
                     int numThreads) {
   int rowMax = grid.getRow();
   int colMax = grid.getCol();
-  double deltaRow = grid.getDeltaRow();
+  T deltaRow = grid.getDeltaRow();
-  double deltaCol = grid.getDeltaCol();
+  T deltaCol = grid.getDeltaCol();
 
   // matrix for concentrations at time t+1
-  Eigen::MatrixXd concentrations_t1 =
+  Eigen::MatrixX<T> concentrations_t1 =
-      Eigen::MatrixXd::Constant(rowMax, colMax, 0);
+      Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
 
   // inner cells
   // these are independent of the boundary condition type
@@ -369,7 +384,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, double &timestep,
 }
 
 // entry point; differentiate between 1D and 2D grid
-void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
+template <class T>
+void FTCS(Grid<T> &grid, Boundary<T> &bc, T timestep, int &numThreads) {
   if (grid.getDim() == 1) {
     FTCS_1D(grid, bc, timestep);
   } else if (grid.getDim() == 2) {
@@ -379,3 +395,8 @@ void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads) {
         "Error: Only 1- and 2-dimensional grids are defined!");
   }
 }
+
+template void FTCS(Grid<double> &grid, Boundary<double> &bc, double timestep,
+                   int &numThreads);
+template void FTCS(Grid<float> &grid, Boundary<float> &bc, float timestep,
+                   int &numThreads);

src/Grid.cpp

@@ -1,167 +0,0 @@
-#include "TugUtils.hpp"
-
-#include <iostream>
-#include <tug/Grid.hpp>
-
-constexpr double MAT_INIT_VAL = 0;
-
-Grid::Grid(int length) : col(length), domainCol(length) {
-  if (length <= 3) {
-    throw_invalid_argument(
-        "Given grid length too small. Must be greater than 3.");
-  }
-
-  this->dim = 1;
-  this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
-
-  this->concentrations = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL);
-  this->alphaX = Eigen::MatrixXd::Constant(1, col, MAT_INIT_VAL);
-}
-
-Grid::Grid(int _row, int _col)
-    : row(_row), col(_col), domainRow(_row), domainCol(_col) {
-  if (row <= 3 || col <= 3) {
-    throw_invalid_argument(
-        "Given grid dimensions too small. Must each be greater than 3.");
-  }
-
-  this->dim = 2;
-  this->deltaRow = double(this->domainRow) / double(this->row); // -> 1
-  this->deltaCol = double(this->domainCol) / double(this->col); // -> 1
-
-  this->concentrations = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
-  this->alphaX = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
-  this->alphaY = Eigen::MatrixXd::Constant(row, col, MAT_INIT_VAL);
-}
-
-void Grid::setConcentrations(const Eigen::MatrixXd &concentrations) {
-  if (concentrations.rows() != this->row ||
-      concentrations.cols() != this->col) {
-    throw_invalid_argument(
-        "Given matrix of concentrations mismatch with Grid dimensions!");
-  }
-
-  this->concentrations = concentrations;
-}
-
-void Grid::setAlpha(const Eigen::MatrixXd &alpha) {
-  if (dim != 1) {
-    throw_invalid_argument("Grid is not one dimensional, you should probably "
-                           "use 2D setter function!");
-  }
-  if (alpha.rows() != 1 || alpha.cols() != this->col) {
-    throw_invalid_argument(
-        "Given matrix of alpha coefficients mismatch with Grid dimensions!");
-  }
-
-  this->alphaX = alpha;
-}
-
-void Grid::setAlpha(const Eigen::MatrixXd &alphaX,
-                    const Eigen::MatrixXd &alphaY) {
-  if (dim != 2) {
-    throw_invalid_argument("Grid is not two dimensional, you should probably "
-                           "use 1D setter function!");
-  }
-  if (alphaX.rows() != this->row || alphaX.cols() != this->col) {
-    throw_invalid_argument("Given matrix of alpha coefficients in x-direction "
-                           "mismatch with GRid dimensions!");
-  }
-  if (alphaY.rows() != this->row || alphaY.cols() != this->col) {
-    throw_invalid_argument("Given matrix of alpha coefficients in y-direction "
-                           "mismatch with GRid dimensions!");
-  }
-
-  this->alphaX = alphaX;
-  this->alphaY = alphaY;
-}
-
-const Eigen::MatrixXd &Grid::getAlpha() {
-  if (dim != 1) {
-    throw_invalid_argument("Grid is not one dimensional, you should probably "
-                           "use either getAlphaX() or getAlphaY()!");
-  }
-
-  return this->alphaX;
-}
-
-const Eigen::MatrixXd &Grid::getAlphaX() {
-  if (dim != 2) {
-    throw_invalid_argument(
-        "Grid is not two dimensional, you should probably use getAlpha()!");
-  }
-
-  return this->alphaX;
-}
-
-const Eigen::MatrixXd &Grid::getAlphaY() {
-  if (dim != 2) {
-    throw_invalid_argument(
-        "Grid is not two dimensional, you should probably use getAlpha()!");
-  }
-
-  return this->alphaY;
-}
-
-int Grid::getDim() { return dim; }
-
-int Grid::getLength() {
-  if (dim != 1) {
-    throw_invalid_argument("Grid is not one dimensional, you should probably "
-                           "use getRow() or getCol()!");
-  }
-
-  return col;
-}
-
-int Grid::getRow() { return row; }
-
-int Grid::getCol() { return col; }
-
-void Grid::setDomain(double domainLength) {
-  if (dim != 1) {
-    throw_invalid_argument("Grid is not one dimensional, you should probaly "
-                           "use the 2D domain setter!");
-  }
-  if (domainLength <= 0) {
-    throw_invalid_argument("Given domain length is not positive!");
-  }
-
-  this->domainCol = domainLength;
-  this->deltaCol = double(this->domainCol) / double(this->col);
-}
-
-void Grid::setDomain(double domainRow, double domainCol) {
-  if (dim != 2) {
-    throw_invalid_argument("Grid is not two dimensional, you should probably "
-                           "use the 1D domain setter!");
-  }
-  if (domainRow <= 0 || domainCol <= 0) {
-    throw_invalid_argument("Given domain size is not positive!");
-  }
-
-  this->domainRow = domainRow;
-  this->domainCol = domainCol;
-  this->deltaRow = double(this->domainRow) / double(this->row);
-  this->deltaCol = double(this->domainCol) / double(this->col);
-}
-
-double Grid::getDelta() {
-  if (dim != 1) {
-    throw_invalid_argument("Grid is not one dimensional, you should probably "
-                           "use the 2D delta getters");
-  }
-
-  return this->deltaCol;
-}
-
-double Grid::getDeltaCol() { return this->deltaCol; }
-
-double Grid::getDeltaRow() {
-  if (dim != 2) {
-    throw_invalid_argument("Grid is not two dimensional, meaning there is no "
-                           "delta in y-direction!");
-  }
-
-  return this->deltaRow;
-}

src/Schemes.hpp

@@ -16,13 +16,16 @@
 #include <tug/Grid.hpp>
 
 // entry point; differentiate between 1D and 2D grid
-extern void FTCS(Grid &grid, Boundary &bc, double &timestep, int &numThreads);
+template <class T>
+extern void FTCS(Grid<T> &grid, Boundary<T> &bc, T timestep, int &numThreads);
 
 // entry point for EigenLU solver; differentiate between 1D and 2D grid
-extern void BTCS_LU(Grid &grid, Boundary &bc, double timestep, int numThreads);
+template <class T>
+extern void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads);
 
 // entry point for Thomas algorithm solver; differentiate 1D and 2D grid
-extern void BTCS_Thomas(Grid &grid, Boundary &bc, double timestep,
+template <class T>
+extern void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep,
                         int numThreads);
 
 #endif // SCHEMES_H_

src/Simulation.cpp

@@ -12,56 +12,28 @@
 #include "Schemes.hpp"
 #include "TugUtils.hpp"
 
-Simulation::Simulation(Grid &_grid, Boundary &_bc, APPROACH _approach)
+template <class T> void Simulation<T>::setTimestep(T timestep) {
-    : grid(_grid), bc(_bc), approach(_approach) {}
-
-void Simulation::setOutputCSV(CSV_OUTPUT csv_output) {
-  if (csv_output < CSV_OUTPUT_OFF && csv_output > CSV_OUTPUT_VERBOSE) {
-    throw_invalid_argument("Invalid CSV output option given!");
-  }
-
-  this->csv_output = csv_output;
-}
-
-void Simulation::setOutputConsole(CONSOLE_OUTPUT console_output) {
-  if (console_output < CONSOLE_OUTPUT_OFF &&
-      console_output > CONSOLE_OUTPUT_VERBOSE) {
-    throw_invalid_argument("Invalid console output option given!");
-  }
-
-  this->console_output = console_output;
-}
-
-void Simulation::setTimeMeasure(TIME_MEASURE time_measure) {
-  if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) {
-    throw_invalid_argument("Invalid time measure option given!");
-  }
-
-  this->time_measure = time_measure;
-}
-
-void Simulation::setTimestep(double timestep) {
   if (timestep <= 0) {
     throw_invalid_argument("Timestep has to be greater than zero.");
   }
 
   if (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) {
 
-    const double deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
+    const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
     // will be 0 if 1D, else ...
-    const double deltaRowSquare = grid.getDim() != 1
+    const T deltaRowSquare = grid.getDim() != 1
-                                      ? grid.getDeltaRow() * grid.getDeltaRow()
+                                 ? grid.getDeltaRow() * grid.getDeltaRow()
-                                      : deltaColSquare;
+                                 : deltaColSquare;
-    const double minDeltaSquare =
+    const T minDeltaSquare =
         (deltaRowSquare < deltaColSquare) ? deltaRowSquare : deltaColSquare;
 
-    double maxAlpha = std::numeric_limits<double>::quiet_NaN();
+    T maxAlpha = std::numeric_limits<T>::quiet_NaN();
 
     // determine maximum alpha
     if (grid.getDim() == 2) {
 
-      const double maxAlphaX = grid.getAlphaX().maxCoeff();
+      const T maxAlphaX = grid.getAlphaX().maxCoeff();
-      const double maxAlphaY = grid.getAlphaY().maxCoeff();
+      const T maxAlphaY = grid.getAlphaY().maxCoeff();
       maxAlpha = (maxAlphaX > maxAlphaY) ? maxAlphaX : maxAlphaY;
 
     } else if (grid.getDim() == 1) {
@@ -74,7 +46,7 @@ void Simulation::setTimestep(double timestep) {
     const std::string dim = std::to_string(grid.getDim()) + "D";
 
     // Courant-Friedrichs-Lewy condition
-    double cfl = minDeltaSquare / (4 * maxAlpha);
+    T cfl = minDeltaSquare / (4 * maxAlpha);
 
     // stability equation from Wikipedia; might be useful if applied cfl does
     // not work in some cases double CFL_Wiki = 1 / (4 * maxAlpha *
@@ -112,109 +84,7 @@ void Simulation::setTimestep(double timestep) {
   }
 }
 
-double Simulation::getTimestep() { return this->timestep; }
+template <class T> void Simulation<T>::run() {
-
-void Simulation::setIterations(int iterations) {
-  if (iterations <= 0) {
-    throw_invalid_argument("Number of iterations must be greater than zero.");
-  }
-  this->iterations = iterations;
-}
-
-void Simulation::setSolver(SOLVER solver) {
-  if (this->approach == FTCS_APPROACH) {
-    std::cerr
-        << "Warning: Solver was set, but FTCS approach initialized. Setting "
-           "the solver "
-           "is thus without effect."
-        << std::endl;
-  }
-
-  this->solver = solver;
-}
-
-void Simulation::setNumberThreads(int numThreads) {
-  if (numThreads > 0 && numThreads <= omp_get_num_procs()) {
-    this->numThreads = numThreads;
-  } else {
-    int maxThreadNumber = omp_get_num_procs();
-    std::string outputMessage =
-        "Number of threads exceeds the number of processor cores (" +
-        std::to_string(maxThreadNumber) + ") or is less than 1.";
-
-    throw_invalid_argument(outputMessage);
-  }
-}
-
-int Simulation::getIterations() { return this->iterations; }
-
-void Simulation::printConcentrationsConsole() {
-  std::cout << grid.getConcentrations() << std::endl;
-  std::cout << std::endl;
-}
-
-std::string Simulation::createCSVfile() {
-  std::ofstream file;
-  int appendIdent = 0;
-  std::string appendIdentString;
-
-  // string approachString = (approach == 0) ? "FTCS" : "BTCS";
-  const std::string &approachString = this->approach_names[approach];
-  std::string row = std::to_string(grid.getRow());
-  std::string col = std::to_string(grid.getCol());
-  std::string numIterations = std::to_string(iterations);
-
-  std::string filename =
-      approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
-
-  while (std::filesystem::exists(filename)) {
-    appendIdent += 1;
-    appendIdentString = std::to_string(appendIdent);
-    filename = approachString + "_" + row + "_" + col + "_" + numIterations +
-               "-" + appendIdentString + ".csv";
-  }
-
-  file.open(filename);
-  if (!file) {
-    exit(1);
-  }
-
-  // adds lines at the beginning of verbose output csv that represent the
-  // boundary conditions and their values -1 in case of closed boundary
-  if (csv_output == CSV_OUTPUT_XTREME) {
-    Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "",
-                            " ");
-    file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row)
-         << std::endl; // boundary left
-    file << bc.getBoundarySideValues(BC_SIDE_RIGHT).format(one_row)
-         << std::endl; // boundary right
-    file << bc.getBoundarySideValues(BC_SIDE_TOP).format(one_row)
-         << std::endl; // boundary top
-    file << bc.getBoundarySideValues(BC_SIDE_BOTTOM).format(one_row)
-         << std::endl; // boundary bottom
-    file << std::endl << std::endl;
-  }
-
-  file.close();
-
-  return filename;
-}
-
-void Simulation::printConcentrationsCSV(const std::string &filename) {
-  std::ofstream file;
-
-  file.open(filename, std::ios_base::app);
-  if (!file) {
-    exit(1);
-  }
-
-  Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols);
-  file << grid.getConcentrations().format(do_not_align) << std::endl;
-  file << std::endl << std::endl;
-  file.close();
-}
-
-void Simulation::run() {
   if (this->timestep == -1) {
     throw_invalid_argument("Timestep is not set!");
   }
@@ -280,14 +150,14 @@ void Simulation::run() {
 
   } else if (approach == CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case
 
-    double beta = 0.5;
+    constexpr T beta = 0.5;
 
     // TODO this implementation is very inefficient!
     // a separate implementation that sets up a specific tridiagonal matrix for
     // Crank-Nicolson would be better
-    Eigen::MatrixXd concentrations;
+    Eigen::MatrixX<T> concentrations;
-    Eigen::MatrixXd concentrationsFTCS;
+    Eigen::MatrixX<T> concentrationsFTCS;
-    Eigen::MatrixXd concentrationsResult;
+    Eigen::MatrixX<T> concentrationsResult;
     for (int i = 0; i < iterations * innerIterations; i++) {
       if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
         printConcentrationsConsole();
@@ -324,3 +194,9 @@ void Simulation::run() {
               << milliseconds.count() << "ms" << std::endl;
   }
 }
+
+template void Simulation<double>::setTimestep(double timestep);
+template void Simulation<float>::setTimestep(float timestep);
+
+template void Simulation<double>::run();
+template void Simulation<float>::run();

test/testBoundary.cpp

@@ -6,12 +6,11 @@
 #include <typeinfo>
 
 using namespace std;
-
 TEST_CASE("BoundaryElement") {
 
   SUBCASE("Closed case") {
-    BoundaryElement boundaryElementClosed = BoundaryElement();
+    BoundaryElement boundaryElementClosed = BoundaryElement<double>();
-    CHECK_NOTHROW(BoundaryElement());
+    CHECK_NOTHROW(BoundaryElement<double>());
     CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED);
     CHECK_EQ(boundaryElementClosed.getValue(), -1);
     CHECK_THROWS(boundaryElementClosed.setValue(0.2));
@@ -28,11 +27,11 @@ TEST_CASE("BoundaryElement") {
 }
 
 TEST_CASE("Boundary Class") {
-  Grid grid1D = Grid(10);
+  Grid grid1D = Grid64(10);
-  Grid grid2D = Grid(10, 12);
+  Grid grid2D = Grid64(10, 12);
   Boundary boundary1D = Boundary(grid1D);
   Boundary boundary2D = Boundary(grid2D);
-  vector<BoundaryElement> boundary1DVector(1, BoundaryElement(1.0));
+  vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));
 
   SUBCASE("Boundaries 1D case") {
     CHECK_NOTHROW(Boundary boundary(grid1D));

test/testGrid.cpp

@@ -8,22 +8,22 @@ using namespace std;
 
 TEST_CASE("1D Grid, too small length") {
   int l = 2;
-  CHECK_THROWS(Grid(l));
+  CHECK_THROWS(Grid64(l));
 }
 
-TEST_CASE("2D Grid, too small side") {
+TEST_CASE("2D Grid64, too small side") {
   int r = 2;
   int c = 4;
-  CHECK_THROWS(Grid(r, c));
+  CHECK_THROWS(Grid64(r, c));
 
   r = 4;
   c = 2;
-  CHECK_THROWS(Grid(r, c));
+  CHECK_THROWS(Grid64(r, c));
 }
 
-TEST_CASE("1D Grid") {
+TEST_CASE("1D Grid64") {
   int l = 12;
-  Grid grid(l);
+  Grid64 grid(l);
 
   SUBCASE("correct construction") {
     CHECK_EQ(grid.getDim(), 1);
@@ -89,9 +89,9 @@ TEST_CASE("1D Grid") {
   }
 }
 
-TEST_CASE("2D Grid quadratic") {
+TEST_CASE("2D Grid64 quadratic") {
   int rc = 12;
-  Grid grid(rc, rc);
+  Grid64 grid(rc, rc);
 
   SUBCASE("correct construction") {
     CHECK_EQ(grid.getDim(), 2);
@@ -170,10 +170,10 @@ TEST_CASE("2D Grid quadratic") {
   }
 }
 
-TEST_CASE("2D Grid non-quadratic") {
+TEST_CASE("2D Grid64 non-quadratic") {
   int r = 12;
   int c = 15;
-  Grid grid(r, c);
+  Grid64 grid(r, c);
 
   SUBCASE("correct construction") {
     CHECK_EQ(grid.getDim(), 2);

test/testSimulation.cpp

@@ -11,7 +11,7 @@
 using namespace Eigen;
 using namespace std;
 
-static Grid setupSimulation(APPROACH approach, double timestep,
+static Grid64 setupSimulation(APPROACH approach, double timestep,
                             int iterations) {
   int row = 11;
   int col = 11;
@@ -19,7 +19,7 @@ static Grid setupSimulation(APPROACH approach, double timestep,
   int domain_col = 10;
 
   // Grid
-  Grid grid = Grid(row, col);
+  Grid grid = Grid64(row, col);
   grid.setDomain(domain_row, domain_col);
 
   MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
@@ -80,7 +80,7 @@ TEST_CASE("equality to reference matrix with BTCS") {
 
 TEST_CASE("Initialize environment") {
   int rc = 12;
-  Grid grid(rc, rc);
+  Grid64 grid(rc, rc);
   Boundary boundary(grid);
 
   CHECK_NOTHROW(Simulation sim(grid, boundary, BTCS_APPROACH));
@@ -88,7 +88,7 @@ TEST_CASE("Initialize environment") {
 
 TEST_CASE("Simulation environment") {
   int rc = 12;
-  Grid grid(rc, rc);
+  Grid64 grid(rc, rc);
   Boundary boundary(grid);
   Simulation sim(grid, boundary, FTCS_APPROACH);