Merge pull request 'BREAKING CHANGE: Refactor simulation grid' (#1) from ml/refactor-simulation-grid into main

.gitlab-ci.yml

@@ -12,7 +12,7 @@ test:
     stage: test
     script:
       - mkdir build && cd build
-      - cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON -G Ninja ..
+      - cmake -DCMAKE_BUILD_TYPE=Debug -DTUG_ENABLE_TESTING=ON -G Ninja ..
       - ninja 
       - ctest --output-junit test_results.xml
     artifacts:
@@ -27,7 +27,7 @@ pages:
   image: python:slim
   before_script:
     - apt-get update && apt-get install --no-install-recommends -y graphviz imagemagick doxygen make
-    - pip install --upgrade pip && pip install Sphinx Pillow breathe sphinx-rtd-theme m2r2
+    - pip install --upgrade pip && pip install Sphinx Pillow breathe sphinx-rtd-theme sphinx-mdinclude
     - mkdir public
   script:
     - pushd docs_sphinx

docs_sphinx/Diffusion.rst

@@ -1,4 +1,4 @@
-Simulation
+Diffusion
 ==========
 
 .. doxygenenum:: tug::APPROACH
@@ -7,4 +7,4 @@ Simulation
 .. doxygenenum:: tug::CONSOLE_OUTPUT
 .. doxygenenum:: tug::TIME_MEASURE
 
-.. doxygenclass:: tug::Simulation
+.. doxygenclass:: tug::Diffusion

docs_sphinx/Grid.rst

@@ -1,4 +0,0 @@
-Grid
-====
-
-.. doxygenclass:: tug::Grid

docs_sphinx/conf.py

@@ -42,7 +42,7 @@ extensions = [
     'sphinx.ext.viewcode',
     'sphinx.ext.inheritance_diagram',
     'breathe',
-    'm2r2'
+    'sphinx_mdinclude'
 ]
 
 html_baseurl = "/index.html"

docs_sphinx/user.rst

@@ -3,6 +3,5 @@ User API
 
 .. toctree::    
   
-    Grid
     Boundary
-    Simulation
+    Diffusion

examples/BTCS_1D_proto_example.cpp

@@ -1,51 +0,0 @@
-#include <Eigen/Eigen>
-#include <tug/Simulation.hpp>
-
-using namespace Eigen;
-using namespace tug;
-
-int main(int argc, char *argv[]) {
-  // **************
-  // **** GRID ****
-  // **************
-
-  // create a linear grid with 20 cells
-  int cells = 20;
-  Grid64 grid(cells);
-
-  MatrixXd concentrations = MatrixXd::Constant(1, 20, 0);
-  concentrations(0, 0) = 2000;
-  // TODO add option to set concentrations with a vector in 1D case
-  grid.setConcentrations(concentrations);
-
-  // ******************
-  // **** BOUNDARY ****
-  // ******************
-
-  // create a boundary with constant values
-  Boundary bc = Boundary(grid);
-  bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
-  bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
-
-  // ************************
-  // **** SIMULATION ENV ****
-  // ************************
-
-  // set up a simulation environment
-  Simulation simulation = Simulation(grid, bc); // grid,boundary
-
-  // set the timestep of the simulation
-  simulation.setTimestep(0.1); // timestep
-
-  // set the number of iterations
-  simulation.setIterations(100);
-
-  // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
-  // CSV_OUTPUT_VERBOSE]
-  simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
-
-  // **** RUN SIMULATION ****
-
-  // run the simulation
-  simulation.run();
-}

examples/BTCS_2D_proto_example.cpp

@@ -1,5 +1,5 @@
 #include <Eigen/Eigen>
-#include <tug/Simulation.hpp>
+#include <tug/Diffusion.hpp>
 
 using namespace Eigen;
 using namespace tug;
@@ -14,7 +14,6 @@ int main(int argc, char *argv[]) {
   // create a grid with a 20 x 20 field
   int row = 40;
   int col = 50;
-  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   // grid.setDomain(20, 20);
@@ -24,7 +23,7 @@ int main(int argc, char *argv[]) {
   // #row,#col,value grid.setConcentrations(concentrations);
   MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
   concentrations(10, 10) = 2000;
-  grid.setConcentrations(concentrations);
+  Grid64 grid(concentrations);
 
   // (optional) set alphas of the grid, e.g.:
   // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value
@@ -61,8 +60,7 @@ int main(int argc, char *argv[]) {
   // ************************
 
   // set up a simulation environment
-  Simulation simulation =
+  Diffusion simulation(grid, bc); // grid,boundary
-      Simulation(grid, bc); // grid,boundary
 
   // set the timestep of the simulation
   simulation.setTimestep(0.1); // timestep

examples/CMakeLists.txt

@@ -1,20 +1,7 @@
-add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
-add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
-add_executable(BTCS_1D_proto_example BTCS_1D_proto_example.cpp)
 add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
-add_executable(CRNI_2D_proto_example CRNI_2D_proto_example.cpp)
-add_executable(reference-FTCS_2D_closed reference-FTCS_2D_closed.cpp)
-add_executable(profiling_openmp profiling_openmp.cpp)
-
-target_link_libraries(FTCS_1D_proto_example tug)
-target_link_libraries(FTCS_2D_proto_example tug)
-target_link_libraries(BTCS_1D_proto_example tug)
-target_link_libraries(BTCS_2D_proto_example tug)
-target_link_libraries(CRNI_2D_proto_example tug)
-target_link_libraries(reference-FTCS_2D_closed tug)
-target_link_libraries(profiling_openmp tug)
-
 add_executable(FTCS_2D_proto_example_mdl FTCS_2D_proto_example_mdl.cpp)
 add_executable(FTCS_2D_proto_closed_mdl FTCS_2D_proto_closed_mdl.cpp)
+
+target_link_libraries(BTCS_2D_proto_example tug)
 target_link_libraries(FTCS_2D_proto_closed_mdl tug)
 target_link_libraries(FTCS_2D_proto_example_mdl tug)

examples/CRNI_2D_proto_example.cpp

@@ -1,29 +0,0 @@
-#include <Eigen/Eigen>
-#include <tug/Simulation.hpp>
-
-using namespace Eigen;
-using namespace tug;
-
-int main(int argc, char *argv[]) {
-  int row = 20;
-  int col = 20;
-  Grid64 grid(row, col);
-
-  MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
-  concentrations(10, 10) = 2000;
-  grid.setConcentrations(concentrations);
-
-  Boundary bc = Boundary(grid);
-  bc.setBoundarySideClosed(BC_SIDE_LEFT);
-  bc.setBoundarySideClosed(BC_SIDE_RIGHT);
-  bc.setBoundarySideClosed(BC_SIDE_TOP);
-  bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
-
-  Simulation simulation =
-      Simulation<double, tug::CRANK_NICOLSON_APPROACH>(grid, bc);
-  simulation.setTimestep(0.1);
-  simulation.setIterations(50);
-  simulation.setOutputCSV(CSV_OUTPUT_XTREME);
-
-  simulation.run();
-}

examples/FTCS_1D_proto_example.cpp

@@ -1,51 +0,0 @@
-#include "tug/Boundary.hpp"
-#include <Eigen/Eigen>
-#include <tug/Simulation.hpp>
-
-using namespace Eigen;
-using namespace tug;
-
-int main(int argc, char *argv[]) {
-  // **************
-  // **** GRID ****
-  // **************
-
-  // create a linear grid with 20 cells
-  int cells = 20;
-  Grid64 grid(cells);
-
-  MatrixXd concentrations = MatrixXd::Constant(1, 20, 20);
-  grid.setConcentrations(concentrations);
-
-  // ******************
-  // **** BOUNDARY ****
-  // ******************
-
-  // create a boundary with constant values
-  Boundary bc = Boundary(grid);
-  bc.setBoundarySideConstant(BC_SIDE_LEFT, 1);
-  bc.setBoundarySideConstant(BC_SIDE_RIGHT, 1);
-
-  // ************************
-  // **** SIMULATION ENV ****
-  // ************************
-
-  // set up a simulation environment
-  Simulation simulation =
-      Simulation<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
-
-  // (optional) set the timestep of the simulation
-  simulation.setTimestep(0.1); // timestep
-
-  // (optional) set the number of iterations
-  simulation.setIterations(100);
-
-  // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
-  // CSV_OUTPUT_VERBOSE]
-  simulation.setOutputCSV(CSV_OUTPUT_OFF);
-
-  // **** RUN SIMULATION ****
-
-  // run the simulation
-  simulation.run();
-}

examples/FTCS_2D_proto_closed_mdl.cpp

@@ -9,7 +9,7 @@
 #include <Eigen/Eigen>
 #include <cstdlib>
 #include <iostream>
-#include <tug/Simulation.hpp>
+#include <tug/Diffusion.hpp>
 
 using namespace Eigen;
 using namespace tug;
@@ -31,7 +31,6 @@ int main(int argc, char *argv[]) {
 
   // create a grid with a 20 x 20 field
   int n2 = row / 2 - 1;
-  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   // grid.setDomain(20, 20);
@@ -44,7 +43,7 @@ int main(int argc, char *argv[]) {
   concentrations(n2, n2 + 1) = 1;
   concentrations(n2 + 1, n2) = 1;
   concentrations(n2 + 1, n2 + 1) = 1;
-  grid.setConcentrations(concentrations);
+  Grid64 grid(concentrations);
 
   // (optional) set alphas of the grid, e.g.:
   MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value
@@ -69,8 +68,8 @@ int main(int argc, char *argv[]) {
   // ************************
 
   // set up a simulation environment
-  Simulation simulation =
+  Diffusion<double, FTCS_APPROACH> simulation(
-      Simulation<double, FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
+      grid, bc); // grid,boundary,simulation-approach
 
   // set the timestep of the simulation
   simulation.setTimestep(10000); // timestep

examples/FTCS_2D_proto_example.cpp

@@ -1,92 +0,0 @@
-/**
- * @file FTCS_2D_proto_example.cpp
- * @author Hannes Signer, Philipp Ungrund
- * @brief Creates a prototypical standard TUG simulation in 2D with FTCS
- * approach and constant boundary condition
- *
- */
-
-#include <Eigen/Eigen>
-#include <tug/Simulation.hpp>
-
-using namespace Eigen;
-using namespace tug;
-// #include <easy/profiler.h>
-// #define EASY_PROFILER_ENABLE ::profiler::setEnabled(true);
-
-int main(int argc, char *argv[]) {
-  // EASY_PROFILER_ENABLE;
-  // profiler::startListen();
-  // **************
-  // **** GRID ****
-  // **************
-  // profiler::startListen();
-  // create a grid with a 20 x 20 field
-  int row = 20;
-  int col = 20;
-  Grid64 grid(row, col);
-
-  // (optional) set the domain, e.g.:
-  // grid.setDomain(20, 20);
-
-  // (optional) set the concentrations, e.g.:
-  // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); //
-  // #row,#col,value grid.setConcentrations(concentrations);
-  MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
-  concentrations(0, 0) = 1999;
-  grid.setConcentrations(concentrations);
-
-  // (optional) set alphas of the grid, e.g.:
-  // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value
-  // MatrixXd alphay = MatrixXd::Constant(20,20,1); // row,col,value
-  // grid.setAlpha(alphax, alphay);
-
-  // ******************
-  // **** BOUNDARY ****
-  // ******************
-
-  // create a boundary with constant values
-  Boundary bc = Boundary(grid);
-  bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
-  bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
-  bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
-  bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
-
-  // (optional) set boundary condition values for one side, e.g.:
-  // VectorXd bc_left_values = VectorXd::Constant(20,1); // length,value
-  // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_left_values); // side,values
-  // VectorXd bc_zero_values = VectorXd::Constant(20,0);
-  // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
-  // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
-  // VectorXd bc_front_values = VectorXd::Constant(20,2000);
-  // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
-  // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
-
-  // ************************
-  // **** SIMULATION ENV ****
-  // ************************
-
-  // set up a simulation environment
-  Simulation simulation =
-      Simulation<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
-
-  // set the timestep of the simulation
-  simulation.setTimestep(0.1); // timestep
-
-  // set the number of iterations
-  simulation.setIterations(10000);
-
-  // set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON,
-  // CSV_OUTPUT_VERBOSE]
-  simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
-
-  // **** RUN SIMULATION ****
-
-  // run the simulation
-
-  // EASY_BLOCK("SIMULATION")
-  simulation.run();
-  // EASY_END_BLOCK;
-  // profiler::dumpBlocksToFile("test_profile.prof");
-  // profiler::stopListen();
-}

examples/FTCS_2D_proto_example_mdl.cpp

@@ -7,7 +7,7 @@
  */
 
 #include <Eigen/Eigen>
-#include <tug/Simulation.hpp>
+#include <tug/Diffusion.hpp>
 
 using namespace Eigen;
 using namespace tug;
@@ -22,7 +22,6 @@ int main(int argc, char *argv[]) {
   int row = 64;
   int col = 64;
   int n2 = row / 2 - 1;
-  Grid64 grid(row, col);
 
   // (optional) set the domain, e.g.:
   // grid.setDomain(20, 20);
@@ -35,7 +34,7 @@ int main(int argc, char *argv[]) {
   concentrations(n2, n2 + 1) = 1;
   concentrations(n2 + 1, n2) = 1;
   concentrations(n2 + 1, n2 + 1) = 1;
-  grid.setConcentrations(concentrations);
+  Grid64 grid(concentrations);
 
   // (optional) set alphas of the grid, e.g.:
   MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value
@@ -60,8 +59,8 @@ int main(int argc, char *argv[]) {
   // ************************
 
   // set up a simulation environment
-  Simulation simulation =
+  Diffusion<double, tug::FTCS_APPROACH> simulation(
-      Simulation<double, tug::FTCS_APPROACH>(grid, bc); // grid,boundary,simulation-approach
+      grid, bc); // grid,boundary,simulation-approach
 
   // (optional) set the timestep of the simulation
   simulation.setTimestep(1000); // timestep

examples/profiling_openmp.cpp

@@ -1,70 +0,0 @@
-#include <Eigen/Eigen>
-#include <chrono>
-#include <fstream>
-#include <iostream>
-#include <string>
-#include <tug/Simulation.hpp>
-
-using namespace Eigen;
-using namespace std;
-using namespace tug;
-
-int main(int argc, char *argv[]) {
-
-  int n[] = {2000};
-  int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
-  int iterations[1] = {1};
-  int repetition = 10;
-
-  for (int l = 0; l < size(threads); l++) {
-    // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv";
-    ofstream myfile;
-    myfile.open("speedup_1000.csv", std::ios::app);
-    myfile << "Number threads: " << threads[l] << endl;
-
-    for (int i = 0; i < size(n); i++) {
-      cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
-      // myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
-      for (int j = 0; j < size(iterations); j++) {
-        cout << "Iterations: " << iterations[j] << endl;
-        // myfile << "Iterations: " << iterations[j] << endl;
-        for (int k = 0; k < repetition; k++) {
-          cout << "Wiederholung: " << k << endl;
-          Grid64 grid(n[i], n[i]);
-          grid.setDomain(1, 1);
-
-          MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);
-          concentrations(n[i] / 2, n[i] / 2) = 1;
-          grid.setConcentrations(concentrations);
-          MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5);
-
-          Boundary bc = Boundary(grid);
-
-          Simulation sim = Simulation(grid, bc);
-
-          if (argc == 2) {
-            int numThreads = atoi(argv[1]);
-            sim.setNumberThreads(numThreads);
-          } else {
-            sim.setNumberThreads(threads[l]);
-          }
-
-          sim.setTimestep(0.01);
-          sim.setIterations(iterations[j]);
-          sim.setOutputCSV(CSV_OUTPUT_OFF);
-
-          auto begin = std::chrono::high_resolution_clock::now();
-          sim.run();
-          auto end = std::chrono::high_resolution_clock::now();
-          auto milliseconds =
-              std::chrono::duration_cast<std::chrono::milliseconds>(end -
-                                                                    begin);
-          myfile << milliseconds.count() << endl;
-        }
-      }
-      cout << endl;
-      myfile << endl;
-    }
-    myfile.close();
-  }
-}

examples/profiling_speedup.cpp

@@ -1,70 +0,0 @@
-#include <Eigen/Eigen>
-#include <chrono>
-#include <fstream>
-#include <iostream>
-#include <string>
-#include <tug/Simulation.hpp>
-
-using namespace Eigen;
-using namespace std;
-using namespace tug;
-
-int main(int argc, char *argv[]) {
-
-  int n[] = {2000};
-  int threads[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
-  int iterations[1] = {1};
-  int repetition = 10;
-
-  for (int l = 0; l < size(threads); l++) {
-    // string filename = "ftcs_openmp_" + to_string(threads[l]) + ".csv";
-    ofstream myfile;
-    myfile.open("speedup_1000.csv", std::ios::app);
-    myfile << "Number threads: " << threads[l] << endl;
-
-    for (int i = 0; i < size(n); i++) {
-      cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
-      // myfile << "Grid size: " << n[i] << " x " << n[i] << endl << endl;
-      for (int j = 0; j < size(iterations); j++) {
-        cout << "Iterations: " << iterations[j] << endl;
-        // myfile << "Iterations: " << iterations[j] << endl;
-        for (int k = 0; k < repetition; k++) {
-          cout << "Wiederholung: " << k << endl;
-          Grid64 grid(n[i], n[i]);
-          grid.setDomain(1, 1);
-
-          MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0);
-          concentrations(n[i] / 2, n[i] / 2) = 1;
-          grid.setConcentrations(concentrations);
-          MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5);
-
-          Boundary bc = Boundary(grid);
-
-          Simulation sim = Simulation(grid, bc);
-
-          if (argc == 2) {
-            int numThreads = atoi(argv[1]);
-            sim.setNumberThreads(numThreads);
-          } else {
-            sim.setNumberThreads(threads[l]);
-          }
-
-          sim.setTimestep(0.01);
-          sim.setIterations(iterations[j]);
-          sim.setOutputCSV(CSV_OUTPUT_OFF);
-
-          auto begin = std::chrono::high_resolution_clock::now();
-          sim.run();
-          auto end = std::chrono::high_resolution_clock::now();
-          auto milliseconds =
-              std::chrono::duration_cast<std::chrono::milliseconds>(end -
-                                                                    begin);
-          myfile << milliseconds.count() << endl;
-        }
-      }
-      cout << endl;
-      myfile << endl;
-    }
-    myfile.close();
-  }
-}

examples/reference-FTCS_2D_closed.cpp

@@ -1,53 +0,0 @@
-#include "Eigen/Core"
-#include <iostream>
-#include <tug/Simulation.hpp>
-
-using namespace std;
-using namespace Eigen;
-using namespace tug;
-
-int main(int argc, char *argv[]) {
-  int row = 50;
-  int col = 50;
-  int domain_row = 10;
-  int domain_col = 10;
-
-  // Grid
-  Grid64 grid(row, col);
-  grid.setDomain(domain_row, domain_col);
-
-  MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
-  concentrations(5, 5) = 1;
-  grid.setConcentrations(concentrations);
-
-  MatrixXd alpha = MatrixXd::Constant(row, col, 1);
-  for (int i = 0; i < 5; i++) {
-    for (int j = 0; j < 6; j++) {
-      alpha(i, j) = 0.01;
-    }
-  }
-  for (int i = 0; i < 5; i++) {
-    for (int j = 6; j < 11; j++) {
-      alpha(i, j) = 0.001;
-    }
-  }
-  for (int i = 5; i < 11; i++) {
-    for (int j = 6; j < 11; j++) {
-      alpha(i, j) = 0.1;
-    }
-  }
-  grid.setAlpha(alpha, alpha);
-
-  // Boundary
-  Boundary bc = Boundary(grid);
-
-  // Simulation
-  Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc);
-  sim.setTimestep(0.001);
-  sim.setIterations(10000);
-  sim.setOutputCSV(CSV_OUTPUT_OFF);
-  sim.setOutputConsole(CONSOLE_OUTPUT_OFF);
-
-  // RUN
-  sim.run();
-}

include/tug/Boundary.hpp

@@ -7,11 +7,13 @@
 #ifndef BOUNDARY_H_
 #define BOUNDARY_H_
 
-#include "Grid.hpp"
+#include "tug/Core/TugUtils.hpp"
 
+#include <Eigen/Dense>
 #include <cstddef>
 #include <cstdint>
 #include <map>
+#include <stdexcept>
 #include <utility>
 #include <vector>
 
@@ -114,7 +116,7 @@ public:
    *
    * @param length Length of the grid
    */
-  Boundary(std::uint32_t length) : Boundary(Grid<T>(length)){};
+  Boundary(std::uint32_t length) : Boundary(1, length) {};
 
   /**
    * @brief Creates a boundary object for a 2D grid
@@ -123,17 +125,7 @@ public:
    * @param cols Number of columns of the grid
    */
   Boundary(std::uint32_t rows, std::uint32_t cols)
-      : Boundary(Grid<T>(rows, cols)){};
+      : dim(rows == 1 ? 1 : 2), cols(cols), rows(rows) {
-
-  /**
-   * @brief Creates a boundary object based on the passed grid object and
-   *        initializes the boundaries as closed.
-   *
-   * @param grid Grid object on the basis of which the simulation takes place
-   *             and from which the dimensions (in 2D case) are taken.
-   */
-  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
@@ -152,8 +144,37 @@ public:
       this->boundaries[BC_SIDE_BOTTOM] =
           std::vector<BoundaryElement<T>>(this->cols, BoundaryElement<T>());
     }
-  }
+  };
 
+  /**
+   * @brief Creates a boundary object based on the passed grid object and
+   *        initializes the boundaries as closed.
+   *
+   * @param grid Grid object on the basis of which the simulation takes place
+   *             and from which the dimensions (in 2D case) are taken.
+   */
+  // 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
    *        condition closed.
@@ -161,13 +182,11 @@ public:
    * @param side Side to be set to closed, e.g. BC_SIDE_LEFT.
    */
   void setBoundarySideClosed(BC_SIDE side) {
-    if (this->dim == 1) {
+    tug_assert((this->dim > 1) ||
-      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
+                   ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
-        throw std::invalid_argument(
+               "For the "
-            "For the one-dimensional case, only the BC_SIDE_LEFT and "
+               "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
-            "BC_SIDE_RIGHT borders exist.");
+               "borders exist.");
-      }
-    }
 
     const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
     const int n = is_vertical ? this->rows : this->cols;
@@ -186,13 +205,11 @@ public:
    * page.
    */
   void setBoundarySideConstant(BC_SIDE side, double value) {
-    if (this->dim == 1) {
+    tug_assert((this->dim > 1) ||
-      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
+                   ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
-        throw std::invalid_argument(
+               "For the "
-            "For the one-dimensional case, only the BC_SIDE_LEFT and "
+               "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
-            "BC_SIDE_RIGHT borders exist.");
+               "borders exist.");
-      }
-    }
 
     const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
     const int n = is_vertical ? this->rows : this->cols;
@@ -212,10 +229,9 @@ public:
    */
   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) {
+    tug_assert(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);
   }
 
@@ -233,10 +249,8 @@ public:
    */
   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) {
+    tug_assert(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);
   }
@@ -251,13 +265,11 @@ public:
    * BoundaryElement<T> objects.
    */
   const std::vector<BoundaryElement<T>> &getBoundarySide(BC_SIDE side) const {
-    if (this->dim == 1) {
+    tug_assert((this->dim > 1) ||
-      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
+                   ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
-        throw std::invalid_argument(
+               "For the "
-            "For the one-dimensional trap, only the BC_SIDE_LEFT and "
+               "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
-            "BC_SIDE_RIGHT borders exist.");
+               "borders exist.");
-      }
-    }
     return this->boundaries[side];
   }
 
@@ -295,10 +307,8 @@ public:
    * object.
    */
   BoundaryElement<T> getBoundaryElement(BC_SIDE side, int index) const {
-    if ((boundaries[side].size() < index) || index < 0) {
+    tug_assert(boundaries[side].size() > index && index >= 0,
-      throw std::invalid_argument(
                "Index is selected either too large or too small.");
-    }
     return this->boundaries[side][index];
   }
 
@@ -313,10 +323,8 @@ public:
    * @return Boundary Type of the corresponding boundary condition.
    */
   BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const {
-    if ((boundaries[side].size() < index) || index < 0) {
+    tug_assert(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();
   }
 
@@ -333,14 +341,12 @@ public:
    * object.
    */
   T getBoundaryElementValue(BC_SIDE side, int index) const {
-    if ((boundaries[side].size() < index) || index < 0) {
+    tug_assert(boundaries[side].size() > index && index >= 0,
-      throw std::invalid_argument(
                "Index is selected either too large or too small.");
-    }
+    tug_assert(
-    if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) {
+        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();
   }
 
@@ -382,13 +388,8 @@ public:
    * @param value Value of the inner constant boundary condition
    */
   void setInnerBoundary(std::uint32_t index, T value) {
-    if (this->dim != 1) {
+    tug_assert(this->dim == 1, "This function is only available for 1D grids.");
-      throw std::invalid_argument(
+    tug_assert(index < this->cols, "Index is out of bounds.");
-          "This function is only available for 1D grids.");
-    }
-    if (index >= this->cols) {
-      throw std::invalid_argument("Index is out of bounds.");
-    }
 
     this->inner_boundary[std::make_pair(0, index)] = value;
   }
@@ -401,13 +402,8 @@ public:
    * @param value Value of the inner constant boundary condition
    */
   void setInnerBoundary(std::uint32_t row, std::uint32_t col, T value) {
-    if (this->dim != 2) {
+    tug_assert(this->dim == 2, "This function is only available for 2D grids.");
-      throw std::invalid_argument(
+    tug_assert(row < this->rows && col < this->cols, "Index is out of bounds.");
-          "This function is only available for 2D grids.");
-    }
-    if (row >= this->rows || col >= this->cols) {
-      throw std::invalid_argument("Index is out of bounds.");
-    }
 
     this->inner_boundary[std::make_pair(row, col)] = value;
   }
@@ -420,13 +416,8 @@ public:
    * set or not) and value of the inner constant boundary condition
    */
   std::pair<bool, T> getInnerBoundary(std::uint32_t index) const {
-    if (this->dim != 1) {
+    tug_assert(this->dim == 1, "This function is only available for 1D grids.");
-      throw std::invalid_argument(
+    tug_assert(index < this->cols, "Index is out of bounds.");
-          "This function is only available for 1D grids.");
-    }
-    if (index >= this->cols) {
-      throw std::invalid_argument("Index is out of bounds.");
-    }
 
     auto it = this->inner_boundary.find(std::make_pair(0, index));
     if (it == this->inner_boundary.end()) {
@@ -445,13 +436,8 @@ public:
    */
   std::pair<bool, T> getInnerBoundary(std::uint32_t row,
                                       std::uint32_t col) const {
-    if (this->dim != 2) {
+    tug_assert(this->dim == 2, "This function is only available for 2D grids.");
-      throw std::invalid_argument(
+    tug_assert(row < this->rows && col < this->cols, "Index is out of bounds.");
-          "This function is only available for 2D grids.");
-    }
-    if (row >= this->rows || col >= this->cols) {
-      throw std::invalid_argument("Index is out of bounds.");
-    }
 
     auto it = this->inner_boundary.find(std::make_pair(row, col));
     if (it == this->inner_boundary.end()) {
@@ -471,9 +457,7 @@ public:
    * condition
    */
   std::vector<std::pair<bool, T>> getInnerBoundaryRow(std::uint32_t row) const {
-    if (row >= this->rows) {
+    tug_assert(row < this->rows, "Index is out of bounds.");
-      throw std::invalid_argument("Index is out of bounds.");
-    }
 
     if (inner_boundary.empty()) {
       return std::vector<std::pair<bool, T>>(this->cols,
@@ -499,14 +483,8 @@ public:
    * condition
    */
   std::vector<std::pair<bool, T>> getInnerBoundaryCol(std::uint32_t col) const {
-    if (this->dim != 2) {
+    tug_assert(this->dim == 2, "This function is only available for 2D grids.");
-      throw std::invalid_argument(
+    tug_assert(col < this->cols, "Index is out of bounds.");
-          "This function is only available for 2D grids.");
-    }
-
-    if (col >= this->cols) {
-      throw std::invalid_argument("Index is out of bounds.");
-    }
 
     if (inner_boundary.empty()) {
       return std::vector<std::pair<bool, T>>(this->rows,

include/tug/Core/BaseSimulation.hpp

@@ -0,0 +1,388 @@
+#pragma once
+
+#include "tug/Boundary.hpp"
+#include <cstddef>
+#include <cstdint>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Core/TugUtils.hpp>
+
+namespace tug {
+
+/**
+ * @brief Enum holding different options for .csv output.
+ *
+ */
+enum class CSV_OUTPUT {
+  OFF,     /*!< do not produce csv output */
+  ON,      /*!< produce csv output with last concentration matrix */
+  VERBOSE, /*!< produce csv output with all concentration matrices */
+  XTREME   /*!< csv output like VERBOSE but additional boundary
+                         conditions at beginning */
+};
+
+/**
+ * @brief Enum holding different options for console output.
+ *
+ */
+enum class CONSOLE_OUTPUT {
+  OFF,    /*!< do not print any output to console */
+  ON,     /*!< print before and after concentrations to console */
+  VERBOSE /*!< print all concentration matrices to console */
+};
+
+/**
+ * @brief Enum holding different options for time measurement.
+ *
+ */
+enum class TIME_MEASURE {
+  OFF, /*!< do not print any time measures */
+  ON   /*!< print time measure after last iteration */
+};
+
+/**
+ * @brief A base class for simulation grids.
+ *
+ * This class provides a base implementation for simulation grids, including
+ * methods for setting and getting grid dimensions, domain sizes, and output
+ * options. It also includes methods for running simulations, which must be
+ * implemented by derived classes.
+ *
+ * @tparam T The type of the elements in the grid.
+ */
+template <typename T> class BaseSimulationGrid {
+private:
+  CSV_OUTPUT csv_output{CSV_OUTPUT::OFF};
+  CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT::OFF};
+  TIME_MEASURE time_measure{TIME_MEASURE::OFF};
+
+  int iterations{1};
+  RowMajMatMap<T> concentrationMatrix;
+  Boundary<T> boundaryConditions;
+
+  const std::uint8_t dim;
+
+  T delta_col;
+  T delta_row;
+
+public:
+  /**
+   * @brief Constructs a BaseSimulationGrid from a given RowMajMat object.
+   *
+   * This constructor initializes a BaseSimulationGrid using the data, number of
+   * rows, and number of columns from the provided RowMajMat object.
+   *
+   * @tparam T The type of the elements in the RowMajMat.
+   * @param origin The RowMajMat object from which to initialize the
+   * BaseSimulationGrid.
+   */
+  BaseSimulationGrid(RowMajMat<T> &origin)
+      : BaseSimulationGrid(origin.data(), origin.rows(), origin.cols()) {}
+
+  /**
+   * @brief Constructs a BaseSimulationGrid object.
+   *
+   * @tparam T The type of the data elements.
+   * @param data Pointer to the data array.
+   * @param rows Number of rows in the grid.
+   * @param cols Number of columns in the grid.
+   *
+   * Initializes the concentration_matrix with the provided data, rows, and
+   * columns. Sets delta_col and delta_row to 1. Determines the dimension (dim)
+   * based on the number of rows: if rows == 1, dim is set to 1; otherwise, dim
+   * is set to 2.
+   */
+  BaseSimulationGrid(T *data, std::size_t rows, std::size_t cols)
+      : concentrationMatrix(data, rows, cols), boundaryConditions(rows, cols),
+        delta_col(1), delta_row(1), dim(rows == 1 ? 1 : 2) {}
+
+  /**
+   * @brief Constructs a BaseSimulationGrid with a single dimension.
+   *
+   * This constructor initializes a BaseSimulationGrid object with the provided
+   * data and length. It assumes the grid has only one dimension.
+   *
+   * @param data Pointer to the data array.
+   * @param length The length of the data array.
+   */
+  BaseSimulationGrid(T *data, std::size_t length)
+      : BaseSimulationGrid(data, 1, length) {}
+
+  /**
+   * @brief Overloaded function call operator to access elements in a
+   * one-dimensional grid.
+   *
+   * This operator provides access to elements in the concentration matrix using
+   * a single index. It asserts that the grid is one-dimensional before
+   * accessing the element.
+   *
+   * @tparam T The type of elements in the concentration matrix.
+   * @param index The index of the element to access.
+   * @return A reference to the element at the specified index in the
+   * concentration matrix.
+   */
+  constexpr T &operator()(std::size_t index) {
+    tug_assert(dim == 1, "Grid is not one dimensional, use 2D index operator!");
+
+    return concentrationMatrix(index);
+  }
+
+  /**
+   * @brief Overloaded function call operator to access elements in a 2D
+   * concentration matrix.
+   *
+   * This operator allows accessing elements in the concentration matrix using
+   * row and column indices. It asserts that the grid is two-dimensional before
+   * accessing the element.
+   *
+   * @param row The row index of the element to access.
+   * @param col The column index of the element to access.
+   * @return A reference to the element at the specified row and column in the
+   * concentration matrix.
+   */
+  constexpr T &operator()(std::size_t row, std::size_t col) {
+    tug_assert(dim == 2, "Grid is not two dimensional, use 1D index operator!");
+
+    return concentrationMatrix(row, col);
+  }
+
+  /**
+   * @brief Retrieves the concentration matrix.
+   *
+   * @tparam T The data type of the elements in the concentration matrix.
+   * @return RowMajMat<T>& Reference to the concentration matrix.
+   */
+  RowMajMatMap<T> &getConcentrationMatrix() { return concentrationMatrix; }
+
+  const RowMajMatMap<T> &getConcentrationMatrix() const {
+    return concentrationMatrix;
+  }
+
+  /**
+   * @brief Retrieves the boundary conditions for the simulation.
+   *
+   * @tparam T The type parameter for the Boundary class.
+   * @return Boundary<T>& A reference to the boundary conditions.
+   */
+  Boundary<T> &getBoundaryConditions() { return boundaryConditions; }
+
+  const Boundary<T> &getBoundaryConditions() const {
+    return boundaryConditions;
+  }
+
+  /**
+   * @brief Retrieves the dimension value.
+   *
+   * @return The dimension value as an 8-bit unsigned integer.
+   */
+  std::uint8_t getDim() const { return dim; }
+
+  /**
+   * @brief Returns the number of rows in the concentration matrix.
+   *
+   * @return std::size_t The number of rows in the concentration matrix.
+   */
+  std::size_t rows() const { return concentrationMatrix.rows(); }
+
+  /**
+   * @brief Get the number of columns in the concentration matrix.
+   *
+   * @return std::size_t The number of columns in the concentration matrix.
+   */
+  std::size_t cols() const { return concentrationMatrix.cols(); }
+
+  /**
+   * @brief Returns the cell size in meter of the x-direction.
+   *
+   * This function returns the value of the delta column, which is used
+   * to represent the difference or change in the column value.
+   *
+   * @return T The cell size in meter of the x-direction.
+   */
+  T deltaCol() const { return delta_col; }
+
+  /**
+   * @brief Returns the cell size in meter of the y-direction.
+   *
+   * This function asserts that the grid is two-dimensional. If the grid is not
+   * two-dimensional, an assertion error is raised with the message "Grid is not
+   * two dimensional, there is no delta in y-direction!".
+   *
+   * @return The cell size in meter of the y-direction.
+   */
+  T deltaRow() const {
+    tug_assert(
+        dim == 2,
+        "Grid is not two dimensional, there is no delta in y-direction!");
+
+    return delta_row;
+  }
+
+  /**
+   * @brief Computes the domain size in the X direction.
+   *
+   * This function calculates the size of the domain in the X direction by
+   * multiplying the column spacing (delta_col) by the number of columns (cols).
+   *
+   * @return The size of the domain in the X direction.
+   */
+  T domainX() const { return delta_col * cols(); }
+
+  /**
+   * @brief Returns the size of the domain in the y-direction.
+   *
+   * This function calculates the size of the domain in the y-direction
+   * by multiplying the row spacing (delta_row) by the number of rows.
+   * It asserts that the grid is two-dimensional before performing the
+   * calculation.
+   *
+   * @return The size of the domain in the y-direction.
+   */
+  T domainY() const {
+    tug_assert(
+        dim == 2,
+        "Grid is not two dimensional, there is no domain in y-direction!");
+
+    return delta_row * rows();
+  }
+
+  /**
+   * @brief Sets the domain length for a one-dimensional grid.
+   *
+   * This function sets the domain length for a one-dimensional grid and
+   * calculates the column width (delta_col) based on the given domain length
+   * and the number of columns. It asserts that the grid is one-dimensional and
+   * that the given domain length is positive.
+   *
+   * @param domain_length The length of the domain. Must be positive.
+   */
+  void setDomain(T domain_length) {
+    tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!");
+    tug_assert(domain_length > 0, "Given domain length is not positive!");
+
+    delta_col = domain_length / cols();
+  }
+
+  /**
+   * @brief Sets the domain size for a 2D grid simulation.
+   *
+   * This function sets the domain size in the x and y directions for a
+   * two-dimensional grid simulation. It asserts that the grid is indeed
+   * two-dimensional and that the provided domain sizes are positive.
+   *
+   * @tparam T The type of the domain size parameters.
+   * @param domain_row The size of the domain in the y-direction.
+   * @param domain_col The size of the domain in the x-direction.
+   */
+  void setDomain(T domain_row, T domain_col) {
+    tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!");
+    tug_assert(domain_col > 0,
+               "Given domain size in x-direction is not positive!");
+    tug_assert(domain_row > 0,
+               "Given domain size in y-direction is not positive!");
+
+    delta_row = domain_row / rows();
+    delta_col = domain_col / cols();
+  }
+
+  /**
+   * @brief Set the option to output the results to a CSV file. Off by default.
+   *
+   *
+   * @param csv_output Valid output option. The following options can be set
+   *                   here:
+   *                     - CSV_OUTPUT_OFF: do not produce csv output
+   *                     - CSV_OUTPUT_ON: produce csv output with last
+   *                       concentration matrix
+   *                     - CSV_OUTPUT_VERBOSE: produce csv output with all
+   *                       concentration matrices
+   *                     - CSV_OUTPUT_XTREME: produce csv output with all
+   *                       concentration matrices and simulation environment
+   */
+  void setOutputCSV(CSV_OUTPUT csv_output) { this->csv_output = csv_output; }
+
+  /**
+   * @brief Retrieves the CSV output.
+   *
+   * This function returns the CSV output associated with the simulation.
+   *
+   * @return CSV_OUTPUT The CSV output of the simulation.
+   */
+  constexpr CSV_OUTPUT getOutputCSV() const { return this->csv_output; }
+
+  /**
+   * @brief Set the options for outputting information to the console. Off by
+   * default.
+   *
+   * @param console_output Valid output option. The following options can be set
+   *                       here:
+   *                        - CONSOLE_OUTPUT_OFF: do not print any output to
+   * console
+   *                        - CONSOLE_OUTPUT_ON: print before and after
+   * concentrations to console
+   *                        - CONSOLE_OUTPUT_VERBOSE: print all concentration
+   * matrices to console
+   */
+  void setOutputConsole(CONSOLE_OUTPUT console_output) {
+    this->console_output = console_output;
+  }
+
+  /**
+   * @brief Retrieves the console output.
+   *
+   * This function returns the current state of the console output.
+   *
+   * @return CONSOLE_OUTPUT The current console output.
+   */
+  constexpr CONSOLE_OUTPUT getOutputConsole() const {
+    return this->console_output;
+  }
+
+  /**
+   * @brief Set the Time Measure option. Off by default.
+   *
+   * @param time_measure The following options are allowed:
+   *                     - TIME_MEASURE_OFF: Time of simulation is not printed
+   * to console
+   *                     - TIME_MEASURE_ON: Time of simulation run is printed to
+   * console
+   */
+  void setTimeMeasure(TIME_MEASURE time_measure) {
+    this->time_measure = time_measure;
+  }
+
+  /**
+   * @brief Retrieves the current time measurement.
+   *
+   * @return TIME_MEASURE The current time measurement.
+   */
+  constexpr TIME_MEASURE getTimeMeasure() const { return this->time_measure; }
+
+  /**
+   * @brief Set the desired iterations to be calculated. A value greater
+   *        than zero must be specified here. Setting iterations is required.
+   *
+   * @param iterations Number of iterations to be simulated.
+   */
+  void setIterations(int iterations) {
+    tug_assert(iterations > 0,
+               "Number of iterations must be greater than zero.");
+
+    this->iterations = iterations;
+  }
+
+  /**
+   * @brief Return the currently set iterations to be calculated.
+   *
+   * @return int Number of iterations.
+   */
+  int getIterations() const { return this->iterations; }
+
+  /**
+   * @brief Method starts the simulation process with the previously set
+   *        parameters.
+   */
+  virtual void run() = 0;
+
+  virtual void setTimestep(T timestep) = 0;
+};
+} // namespace tug

include/tug/Core/FTCS.hpp

@@ -1,402 +0,0 @@
-/**
- * @file FTCS.hpp
- * @brief Implementation of heterogenous FTCS (forward time-centered space)
- * solution of diffusion equation in 1D and 2D space.
- *
- */
-
-#ifndef FTCS_H_
-#define FTCS_H_
-
-#include "TugUtils.hpp"
-
-#include <cstddef>
-#include <iostream>
-#include <tug/Boundary.hpp>
-
-#ifdef _OPENMP
-#include <omp.h>
-#else
-#define omp_get_thread_num() 0
-#endif
-
-namespace tug {
-
-// calculates horizontal change on one cell independent of boundary type
-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)) *
-             grid.getConcentrations()(row, col + 1) -
-         (calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
-                             grid.getAlphaX()(row, col)) +
-          calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
-                             grid.getAlphaX()(row, col))) *
-             grid.getConcentrations()(row, col) +
-         calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
-                            grid.getAlphaX()(row, col)) *
-             grid.getConcentrations()(row, col - 1);
-}
-
-// calculates vertical change on one cell independent of boundary type
-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)) *
-             grid.getConcentrations()(row + 1, col) -
-         (calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
-                             grid.getAlphaY()(row, col)) +
-          calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
-                             grid.getAlphaY()(row, col))) *
-             grid.getConcentrations()(row, col) +
-         calcAlphaIntercell(grid.getAlphaY()(row - 1, col),
-                            grid.getAlphaY()(row, col)) *
-             grid.getConcentrations()(row - 1, col);
-}
-
-// calculates horizontal change on one cell with a constant left boundary
-template <class T>
-static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid,
-                                                         Boundary<T> &bc,
-                                                         int &row, int &col) {
-
-  return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
-                            grid.getAlphaX()(row, col)) *
-             grid.getConcentrations()(row, col + 1) -
-         (calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
-                             grid.getAlphaX()(row, col)) +
-          2 * grid.getAlphaX()(row, col)) *
-             grid.getConcentrations()(row, col) +
-         2 * grid.getAlphaX()(row, col) *
-             bc.getBoundaryElementValue(BC_SIDE_LEFT, row);
-}
-
-// calculates horizontal change on one cell with a closed left boundary
-template <class T>
-static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &grid, int &row,
-                                                       int &col) {
-
-  return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
-                            grid.getAlphaX()(row, col)) *
-         (grid.getConcentrations()(row, col + 1) -
-          grid.getConcentrations()(row, col));
-}
-
-// checks boundary condition type for a cell on the left edge of grid
-template <class T>
-static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &bc,
-                                                 int &row, int &col) {
-  if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) {
-    return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
-  } else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED) {
-    return calcHorizontalChangeLeftBoundaryClosed(grid, row, col);
-  } else {
-    throw_invalid_argument("Undefined Boundary Condition Type!");
-  }
-}
-
-// calculates horizontal change on one cell with a constant right boundary
-template <class T>
-static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid,
-                                                          Boundary<T> &bc,
-                                                          int &row, int &col) {
-
-  return 2 * grid.getAlphaX()(row, col) *
-             bc.getBoundaryElementValue(BC_SIDE_RIGHT, row) -
-         (calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
-                             grid.getAlphaX()(row, col)) +
-          2 * grid.getAlphaX()(row, col)) *
-             grid.getConcentrations()(row, col) +
-         calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
-                            grid.getAlphaX()(row, col)) *
-             grid.getConcentrations()(row, col - 1);
-}
-
-// calculates horizontal change on one cell with a closed right boundary
-template <class T>
-static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &grid, int &row,
-                                                        int &col) {
-
-  return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
-                              grid.getAlphaX()(row, col)) *
-           (grid.getConcentrations()(row, col) -
-            grid.getConcentrations()(row, col - 1)));
-}
-
-// checks boundary condition type for a cell on the right edge of grid
-template <class T>
-static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid,
-                                                  Boundary<T> &bc, int &row,
-                                                  int &col) {
-  if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) {
-    return calcHorizontalChangeRightBoundaryConstant(grid, bc, row, col);
-  } else if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CLOSED) {
-    return calcHorizontalChangeRightBoundaryClosed(grid, row, col);
-  } else {
-    throw_invalid_argument("Undefined Boundary Condition Type!");
-  }
-}
-
-// calculates vertical change on one cell with a constant top boundary
-template <class T>
-static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid,
-                                                      Boundary<T> &bc, int &row,
-                                                      int &col) {
-
-  return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
-                            grid.getAlphaY()(row, col)) *
-             grid.getConcentrations()(row + 1, col) -
-         (calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
-                             grid.getAlphaY()(row, col)) +
-          2 * grid.getAlphaY()(row, col)) *
-             grid.getConcentrations()(row, col) +
-         2 * grid.getAlphaY()(row, col) *
-             bc.getBoundaryElementValue(BC_SIDE_TOP, col);
-}
-
-// calculates vertical change on one cell with a closed top boundary
-template <class T>
-static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &grid, int &row,
-                                                    int &col) {
-
-  return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
-                            grid.getAlphaY()(row, col)) *
-         (grid.getConcentrations()(row + 1, col) -
-          grid.getConcentrations()(row, col));
-}
-
-// checks boundary condition type for a cell on the top edge of grid
-template <class T>
-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) {
-    return calcVerticalChangeTopBoundaryClosed(grid, row, col);
-  } else {
-    throw_invalid_argument("Undefined Boundary Condition Type!");
-  }
-}
-
-// calculates vertical change on one cell with a constant bottom boundary
-template <class T>
-static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid,
-                                                         Boundary<T> &bc,
-                                                         int &row, int &col) {
-
-  return 2 * grid.getAlphaY()(row, col) *
-             bc.getBoundaryElementValue(BC_SIDE_BOTTOM, col) -
-         (calcAlphaIntercell(grid.getAlphaY()(row, col),
-                             grid.getAlphaY()(row - 1, col)) +
-          2 * grid.getAlphaY()(row, col)) *
-             grid.getConcentrations()(row, col) +
-         calcAlphaIntercell(grid.getAlphaY()(row, col),
-                            grid.getAlphaY()(row - 1, col)) *
-             grid.getConcentrations()(row - 1, col);
-}
-
-// calculates vertical change on one cell with a closed bottom boundary
-template <class T>
-static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &grid, int &row,
-                                                       int &col) {
-
-  return -(calcAlphaIntercell(grid.getAlphaY()(row, col),
-                              grid.getAlphaY()(row - 1, col)) *
-           (grid.getConcentrations()(row, col) -
-            grid.getConcentrations()(row - 1, col)));
-}
-
-// checks boundary condition type for a cell on the bottom edge of grid
-template <class T>
-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) {
-    return calcVerticalChangeBottomBoundaryClosed(grid, row, col);
-  } else {
-    throw_invalid_argument("Undefined Boundary Condition Type!");
-  }
-}
-
-// FTCS solution for 1D grid
-template <class T>
-static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
-  int colMax = grid.getCol();
-  T deltaCol = grid.getDeltaCol();
-
-  // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(1, colMax, 0);
-
-  // only one row in 1D case -> row constant at index 0
-  int row = 0;
-
-  // inner cells
-  // independent of boundary condition type
-  for (int col = 1; col < colMax - 1; col++) {
-    concentrations_t1(row, col) = grid.getConcentrations()(row, col) +
-                                  timestep / (deltaCol * deltaCol) *
-                                      (calcHorizontalChange(grid, row, col));
-  }
-
-  // left boundary; hold column constant at index 0
-  int col = 0;
-  concentrations_t1(row, col) =
-      grid.getConcentrations()(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeLeftBoundary(grid, bc, row, col));
-
-  // right boundary; hold column constant at max index
-  col = colMax - 1;
-  concentrations_t1(row, col) =
-      grid.getConcentrations()(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeRightBoundary(grid, bc, row, col));
-
-  // overwrite obsolete concentrations
-  grid.setConcentrations(concentrations_t1);
-}
-
-// FTCS solution for 2D grid
-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();
-  T deltaRow = grid.getDeltaRow();
-  T deltaCol = grid.getDeltaCol();
-
-  // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
-
-  // inner cells
-  // these are independent of the boundary condition type
-  // omp_set_num_threads(10);
-#pragma omp parallel for num_threads(numThreads)
-  for (int row = 1; row < rowMax - 1; row++) {
-    for (int col = 1; col < colMax - 1; col++) {
-      concentrations_t1(row, col) = grid.getConcentrations()(row, col) +
-                                    timestep / (deltaRow * deltaRow) *
-                                        (calcVerticalChange(grid, row, col)) +
-                                    timestep / (deltaCol * deltaCol) *
-                                        (calcHorizontalChange(grid, row, col));
-    }
-  }
-
-  // boundary conditions
-  // left without corners / looping over rows
-  // hold column constant at index 0
-  int col = 0;
-#pragma omp parallel for num_threads(numThreads)
-  for (int row = 1; row < rowMax - 1; row++) {
-    concentrations_t1(row, col) =
-        grid.getConcentrations()(row, col) +
-        timestep / (deltaCol * deltaCol) *
-            (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
-        timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
-  }
-
-  // right without corners / looping over rows
-  // hold column constant at max index
-  col = colMax - 1;
-#pragma omp parallel for num_threads(numThreads)
-  for (int row = 1; row < rowMax - 1; row++) {
-    concentrations_t1(row, col) =
-        grid.getConcentrations()(row, col) +
-        timestep / (deltaCol * deltaCol) *
-            (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
-        timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col));
-  }
-
-  // top without corners / looping over columns
-  // hold row constant at index 0
-  int row = 0;
-#pragma omp parallel for num_threads(numThreads)
-  for (int col = 1; col < colMax - 1; col++) {
-    concentrations_t1(row, col) =
-        grid.getConcentrations()(row, col) +
-        timestep / (deltaRow * deltaRow) *
-            (calcVerticalChangeTopBoundary(grid, bc, row, col)) +
-        timestep / (deltaCol * deltaCol) *
-            (calcHorizontalChange(grid, row, col));
-  }
-
-  // bottom without corners / looping over columns
-  // hold row constant at max index
-  row = rowMax - 1;
-#pragma omp parallel for num_threads(numThreads)
-  for (int col = 1; col < colMax - 1; col++) {
-    concentrations_t1(row, col) =
-        grid.getConcentrations()(row, col) +
-        timestep / (deltaRow * deltaRow) *
-            (calcVerticalChangeBottomBoundary(grid, bc, row, col)) +
-        timestep / (deltaCol * deltaCol) *
-            (calcHorizontalChange(grid, row, col));
-  }
-
-  // corner top left
-  // hold row and column constant at 0
-  row = 0;
-  col = 0;
-  concentrations_t1(row, col) =
-      grid.getConcentrations()(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
-      timestep / (deltaRow * deltaRow) *
-          (calcVerticalChangeTopBoundary(grid, bc, row, col));
-
-  // corner top right
-  // hold row constant at 0 and column constant at max index
-  row = 0;
-  col = colMax - 1;
-  concentrations_t1(row, col) =
-      grid.getConcentrations()(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
-      timestep / (deltaRow * deltaRow) *
-          (calcVerticalChangeTopBoundary(grid, bc, row, col));
-
-  // corner bottom left
-  // hold row constant at max index and column constant at 0
-  row = rowMax - 1;
-  col = 0;
-  concentrations_t1(row, col) =
-      grid.getConcentrations()(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) +
-      timestep / (deltaRow * deltaRow) *
-          (calcVerticalChangeBottomBoundary(grid, bc, row, col));
-
-  // corner bottom right
-  // hold row and column constant at max index
-  row = rowMax - 1;
-  col = colMax - 1;
-  concentrations_t1(row, col) =
-      grid.getConcentrations()(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
-      timestep / (deltaRow * deltaRow) *
-          (calcVerticalChangeBottomBoundary(grid, bc, row, col));
-
-  // overwrite obsolete concentrations
-  grid.setConcentrations(concentrations_t1);
-  // }
-}
-
-// entry point; differentiate between 1D and 2D grid
-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) {
-    FTCS_2D(grid, bc, timestep, numThreads);
-  } else {
-    throw_invalid_argument(
-        "Error: Only 1- and 2-dimensional grids are defined!");
-  }
-}
-} // namespace tug
-
-#endif // FTCS_H_

include/tug/Core/Numeric/BTCS.hpp

@@ -10,15 +10,16 @@
 #ifndef BTCS_H_
 #define BTCS_H_
 
-#include "Matrix.hpp"
-#include "TugUtils.hpp"
-
 #include <cstddef>
 #include <tug/Boundary.hpp>
-#include <tug/Grid.hpp>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Core/Numeric/SimulationInput.hpp>
 #include <utility>
 #include <vector>
 
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
 #ifdef _OPENMP
 #include <omp.h>
 #else
@@ -159,9 +160,9 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc,
 
 // creates a solution vector for next time step from the current state of
 // concentrations
-template <class T>
+template <class T, class EigenType>
 static Eigen::VectorX<T>
-createSolutionVector(const RowMajMat<T> &concentrations,
+createSolutionVector(const EigenType &concentrations,
                      const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY,
                      const std::vector<BoundaryElement<T>> &bcLeft,
                      const std::vector<BoundaryElement<T>> &bcRight,
@@ -351,25 +352,27 @@ static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
 
 // BTCS solution for 1D grid
 template <class T>
-static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
+static void BTCS_1D(SimulationInput<T> &input,
                     Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
                                                     Eigen::VectorX<T> &b)) {
-  int length = grid.getLength();
+  const std::size_t &length = input.colMax;
-  T sx = timestep / (grid.getDelta() * grid.getDelta());
+  T sx = input.timestep / (input.deltaCol * input.deltaCol);
 
   Eigen::VectorX<T> concentrations_t1(length);
 
   Eigen::SparseMatrix<T> A;
   Eigen::VectorX<T> b(length);
 
-  const auto &alpha = grid.getAlpha();
+  const auto &alpha = input.alphaX;
+
+  const auto &bc = input.boundaries;
 
   const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
   const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
   const auto inner_bc = bc.getInnerBoundaryRow(0);
 
-  RowMajMat<T> &concentrations = grid.getConcentrations();
+  RowMajMatMap<T> &concentrations = input.concentrations;
   int rowIndex = 0;
   A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex,
                         sx); // this is exactly same as in 2D
@@ -396,29 +399,31 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
 
 // BTCS solution for 2D grid
 template <class T>
-static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
+static void BTCS_2D(SimulationInput<T> &input,
                     Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
                                                     Eigen::VectorX<T> &b),
                     int numThreads) {
-  int rowMax = grid.getRow();
+  const std::size_t &rowMax = input.rowMax;
-  int colMax = grid.getCol();
+  const std::size_t &colMax = input.colMax;
-  T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
+  const T sx = input.timestep / (2 * input.deltaCol * input.deltaCol);
-  T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
+  const T sy = input.timestep / (2 * input.deltaRow * input.deltaRow);
 
   RowMajMat<T> concentrations_t1(rowMax, colMax);
 
   Eigen::SparseMatrix<T> A;
   Eigen::VectorX<T> b;
 
-  RowMajMat<T> alphaX = grid.getAlphaX();
+  const RowMajMat<T> &alphaX = input.alphaX;
-  RowMajMat<T> alphaY = grid.getAlphaY();
+  const RowMajMat<T> &alphaY = input.alphaY;
+
+  const auto &bc = input.boundaries;
 
   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);
 
-  RowMajMat<T> &concentrations = grid.getConcentrations();
+  RowMajMatMap<T> &concentrations = input.concentrations;
 
 #pragma omp parallel for num_threads(numThreads) private(A, b)
   for (int i = 0; i < rowMax; i++) {
@@ -432,44 +437,43 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
   }
 
   concentrations_t1.transposeInPlace();
-  alphaX.transposeInPlace();
+  const RowMajMat<T> alphaX_t = alphaX.transpose();
-  alphaY.transposeInPlace();
+  const RowMajMat<T> alphaY_t = alphaY.transpose();
 
 #pragma omp parallel for num_threads(numThreads) private(A, b)
   for (int i = 0; i < colMax; i++) {
     auto inner_bc = bc.getInnerBoundaryCol(i);
     // swap alphas, boundary conditions and sx/sy for column-wise calculation
-    A = createCoeffMatrix(alphaY, bcTop, bcBottom, inner_bc, rowMax, i, sy);
+    A = createCoeffMatrix(alphaY_t, bcTop, bcBottom, inner_bc, rowMax, i, sy);
-    b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
+    b = createSolutionVector(concentrations_t1, alphaY_t, alphaX_t, bcTop,
-                             bcLeft, bcRight, inner_bc, rowMax, i, sy, sx);
+                             bcBottom, bcLeft, bcRight, inner_bc, rowMax, i, sy,
+                             sx);
 
     concentrations.col(i) = solverFunc(A, b);
   }
 }
 
 // entry point for EigenLU solver; differentiate between 1D and 2D grid
-template <class T>
+template <class T> void BTCS_LU(SimulationInput<T> &input, int numThreads) {
-void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
+  tug_assert(input.dim <= 2,
-  if (grid.getDim() == 1) {
-    BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
-  } else if (grid.getDim() == 2) {
-    BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads);
-  } else {
-    throw_invalid_argument(
              "Error: Only 1- and 2-dimensional grids are defined!");
+
+  if (input.dim == 1) {
+    BTCS_1D(input, EigenLUAlgorithm);
+  } else {
+    BTCS_2D(input.dim, EigenLUAlgorithm, numThreads);
   }
 }
 
 // entry point for Thomas algorithm solver; differentiate 1D and 2D grid
-template <class T>
+template <class T> void BTCS_Thomas(SimulationInput<T> &input, int numThreads) {
-void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
+  tug_assert(input.dim <= 2,
-  if (grid.getDim() == 1) {
-    BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
-  } else if (grid.getDim() == 2) {
-    BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads);
-  } else {
-    throw_invalid_argument(
              "Error: Only 1- and 2-dimensional grids are defined!");
+
+  if (input.dim == 1) {
+    BTCS_1D(input, ThomasAlgorithm);
+  } else {
+    BTCS_2D(input, ThomasAlgorithm, numThreads);
   }
 }
 } // namespace tug

include/tug/Core/Numeric/FTCS.hpp

@@ -0,0 +1,240 @@
+/**
+ * @file FTCS.hpp
+ * @brief Implementation of heterogenous FTCS (forward time-centered space)
+ * solution of diffusion equation in 1D and 2D space.
+ *
+ */
+
+#ifndef FTCS_H_
+#define FTCS_H_
+
+#include "tug/Core/TugUtils.hpp"
+#include <cstddef>
+#include <cstring>
+#include <tug/Boundary.hpp>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Core/Numeric/SimulationInput.hpp>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_thread_num() 0
+#endif
+
+namespace tug {
+
+template <class T>
+constexpr T calcChangeInner(T conc_c, T conc_left, T conc_right, T alpha_c,
+                            T alpha_left, T alpha_right) {
+  const T alpha_center_left = calcAlphaIntercell(alpha_left, alpha_c);
+  const T alpha_center_right = calcAlphaIntercell(alpha_right, alpha_c);
+
+  return alpha_center_left * conc_left -
+         (alpha_center_left + alpha_center_right) * conc_c +
+         alpha_center_right * conc_right;
+}
+
+template <class T>
+constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center,
+                               T alpha_neighbor, const BoundaryElement<T> &bc) {
+  const T alpha_center_neighbor =
+      calcAlphaIntercell(alpha_center, alpha_neighbor);
+  const T &conc_boundary = bc.getValue();
+
+  switch (bc.getType()) {
+  case BC_TYPE_CONSTANT: {
+    return 2 * alpha_center * conc_boundary -
+           (alpha_center_neighbor + 2 * alpha_center) * conc_c +
+           alpha_center_neighbor * conc_neighbor;
+  }
+  case BC_TYPE_CLOSED: {
+    return (alpha_center_neighbor * (conc_neighbor - conc_c));
+  }
+  }
+
+  tug_assert(false, "Undefined Boundary Condition Type!");
+}
+
+// FTCS solution for 1D grid
+template <class T> static void FTCS_1D(SimulationInput<T> &input) {
+  const std::size_t &colMax = input.colMax;
+  const T &deltaCol = input.deltaCol;
+  const T &timestep = input.timestep;
+
+  RowMajMatMap<T> &concentrations_grid = input.concentrations;
+  // matrix for concentrations at time t+1
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
+
+  const auto &alphaX = input.alphaX;
+  const auto &bc = input.boundaries;
+
+  // only one row in 1D case -> row constant at index 0
+  int row = 0;
+
+  // inner cells
+  // independent of boundary condition type
+  for (int col = 1; col < colMax - 1; col++) {
+    const T &conc_c = concentrations_grid(row, col);
+    const T &conc_left = concentrations_grid(row, col - 1);
+    const T &conc_right = concentrations_grid(row, col + 1);
+
+    const T &alpha_c = alphaX(row, col);
+    const T &alpha_left = alphaX(row, col - 1);
+    const T &alpha_right = alphaX(row, col + 1);
+
+    concentrations_t1(row, col) =
+        concentrations_grid(row, col) +
+        timestep / (deltaCol * deltaCol) *
+            calcChangeInner(conc_c, conc_left, conc_right, alpha_c, alpha_left,
+                            alpha_right);
+  }
+
+  // left boundary; hold column constant at index 0
+  {
+    int col = 0;
+    const T &conc_c = concentrations_grid(row, col);
+    const T &conc_right = concentrations_grid(row, col + 1);
+    const T &alpha_c = alphaX(row, col);
+    const T &alpha_right = alphaX(row, col + 1);
+    const BoundaryElement<T> &bc_element =
+        input.boundaries.getBoundaryElement(BC_SIDE_LEFT, row);
+
+    concentrations_t1(row, col) =
+        concentrations_grid(row, col) +
+        timestep / (deltaCol * deltaCol) *
+            calcChangeBoundary(conc_c, conc_right, alpha_c, alpha_right,
+                               bc_element);
+  }
+
+  // right boundary; hold column constant at max index
+  {
+    int col = colMax - 1;
+    const T &conc_c = concentrations_grid(row, col);
+    const T &conc_left = concentrations_grid(row, col - 1);
+    const T &alpha_c = alphaX(row, col);
+    const T &alpha_left = alphaX(row, col - 1);
+    const BoundaryElement<T> &bc_element =
+        bc.getBoundaryElement(BC_SIDE_RIGHT, row);
+
+    concentrations_t1(row, col) =
+        concentrations_grid(row, col) +
+        timestep / (deltaCol * deltaCol) *
+            calcChangeBoundary(conc_c, conc_left, alpha_c, alpha_left,
+                               bc_element);
+  }
+  // overwrite obsolete concentrations
+  concentrations_grid = concentrations_t1;
+}
+
+// FTCS solution for 2D grid
+template <class T>
+static void FTCS_2D(SimulationInput<T> &input, int numThreads) {
+  const std::size_t &rowMax = input.rowMax;
+  const std::size_t &colMax = input.colMax;
+  const T &deltaRow = input.deltaRow;
+  const T &deltaCol = input.deltaCol;
+  const T &timestep = input.timestep;
+
+  RowMajMatMap<T> &concentrations_grid = input.concentrations;
+
+  // matrix for concentrations at time t+1
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
+
+  const auto &alphaX = input.alphaX;
+  const auto &alphaY = input.alphaY;
+  const auto &bc = input.boundaries;
+
+  const T sx = timestep / (deltaCol * deltaCol);
+  const T sy = timestep / (deltaRow * deltaRow);
+
+#pragma omp parallel for num_threads(numThreads)
+  for (std::size_t row_i = 0; row_i < rowMax; row_i++) {
+    for (std::size_t col_i = 0; col_i < colMax; col_i++) {
+      // horizontal change
+      T horizontal_change;
+      {
+
+        const T &conc_c = concentrations_grid(row_i, col_i);
+        const T &alpha_c = alphaX(row_i, col_i);
+
+        if (col_i == 0 || col_i == colMax - 1) {
+          // left or right boundary
+          const T &conc_neigbor =
+              concentrations_grid(row_i, col_i == 0 ? col_i + 1 : col_i - 1);
+          const T &alpha_neigbor =
+              alphaX(row_i, col_i == 0 ? col_i + 1 : col_i - 1);
+
+          const BoundaryElement<T> &bc_element = bc.getBoundaryElement(
+              col_i == 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i);
+
+          horizontal_change = calcChangeBoundary(conc_c, conc_neigbor, alpha_c,
+                                                 alpha_neigbor, bc_element);
+        } else {
+          // inner cell
+          const T &conc_left = concentrations_grid(row_i, col_i - 1);
+          const T &conc_right = concentrations_grid(row_i, col_i + 1);
+
+          const T &alpha_left = alphaX(row_i, col_i - 1);
+          const T &alpha_right = alphaX(row_i, col_i + 1);
+
+          horizontal_change = calcChangeInner(conc_c, conc_left, conc_right,
+                                              alpha_c, alpha_left, alpha_right);
+        }
+      }
+
+      // vertical change
+      T vertical_change;
+      {
+        const T &conc_c = concentrations_grid(row_i, col_i);
+        const T &alpha_c = alphaY(row_i, col_i);
+
+        if (row_i == 0 || row_i == rowMax - 1) {
+          // top or bottom boundary
+          const T &conc_neigbor =
+              concentrations_grid(row_i == 0 ? row_i + 1 : row_i - 1, col_i);
+
+          const T &alpha_neigbor =
+              alphaY(row_i == 0 ? row_i + 1 : row_i - 1, col_i);
+
+          const BoundaryElement<T> &bc_element = bc.getBoundaryElement(
+              row_i == 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i);
+
+          vertical_change = calcChangeBoundary(conc_c, conc_neigbor, alpha_c,
+                                               alpha_neigbor, bc_element);
+        } else {
+          // inner cell
+          const T &conc_bottom = concentrations_grid(row_i - 1, col_i);
+          const T &conc_top = concentrations_grid(row_i + 1, col_i);
+
+          const T &alpha_bottom = alphaY(row_i - 1, col_i);
+          const T &alpha_top = alphaY(row_i + 1, col_i);
+
+          vertical_change = calcChangeInner(conc_c, conc_bottom, conc_top,
+                                            alpha_c, alpha_bottom, alpha_top);
+        }
+      }
+
+      concentrations_t1(row_i, col_i) = concentrations_grid(row_i, col_i) +
+                                        sx * horizontal_change +
+                                        sy * vertical_change;
+    }
+  }
+
+  // overwrite obsolete concentrations
+  concentrations_grid = concentrations_t1;
+}
+
+// entry point; differentiate between 1D and 2D grid
+template <class T> void FTCS(SimulationInput<T> &input, int &numThreads) {
+  tug_assert(input.dim <= 2,
+             "Error: Only 1- and 2-dimensional grids are defined!");
+
+  if (input.dim == 1) {
+    FTCS_1D(input);
+  } else {
+    FTCS_2D(input, numThreads);
+  }
+}
+} // namespace tug
+
+#endif // FTCS_H_

include/tug/Core/Numeric/SimulationInput.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <tug/Boundary.hpp>
+#include <tug/Core/Matrix.hpp>
+
+namespace tug {
+
+template <typename T> struct SimulationInput {
+  RowMajMatMap<T> &concentrations;
+  const RowMajMat<T> &alphaX;
+  const RowMajMat<T> &alphaY;
+  const Boundary<T> boundaries;
+
+  const std::uint8_t dim;
+  const T timestep;
+  const std::size_t rowMax;
+  const std::size_t colMax;
+  const T deltaRow;
+  const T deltaCol;
+};
+} // namespace tug

include/tug/Core/TugUtils.hpp

@@ -1,9 +1,6 @@
-#ifndef TUGUTILS_H_
+#pragma once
-#define TUGUTILS_H_
 
-#include <chrono>
+#include <cassert>
-#include <stdexcept>
-#include <string>
 
 #define throw_invalid_argument(msg)                                            \
   throw std::invalid_argument(std::string(__FILE__) + ":" +                    \
@@ -24,6 +21,8 @@
     duration.count();                                                          \
   })
 
+#define tug_assert(expr, msg) assert((expr) && msg)
+
 // calculates arithmetic or harmonic mean of alpha between two cells
 template <typename T>
 constexpr T calcAlphaIntercell(T alpha1, T alpha2, bool useHarmonic = true) {
@@ -38,4 +37,3 @@ constexpr T calcAlphaIntercell(T alpha1, T alpha2, bool useHarmonic = true) {
     return 0.5 * (alpha1 + alpha2);
   }
 }
-#endif // TUGUTILS_H_

include/tug/Diffusion.hpp

@@ -1,16 +1,14 @@
 /**
- * @file Simulation.hpp
+ * @file Diffusion.hpp
- * @brief API of Simulation class, that holds all information regarding a
+ * @brief API of Diffusion class, that holds all information regarding a
  * specific simulation run like its timestep, number of iterations and output
- * options. Simulation object also holds a predefined Grid and Boundary object.
+ * options. Diffusion object also holds a predefined Grid and Boundary object.
  *
  */
 
-#ifndef SIMULATION_H_
+#pragma once
-#define SIMULATION_H_
 
-#include "Boundary.hpp"
+#include "tug/Core/Matrix.hpp"
-#include "Grid.hpp"
 #include <algorithm>
 #include <filesystem>
 #include <fstream>
@@ -19,9 +17,9 @@
 #include <string>
 #include <vector>
 
-#include "Core/BTCS.hpp"
+#include <tug/Core/BaseSimulation.hpp>
-#include "Core/FTCS.hpp"
+#include <tug/Core/Numeric/BTCS.hpp>
-#include "Core/TugUtils.hpp"
+#include <tug/Core/Numeric/FTCS.hpp>
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -51,37 +49,6 @@ enum SOLVER {
                              tridiagonal matrices */
 };
 
-/**
- * @brief Enum holding different options for .csv output.
- *
- */
-enum CSV_OUTPUT {
-  CSV_OUTPUT_OFF,     /*!< do not produce csv output */
-  CSV_OUTPUT_ON,      /*!< produce csv output with last concentration matrix */
-  CSV_OUTPUT_VERBOSE, /*!< produce csv output with all concentration matrices */
-  CSV_OUTPUT_XTREME   /*!< csv output like VERBOSE but additional boundary
-                         conditions at beginning */
-};
-
-/**
- * @brief Enum holding different options for console output.
- *
- */
-enum CONSOLE_OUTPUT {
-  CONSOLE_OUTPUT_OFF, /*!< do not print any output to console */
-  CONSOLE_OUTPUT_ON,  /*!< print before and after concentrations to console */
-  CONSOLE_OUTPUT_VERBOSE /*!< print all concentration matrices to console */
-};
-
-/**
- * @brief Enum holding different options for time measurement.
- *
- */
-enum TIME_MEASURE {
-  TIME_MEASURE_OFF, /*!< do not print any time measures */
-  TIME_MEASURE_ON   /*!< print time measure after last iteration */
-};
-
 /**
  * @brief The class forms the interface for performing the diffusion simulations
  * and contains all the methods for controlling the desired parameters, such as
@@ -94,81 +61,94 @@ enum TIME_MEASURE {
  */
 template <class T, APPROACH approach = BTCS_APPROACH,
           SOLVER solver = THOMAS_ALGORITHM_SOLVER>
-class Simulation {
+class Diffusion : public BaseSimulationGrid<T> {
+private:
+  T timestep{-1};
+  int innerIterations{1};
+  int numThreads{omp_get_num_procs()};
+
+  RowMajMat<T> alphaX;
+  RowMajMat<T> alphaY;
+
+  const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
+
+  static constexpr T DEFAULT_ALPHA = 1E-8;
+
+  void init_alpha() {
+    this->alphaX =
+        RowMajMat<T>::Constant(this->rows(), this->cols(), DEFAULT_ALPHA);
+    if (this->getDim() == 2) {
+      this->alphaY =
+          RowMajMat<T>::Constant(this->rows(), this->cols(), DEFAULT_ALPHA);
+    }
+  }
+
 public:
   /**
-   * @brief Set up a simulation environment. The timestep and number of
+   * @brief Construct a new Diffusion object from a given Eigen matrix
-   * iterations must be set. For the BTCS approach, the Thomas algorithm is used
-   * as the default linear equation solver as this is faster for tridiagonal
-   *        matrices. CSV output, console output and time measure are off by
-   * default. Also, the number of cores is set to the maximum number of cores -1
-   * by default.
    *
-   * @param grid Valid grid object
-   * @param bc Valid boundary condition object
-   * @param approach Approach to solving the problem. Either FTCS or BTCS.
    */
-  Simulation(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
+  Diffusion(RowMajMat<T> &origin) : BaseSimulationGrid<T>(origin) {
-
+    init_alpha();
-  /**
-   * @brief Set the option to output the results to a CSV file. Off by default.
-   *
-   *
-   * @param csv_output Valid output option. The following options can be set
-   *                   here:
-   *                     - CSV_OUTPUT_OFF: do not produce csv output
-   *                     - CSV_OUTPUT_ON: produce csv output with last
-   *                       concentration matrix
-   *                     - CSV_OUTPUT_VERBOSE: produce csv output with all
-   *                       concentration matrices
-   *                     - CSV_OUTPUT_XTREME: produce csv output with all
-   *                       concentration matrices and simulation environment
-   */
-  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
+   * @brief Construct a new 2D Diffusion object from a given data pointer and
-   * default.
+   * the dimensions.
    *
-   * @param console_output Valid output option. The following options can be set
-   *                       here:
-   *                        - CONSOLE_OUTPUT_OFF: do not print any output to
-   * console
-   *                        - CONSOLE_OUTPUT_ON: print before and after
-   * concentrations to console
-   *                        - CONSOLE_OUTPUT_VERBOSE: print all concentration
-   * matrices to console
    */
-  void setOutputConsole(CONSOLE_OUTPUT console_output) {
+  Diffusion(T *data, int rows, int cols)
-    if (console_output < CONSOLE_OUTPUT_OFF &&
+      : BaseSimulationGrid<T>(data, rows, cols) {
-        console_output > CONSOLE_OUTPUT_VERBOSE) {
+    init_alpha();
-      throw std::invalid_argument("Invalid console output option given!");
-    }
-
-    this->console_output = console_output;
   }
 
   /**
-   * @brief Set the Time Measure option. Off by default.
+   * @brief Construct a new 1D Diffusion object from a given data pointer and
+   * the length.
    *
-   * @param time_measure The following options are allowed:
-   *                     - TIME_MEASURE_OFF: Time of simulation is not printed
-   * to console
-   *                     - TIME_MEASURE_ON: Time of simulation run is printed to
-   * console
    */
-  void setTimeMeasure(TIME_MEASURE time_measure) {
+  Diffusion(T *data, std::size_t length) : BaseSimulationGrid<T>(data, length) {
-    if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) {
+    init_alpha();
-      throw std::invalid_argument("Invalid time measure option given!");
   }
 
-    this->time_measure = time_measure;
+  /**
+   * @brief Get the alphaX matrix.
+   *
+   * @return RowMajMat<T>& Reference to the alphaX matrix.
+   */
+  RowMajMat<T> &getAlphaX() { return alphaX; }
+
+  /**
+   * @brief Get the alphaY matrix.
+   *
+   * @return RowMajMat<T>& Reference to the alphaY matrix.
+   */
+  RowMajMat<T> &getAlphaY() {
+    tug_assert(
+        this->getDim(),
+        "Grid is not two dimensional, there is no domain in y-direction!");
+
+    return alphaY;
+  }
+
+  /**
+   * @brief Set the alphaX matrix.
+   *
+   * @param alphaX The new alphaX matrix.
+   */
+  void setAlphaX(const RowMajMat<T> &alphaX) { this->alphaX = alphaX; }
+
+  /**
+   * @brief Set the alphaY matrix.
+   *
+   * @param alphaY The new alphaY matrix.
+   */
+  void setAlphaY(const RowMajMat<T> &alphaY) {
+    tug_assert(
+        this->getDim(),
+        "Grid is not two dimensional, there is no domain in y-direction!");
+
+    this->alphaY = alphaY;
   }
 
   /**
@@ -177,33 +157,31 @@ public:
    *
    * @param timestep Valid timestep greater than zero.
    */
-  void setTimestep(T timestep) {
+  void setTimestep(T timestep) override {
-    if (timestep <= 0) {
+    tug_assert(timestep > 0, "Timestep has to be greater than zero.");
-      throw_invalid_argument("Timestep has to be greater than zero.");
-    }
 
     if constexpr (approach == FTCS_APPROACH ||
                   approach == CRANK_NICOLSON_APPROACH) {
       T cfl;
-      if (grid.getDim() == 1) {
+      if (this->getDim() == 1) {
 
-        const T deltaSquare = grid.getDelta();
+        const T deltaSquare = this->deltaCol();
-        const T maxAlpha = grid.getAlpha().maxCoeff();
+        const T maxAlpha = this->alphaX.maxCoeff();
 
         // Courant-Friedrichs-Lewy condition
         cfl = deltaSquare / (4 * maxAlpha);
-      } else if (grid.getDim() == 2) {
+      } else if (this->getDim() == 2) {
-        const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
+        const T deltaColSquare = this->deltaCol() * this->deltaCol();
         // will be 0 if 1D, else ...
-        const T deltaRowSquare = grid.getDeltaRow() * grid.getDeltaRow();
+        const T deltaRowSquare = this->deltaRow() * this->deltaRow();
         const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
 
         const T maxAlpha =
-            std::max(grid.getAlphaX().maxCoeff(), grid.getAlphaY().maxCoeff());
+            std::max(this->alphaX.maxCoeff(), this->alphaY.maxCoeff());
 
         cfl = minDeltaSquare / (4 * maxAlpha);
       }
-      const std::string dim = std::to_string(grid.getDim()) + "D";
+      const std::string dim = std::to_string(this->getDim()) + "D";
 
       const std::string &approachPrefix = this->approach_names[approach];
       std::cout << approachPrefix << "_" << dim << " :: CFL condition: " << cfl
@@ -244,20 +222,6 @@ public:
    */
   T getTimestep() const { return this->timestep; }
 
-  /**
-   * @brief Set the desired iterations to be calculated. A value greater
-   *        than zero must be specified here. Setting iterations is required.
-   *
-   * @param iterations Number of iterations to be simulated.
-   */
-  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 number of desired openMP Threads.
    *
@@ -277,19 +241,12 @@ public:
     }
   }
 
-  /**
-   * @brief Return the currently set iterations to be calculated.
-   *
-   * @return int Number of iterations.
-   */
-  int getIterations() const { return this->iterations; }
-
   /**
    * @brief Outputs the current concentrations of the grid on the console.
    *
    */
-  inline void printConcentrationsConsole() const {
+  void printConcentrationsConsole() const {
-    std::cout << grid.getConcentrations() << std::endl;
+    std::cout << this->getConcentrationMatrix() << std::endl;
     std::cout << std::endl;
   }
 
@@ -309,9 +266,9 @@ public:
 
     // string approachString = (approach == 0) ? "FTCS" : "BTCS";
     const std::string &approachString = this->approach_names[approach];
-    std::string row = std::to_string(grid.getRow());
+    std::string row = std::to_string(this->rows());
-    std::string col = std::to_string(grid.getCol());
+    std::string col = std::to_string(this->cols());
-    std::string numIterations = std::to_string(iterations);
+    std::string numIterations = std::to_string(this->getIterations());
 
     std::string filename =
         approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
@@ -330,7 +287,9 @@ public:
 
     // 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) {
+    if (this->getOutputCSV() == CSV_OUTPUT::XTREME) {
+      const auto &bc = this->getBoundaryConditions();
+
       Eigen::IOFormat one_row(Eigen::StreamPrecision, Eigen::DontAlignCols, "",
                               " ");
       file << bc.getBoundarySideValues(BC_SIDE_LEFT).format(one_row)
@@ -365,7 +324,7 @@ public:
     }
 
     Eigen::IOFormat do_not_align(Eigen::StreamPrecision, Eigen::DontAlignCols);
-    file << grid.getConcentrations().format(do_not_align) << std::endl;
+    file << this->getConcentrationMatrix().format(do_not_align) << std::endl;
     file << std::endl << std::endl;
     file.close();
   }
@@ -374,35 +333,42 @@ public:
    * @brief Method starts the simulation process with the previously set
    *        parameters.
    */
-  void run() {
+  void run() override {
-    if (this->timestep == -1) {
+    tug_assert(this->getTimestep() > 0, "Timestep is not set!");
-      throw_invalid_argument("Timestep is not set!");
+    tug_assert(this->getIterations() > 0, "Number of iterations are not set!");
-    }
-    if (this->iterations == -1) {
-      throw_invalid_argument("Number of iterations are not set!");
-    }
 
     std::string filename;
-    if (this->console_output > CONSOLE_OUTPUT_OFF) {
+    if (this->getOutputConsole() > CONSOLE_OUTPUT::OFF) {
       printConcentrationsConsole();
     }
-    if (this->csv_output > CSV_OUTPUT_OFF) {
+    if (this->getOutputCSV() > CSV_OUTPUT::OFF) {
       filename = createCSVfile();
     }
 
     auto begin = std::chrono::high_resolution_clock::now();
 
-    if constexpr (approach == FTCS_APPROACH) { // FTCS case
+    SimulationInput<T> sim_input = {.concentrations =
+                                        this->getConcentrationMatrix(),
+                                    .alphaX = this->getAlphaX(),
+                                    .alphaY = this->getAlphaY(),
+                                    .boundaries = this->getBoundaryConditions(),
+                                    .dim = this->getDim(),
+                                    .timestep = this->getTimestep(),
+                                    .rowMax = this->rows(),
+                                    .colMax = this->cols(),
+                                    .deltaRow = this->deltaRow(),
+                                    .deltaCol = this->deltaCol()};
 
-      for (int i = 0; i < iterations * innerIterations; i++) {
+    if constexpr (approach == FTCS_APPROACH) { // FTCS case
-        if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
+      for (int i = 0; i < this->getIterations() * innerIterations; i++) {
+        if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
           printConcentrationsConsole();
         }
-        if (csv_output >= CSV_OUTPUT_VERBOSE) {
+        if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
           printConcentrationsCSV(filename);
         }
 
-        FTCS(this->grid, this->bc, this->timestep, this->numThreads);
+        FTCS(sim_input, this->numThreads);
 
         // if (i % (iterations * innerIterations / 100) == 0) {
         //     double percentage = (double)i / ((double)iterations *
@@ -415,29 +381,28 @@ public:
     } else if constexpr (approach == BTCS_APPROACH) { // BTCS case
 
       if constexpr (solver == EIGEN_LU_SOLVER) {
-        for (int i = 0; i < iterations; i++) {
+        for (int i = 0; i < this->getIterations(); i++) {
-          if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
+          if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
             printConcentrationsConsole();
           }
-          if (csv_output >= CSV_OUTPUT_VERBOSE) {
+          if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
             printConcentrationsCSV(filename);
           }
 
-          BTCS_LU(this->grid, this->bc, this->timestep, this->numThreads);
+          BTCS_LU(sim_input, this->numThreads);
         }
       } else if constexpr (solver == THOMAS_ALGORITHM_SOLVER) {
-        for (int i = 0; i < iterations; i++) {
+        for (int i = 0; i < this->getIterations(); i++) {
-          if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
+          if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
             printConcentrationsConsole();
           }
-          if (csv_output >= CSV_OUTPUT_VERBOSE) {
+          if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
             printConcentrationsCSV(filename);
           }
 
-          BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
+          BTCS_Thomas(sim_input, this->numThreads);
         }
       }
-
     } else if constexpr (approach ==
                          CRANK_NICOLSON_APPROACH) { // Crank-Nicolson case
 
@@ -449,22 +414,22 @@ public:
       RowMajMat<T> concentrations;
       RowMajMat<T> concentrationsFTCS;
       RowMajMat<T> concentrationsResult;
-      for (int i = 0; i < iterations * innerIterations; i++) {
+      for (int i = 0; i < this->getIterations() * innerIterations; i++) {
-        if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
+        if (this->getOutputConsole() == CONSOLE_OUTPUT::VERBOSE && i > 0) {
           printConcentrationsConsole();
         }
-        if (csv_output >= CSV_OUTPUT_VERBOSE) {
+        if (this->getOutputCSV() >= CSV_OUTPUT::VERBOSE) {
           printConcentrationsCSV(filename);
         }
 
-        concentrations = grid.getConcentrations();
+        concentrations = this->getConcentrationMatrix();
         FTCS(this->grid, this->bc, this->timestep, this->numThreads);
-        concentrationsFTCS = grid.getConcentrations();
+        concentrationsFTCS = this->getConcentrationMatrix();
-        grid.setConcentrations(concentrations);
+        this->getConcentrationMatrix() = concentrations;
-        BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
+        BTCS_Thomas(sim_input, this->numThreads);
-        concentrationsResult =
+        concentrationsResult = beta * concentrationsFTCS +
-            beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations();
+                               (1 - beta) * this->getConcentrationMatrix();
-        grid.setConcentrations(concentrationsResult);
+        this->getConcentrationMatrix() = concentrationsResult;
       }
     }
 
@@ -472,33 +437,18 @@ public:
     auto milliseconds =
         std::chrono::duration_cast<std::chrono::milliseconds>(end - begin);
 
-    if (this->console_output > CONSOLE_OUTPUT_OFF) {
+    if (this->getOutputConsole() > CONSOLE_OUTPUT::OFF) {
       printConcentrationsConsole();
     }
-    if (this->csv_output > CSV_OUTPUT_OFF) {
+    if (this->getOutputCSV() > CSV_OUTPUT::OFF) {
       printConcentrationsCSV(filename);
     }
-    if (this->time_measure > TIME_MEASURE_OFF) {
+    if (this->getTimeMeasure() > TIME_MEASURE::OFF) {
       const std::string &approachString = this->approach_names[approach];
-      const std::string dimString = std::to_string(grid.getDim()) + "D";
+      const std::string dimString = std::to_string(this->getDim()) + "D";
       std::cout << approachString << dimString << ":: run() finished in "
                 << milliseconds.count() << "ms" << std::endl;
     }
   }
-
-private:
-  T timestep{-1};
-  int iterations{-1};
-  int innerIterations{1};
-  int numThreads{omp_get_num_procs()};
-  CSV_OUTPUT csv_output{CSV_OUTPUT_OFF};
-  CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF};
-  TIME_MEASURE time_measure{TIME_MEASURE_OFF};
-
-  Grid<T> &grid;
-  Boundary<T> &bc;
-
-  const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
 };
 } // namespace tug
-#endif // SIMULATION_H_

include/tug/Grid.hpp

@@ -1,392 +0,0 @@
-#ifndef GRID_H_
-#define GRID_H_
-
-/**
- * @file Grid.hpp
- * @brief API of Grid class, that holds a matrix with concenctrations and a
- *        respective matrix/matrices of alpha coefficients.
- *
- */
-
-#include "Core/Matrix.hpp"
-#include <Eigen/Core>
-#include <Eigen/Sparse>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
-#include <stdexcept>
-
-namespace tug {
-
-/**
- * @brief Holds a matrix with concenctration and respective matrix/matrices of
- * alpha coefficients.
- *
- * @tparam T Type to be used for matrices, e.g. double or float
- */
-template <class T> class Grid {
-public:
-  /**
-   * @brief Constructs a new 1D-Grid object of a given length, which holds a
-   * matrix with concentrations and a respective matrix of alpha coefficients.
-   *        The domain length is per default the same as the length. The
-   * concentrations are all 20 by default and the alpha coefficients are 1.
-   *
-   * @param length Length of the 1D-Grid. Must be greater than 3.
-   */
-  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 =
-        static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
-
-    this->concentrations = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
-    this->alphaX = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
-  }
-
-  /**
-   * @brief Constructs a new 2D-Grid object of given dimensions, which holds a
-   * matrix with concentrations and the respective matrices of alpha coefficient
-   * for each direction. The domain in x- and y-direction is per default equal
-   * to the col length and row length, respectively. The concentrations are all
-   * 20 by default across the entire grid and the alpha coefficients 1 in both
-   * directions.
-   *
-   * @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)
-      : row(_row), col(_col), domainRow(_row), domainCol(_col) {
-    if (row <= 1 || col <= 1) {
-      throw std::invalid_argument(
-          "At least one dimension is 1. Use 1D grid for better results.");
-    }
-
-    this->dim = 2;
-    this->deltaRow =
-        static_cast<T>(this->domainRow) / static_cast<T>(this->row); // -> 1
-    this->deltaCol =
-        static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
-
-    this->concentrations = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
-    this->alphaX = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
-    this->alphaY = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
-  }
-
-  /**
-   * @brief Sets the concentrations matrix for a 1D or 2D-Grid.
-   *
-   * @param concentrations An Eigen3 MatrixX<T> holding the concentrations.
-   * Matrix must have correct dimensions as defined in row and col. (Or length,
-   * in 1D case).
-   */
-  void setConcentrations(const RowMajMat<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 Sets the concentrations matrix for a 1D or 2D-Grid.
-   *
-   * @param concentrations A pointer to an array holding the concentrations.
-   * Array must have correct dimensions as defined in row and col. (Or length,
-   * in 1D case). There is no check for correct dimensions, so be careful!
-   */
-  void setConcentrations(T *concentrations) {
-    tug::RowMajMatMap<T> map(concentrations, this->row, this->col);
-    this->concentrations = map;
-  }
-
-  /**
-   * @brief Gets the concentrations matrix for a Grid.
-   *
-   * @return An Eigen3 matrix holding the concentrations and having
-   * the same dimensions as the grid.
-   */
-  auto &getConcentrations() { return this->concentrations; }
-
-  /**
-   * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
-   * dimensional.
-   *
-   * @param alpha An Eigen3 MatrixX<T> with 1 row holding the alpha
-   * coefficients. Matrix columns must have same size as length of grid.
-   */
-  void setAlpha(const RowMajMat<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 1D-Grid. Grid must be one
-   * dimensional.
-   *
-   * @param alpha A pointer to an array holding the alpha coefficients. Array
-   * must have correct dimensions as defined in length. There is no check for
-   * correct dimensions, so be careful!
-   */
-  void setAlpha(T *alpha) {
-    if (dim != 1) {
-      throw std::invalid_argument(
-          "Grid is not one dimensional, you should probably "
-          "use 2D setter function!");
-    }
-    RowMajMatMap<T> map(alpha, 1, this->col);
-    this->alphaX = map;
-  }
-
-  /**
-   * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
-   * dimensional.
-   *
-   * @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 MatrixX<T> holding the alpha coefficients in
-   * y-direction. Matrix must be of same size as the grid.
-   */
-  void setAlpha(const RowMajMat<T> &alphaX, const RowMajMat<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 Set the alpha coefficients of a 2D-Grid. Grid must be two
-   * dimensional.
-   *
-   * @param alphaX A pointer to an array holding the alpha coefficients in
-   * x-direction. Array must have correct dimensions as defined in row and col.
-   * There is no check for correct dimensions, so be careful!
-   * @param alphaY A pointer to an array holding the alpha coefficients in
-   * y-direction. Array must have correct dimensions as defined in row and col.
-   * There is no check for correct dimensions, so be careful!
-   */
-  void setAlpha(T *alphaX, T *alphaY) {
-    if (dim != 2) {
-      throw std::invalid_argument(
-          "Grid is not two dimensional, you should probably "
-          "use 1D setter function!");
-    }
-    RowMajMatMap<T> mapX(alphaX, this->row, this->col);
-    RowMajMatMap<T> mapY(alphaY, this->row, this->col);
-    this->alphaX = mapX;
-    this->alphaY = mapY;
-  }
-
-  /**
-   * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
-   * dimensional.
-   *
-   * @return A matrix with 1 row holding the alpha coefficients.
-   */
-  const auto &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 A matrix holding the alpha coefficients in x-direction.
-   */
-  const auto &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 A matrix holding the alpha coefficients in y-direction.
-   */
-  const auto &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 Dimensions, either 1 or 2.
-   */
-  int getDim() const { return this->dim; }
-
-  /**
-   * @brief Gets length of 1D grid. Must be one dimensional grid.
-   *
-   * @return Length of 1D grid.
-   */
-  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 Number of rows.
-   */
-  int getRow() const { return this->row; }
-
-  /**
-   * @brief Gets the number of columns of the grid.
-   *
-   * @return Number of columns.
-   */
-  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) {
-    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 positive.
-   * @param domainCol A double value of the domain size in x-direction. Must
-   * be positive.
-   */
-  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 Delta value.
-   */
-  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 Delta value in x-direction.
-   */
-  T getDeltaCol() const { return this->deltaCol; }
-
-  /**
-   * @brief Gets the delta value in y-direction. Must be two dimensional grid.
-   *
-   * @return Delta value in y-direction.
-   */
-  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 row{1};     // number of grid rows
-  int dim;        // 1D or 2D
-  T domainCol;    // number of domain columns
-  T domainRow{0}; // number of domain rows
-  T deltaCol;     // delta in x-direction (between columns)
-  T deltaRow{0};  // delta in y-direction (between rows)
-
-  RowMajMat<T> concentrations; // Matrix holding grid concentrations
-  RowMajMat<T> alphaX; // Matrix holding alpha coefficients in x-direction
-  RowMajMat<T> alphaY; // Matrix holding alpha coefficients in y-direction
-
-  static constexpr T MAT_INIT_VAL = 0;
-};
-
-using Grid64 = Grid<double>;
-using Grid32 = Grid<float>;
-} // namespace tug
-#endif // GRID_H_

naaice/BTCS_2D_NAAICE.cpp

@@ -5,8 +5,7 @@
 #include <ostream>
 #include <stdexcept>
 #include <string>
-#include <string_view>
+#include <tug/Diffusion.hpp>
-#include <tug/Simulation.hpp>
 #include <vector>
 
 #include "files.hpp"
@@ -105,15 +104,14 @@ int main(int argc, char *argv[]) {
   // create a grid with a 5 x 10 field
   constexpr int row = 5;
   constexpr int col = 10;
-  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);
   Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col);
-  grid.setConcentrations(concentrations);
+  Grid64 grid(concentrations);
 
+  grid.setDomain(0.005, 0.01);
   const double sum_init = concentrations.sum();
 
   // // (optional) set alphas of the grid, e.g.:
@@ -141,7 +139,7 @@ int main(int argc, char *argv[]) {
   // // ************************
 
   // set up a simulation environment
-  Simulation simulation = Simulation(grid, bc); // grid,boundary
+  Diffusion simulation(grid, bc); // grid,boundary
 
   // set the timestep of the simulation
   simulation.setTimestep(360); // timestep

naaice/NAAICE_dble_vs_float.cpp

@@ -7,11 +7,11 @@
 #include <stdexcept>
 #include <string>
 #include <string_view>
-#include <tug/Simulation.hpp>
 #include <type_traits>
 #include <vector>
 
-#include "files.hpp"
+#include <files.hpp>
+#include <tug/Diffusion.hpp>
 
 using namespace tug;
 
@@ -114,15 +114,14 @@ template <class T, tug::APPROACH app> int doWork(int ngrid) {
 
   std::cout << name << " grid: " << ngrid << std::endl;
 
-  Grid<T> grid(ngrid, ngrid);
   // Grid64 grid(ngrid, ngrid);
 
   // (optional) set the domain, e.g.:
-  grid.setDomain(0.1, 0.1);
 
   Eigen::MatrixX<T> initConc64 = Eigen::MatrixX<T>::Constant(ngrid, ngrid, 0);
   initConc64(50, 50) = 1E-6;
-  grid.setConcentrations(initConc64);
+  Grid<T> grid(initConc64);
+  grid.setDomain(0.1, 0.1);
 
   const T sum_init64 = initConc64.sum();
 
@@ -142,8 +141,7 @@ template <class T, tug::APPROACH app> int doWork(int ngrid) {
   bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
 
   // set up a simulation environment
-  Simulation Sim =
+  Diffusion Sim(grid, bc); // grid_64,boundary,simulation-approach
-      Simulation<T, app>(grid, bc); // grid_64,boundary,simulation-approach
 
   // Sim64.setSolver(THOMAS_ALGORITHM_SOLVER);
 

test/CMakeLists.txt

@@ -10,12 +10,12 @@ FetchContent_MakeAvailable(googletest)
 
 
 add_executable(testTug 
-testSimulation.cpp 
+setup.cpp
-testGrid.cpp 
+testDiffusion.cpp 
 testFTCS.cpp 
 testBoundary.cpp
 )
-target_link_libraries(testTug tug GTest::gtest_main)
+target_link_libraries(testTug tug GTest::gtest)
   
 include(GoogleTest)
 gtest_discover_tests(testTug)

test/TestUtils.hpp

@@ -1,3 +1,4 @@
+#include "tug/Core/Matrix.hpp"
 #include <Eigen/Core>
 #include <Eigen/Dense>
 #include <fstream>
@@ -8,7 +9,7 @@
 
 #define TUG_TEST(x) TEST(Tug, x)
 
-inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) {
+inline tug::RowMajMat<double> CSV2Eigen(std::string file2Convert) {
 
   std::vector<double> matrixEntries;
 
@@ -31,21 +32,20 @@ inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) {
     }
   }
 
-  return Eigen::Map<
+  return tug::RowMajMatMap<double>(matrixEntries.data(), matrixRowNumber,
-      Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
-      matrixEntries.data(), matrixRowNumber,
                                    matrixEntries.size() / matrixRowNumber);
 }
 
-inline bool checkSimilarity(Eigen::MatrixXd a, Eigen::MatrixXd b,
+inline bool checkSimilarity(tug::RowMajMat<double> &a,
+                            tug::RowMajMatMap<double> &b,
                             double precision = 1e-5) {
   return a.isApprox(b, precision);
 }
 
-inline bool checkSimilarityV2(Eigen::MatrixXd a, Eigen::MatrixXd b,
+inline bool checkSimilarityV2(tug::RowMajMat<double> &a,
-                              double maxDiff) {
+                              tug::RowMajMatMap<double> &b, double maxDiff) {
 
-  Eigen::MatrixXd diff = a - b;
+  tug::RowMajMat<double> diff = a - b;
   double maxCoeff = diff.maxCoeff();
   return abs(maxCoeff) < maxDiff;
 }

test/setup.cpp

@@ -1,2 +1,7 @@
-#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
+#include <gtest/gtest.h>
-#include <doctest/doctest.h>
+
+int main(int argc, char **argv) {
+  ::testing::InitGoogleTest(&argc, argv);
+  GTEST_FLAG_SET(death_test_style, "threadsafe");
+  return RUN_ALL_TESTS();
+}

test/testBoundary.cpp

@@ -1,3 +1,5 @@
+#include "gtest/gtest.h"
+#include <stdexcept>
 #include <tug/Boundary.hpp>
 #include <utility>
 #include <vector>
@@ -15,7 +17,7 @@ BOUNDARY_TEST(Element) {
   EXPECT_NO_THROW(BoundaryElement<double>());
   EXPECT_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED);
   EXPECT_DOUBLE_EQ(boundaryElementClosed.getValue(), -1);
-  EXPECT_ANY_THROW(boundaryElementClosed.setValue(0.2));
+  EXPECT_THROW(boundaryElementClosed.setValue(0.2), std::invalid_argument);
 
   BoundaryElement boundaryElementConstant = BoundaryElement(0.1);
   EXPECT_NO_THROW(BoundaryElement(0.1));
@@ -26,10 +28,8 @@ BOUNDARY_TEST(Element) {
 }
 
 BOUNDARY_TEST(Class) {
-  Grid grid1D = Grid64(10);
+  Boundary<double> boundary1D(10);
-  Grid grid2D = Grid64(10, 12);
+  Boundary<double> boundary2D(10, 12);
-  Boundary boundary1D = Boundary(grid1D);
-  Boundary boundary2D = Boundary(grid2D);
   vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));
 
   constexpr double inner_condition_value = -5;
@@ -43,31 +43,34 @@ BOUNDARY_TEST(Class) {
   col_ibc[0] = innerBoundary;
 
   {
-    EXPECT_NO_THROW(Boundary boundary(grid1D));
     EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1);
     EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1);
     EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
               BC_TYPE_CLOSED);
-    EXPECT_ANY_THROW(boundary1D.getBoundarySide(BC_SIDE_TOP));
+    EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_TOP),
-    EXPECT_ANY_THROW(boundary1D.getBoundarySide(BC_SIDE_BOTTOM));
+                 ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*");
+    EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_BOTTOM),
+                 ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*");
     EXPECT_NO_THROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT));
-    EXPECT_ANY_THROW(boundary1D.setBoundarySideClosed(BC_SIDE_TOP));
+    EXPECT_DEATH(boundary1D.setBoundarySideClosed(BC_SIDE_TOP),
+                 ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*");
     EXPECT_NO_THROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
     EXPECT_DOUBLE_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
-    EXPECT_ANY_THROW(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2));
+    EXPECT_DEATH(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2),
+                 ".*Index is selected either too large or too small.*");
     EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
               BC_TYPE_CONSTANT);
     EXPECT_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
               boundary1DVector[0].getType());
 
     EXPECT_NO_THROW(boundary1D.setInnerBoundary(0, inner_condition_value));
-    EXPECT_ANY_THROW(boundary1D.setInnerBoundary(0, 0, inner_condition_value));
+    EXPECT_DEATH(boundary1D.setInnerBoundary(0, 0, inner_condition_value),
+                 ".*only available for 2D grids.*");
     EXPECT_EQ(boundary1D.getInnerBoundary(0), innerBoundary);
     EXPECT_FALSE(boundary1D.getInnerBoundary(1).first);
   }
 
   {
-    EXPECT_NO_THROW(Boundary boundary(grid1D));
     EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10);
     EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10);
     EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12);
@@ -80,13 +83,15 @@ BOUNDARY_TEST(Class) {
     EXPECT_NO_THROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP));
     EXPECT_NO_THROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
     EXPECT_DOUBLE_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
-    EXPECT_ANY_THROW(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12));
+    EXPECT_DEATH(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12),
+                 ".*too large or too small.*");
     EXPECT_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
               BC_TYPE_CONSTANT);
     EXPECT_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
               boundary1DVector[0].getType());
 
-    EXPECT_ANY_THROW(boundary2D.setInnerBoundary(0, inner_condition_value));
+    EXPECT_DEATH(boundary2D.setInnerBoundary(0, inner_condition_value),
+                 ".* 1D .*");
     EXPECT_NO_THROW(boundary2D.setInnerBoundary(0, 1, inner_condition_value));
     EXPECT_EQ(boundary2D.getInnerBoundary(0, 1), innerBoundary);
     EXPECT_FALSE(boundary2D.getInnerBoundary(0, 2).first);

test/testDiffusion.cpp

@@ -0,0 +1,216 @@
+#include "TestUtils.hpp"
+#include "tug/Core/Matrix.hpp"
+#include "gtest/gtest.h"
+#include <gtest/gtest.h>
+#include <stdexcept>
+#include <tug/Diffusion.hpp>
+
+#include <Eigen/src/Core/Matrix.h>
+#include <string>
+
+// include the configured header file
+#include <testSimulation.hpp>
+
+#define DIFFUSION_TEST(x) TEST(Diffusion, x)
+
+using namespace Eigen;
+using namespace std;
+using namespace tug;
+
+constexpr int row = 11;
+constexpr int col = 11;
+
+template <tug::APPROACH approach>
+Diffusion<double, approach> setupSimulation(RowMajMat<double> &concentrations,
+                                            double timestep, int iterations) {
+  int domain_row = 10;
+  int domain_col = 10;
+
+  // Grid
+  // RowMajMat<double> concentrations = MatrixXd::Constant(row, col, 0);
+  concentrations(5, 5) = 1;
+
+  Diffusion<double, approach> diffusiongrid(concentrations);
+
+  diffusiongrid.getConcentrationMatrix() = concentrations;
+  diffusiongrid.setDomain(domain_row, domain_col);
+
+  diffusiongrid.setTimestep(timestep);
+  diffusiongrid.setIterations(iterations);
+  diffusiongrid.setDomain(domain_row, domain_col);
+
+  MatrixXd alpha = MatrixXd::Constant(row, col, 1);
+  for (int i = 0; i < 5; i++) {
+    for (int j = 0; j < 6; j++) {
+      alpha(i, j) = 0.01;
+    }
+  }
+  for (int i = 0; i < 5; i++) {
+    for (int j = 6; j < 11; j++) {
+      alpha(i, j) = 0.001;
+    }
+  }
+  for (int i = 5; i < 11; i++) {
+    for (int j = 6; j < 11; j++) {
+      alpha(i, j) = 0.1;
+    }
+  }
+  diffusiongrid.setAlphaX(alpha);
+  diffusiongrid.setAlphaY(alpha);
+
+  return diffusiongrid;
+}
+
+constexpr double timestep = 0.001;
+constexpr double iterations = 7000;
+
+DIFFUSION_TEST(EqualityFTCS) {
+  // set string from the header file
+  string test_path = testSimulationCSVDir;
+  RowMajMat<double> reference = CSV2Eigen(test_path);
+  cout << "FTCS Test: " << endl;
+
+  RowMajMat<double> concentrations = MatrixXd::Constant(row, col, 0);
+
+  Diffusion<double, tug::FTCS_APPROACH> sim =
+      setupSimulation<tug::FTCS_APPROACH>(concentrations, timestep, iterations);
+
+  // Boundary bc = Boundary(grid);
+
+  // Simulation
+
+  // Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
+  // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
+  // sim.setTimestep(timestep);
+  // sim.setIterations(iterations);
+  sim.run();
+
+  cout << endl;
+  EXPECT_TRUE(checkSimilarity(reference, sim.getConcentrationMatrix(), 0.1));
+}
+
+DIFFUSION_TEST(EqualityBTCS) {
+  // set string from the header file
+  string test_path = testSimulationCSVDir;
+  RowMajMat<double> reference = CSV2Eigen(test_path);
+  cout << "BTCS Test: " << endl;
+
+  RowMajMat<double> concentrations = MatrixXd::Constant(row, col, 0);
+
+  Diffusion<double, tug::BTCS_APPROACH> sim =
+      setupSimulation<tug::BTCS_APPROACH>(concentrations, timestep,
+                                          iterations); // Boundary
+
+  // Boundary bc = Boundary(grid);
+
+  // Simulation
+  // Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
+  // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
+  // sim.setTimestep(timestep);
+  // sim.setIterations(iterations);
+  sim.run();
+
+  cout << endl;
+  EXPECT_TRUE(checkSimilarityV2(reference, sim.getConcentrationMatrix(), 0.01));
+}
+
+DIFFUSION_TEST(InitializeEnvironment) {
+  int rc = 12;
+  RowMajMat<double> concentrations(rc, rc);
+  // Grid64 grid(concentrations);
+  // Boundary boundary(grid);
+
+  EXPECT_NO_FATAL_FAILURE(Diffusion<double> sim(concentrations));
+}
+
+// DIFFUSION_TEST(SimulationEnvironment) {
+//   int rc = 12;
+//   Eigen::MatrixXd concentrations(rc, rc);
+//   Grid64 grid(concentrations);
+//   grid.initAlpha();
+//   Boundary boundary(grid);
+//   Diffusion<double, tug::FTCS_APPROACH> sim(grid, boundary);
+
+//   EXPECT_EQ(sim.getIterations(), 1);
+
+//   EXPECT_NO_THROW(sim.setIterations(2000));
+//   EXPECT_EQ(sim.getIterations(), 2000);
+//   EXPECT_THROW(sim.setIterations(-300), std::invalid_argument);
+
+//   EXPECT_NO_THROW(sim.setTimestep(0.1));
+//   EXPECT_DOUBLE_EQ(sim.getTimestep(), 0.1);
+//   EXPECT_DEATH(sim.setTimestep(-0.3), ".* greater than zero.*");
+// }
+
+DIFFUSION_TEST(ClosedBoundaries) {
+  constexpr std::uint32_t nrows = 5;
+  constexpr std::uint32_t ncols = 5;
+
+  RowMajMat<double> concentrations =
+      RowMajMat<double>::Constant(nrows, ncols, 1.0);
+  RowMajMat<double> alphax = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
+  RowMajMat<double> alphay = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
+
+  Diffusion<double> sim(concentrations);
+  sim.getAlphaX() = alphax;
+  sim.getAlphaY() = alphay;
+
+  // tug::Grid64 grid(concentrations);
+
+  // grid.setAlpha(alphax, alphay);
+
+  // tug::Boundary bc(grid);
+  auto &bc = sim.getBoundaryConditions();
+  bc.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1.0);
+  bc.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1.0);
+  bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0);
+  bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0);
+
+  // tug::Diffusion<double> sim(grid, bc);
+  sim.setTimestep(1);
+  sim.setIterations(1);
+
+  RowMajMat<double> input_values(concentrations);
+  sim.run();
+
+  EXPECT_TRUE(
+      checkSimilarityV2(input_values, sim.getConcentrationMatrix(), 1E-12));
+}
+
+DIFFUSION_TEST(ConstantInnerCell) {
+  constexpr std::uint32_t nrows = 5;
+  constexpr std::uint32_t ncols = 5;
+
+  RowMajMat<double> concentrations =
+      RowMajMat<double>::Constant(nrows, ncols, 1.0);
+  RowMajMat<double> alphax = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
+  RowMajMat<double> alphay = RowMajMat<double>::Constant(nrows, ncols, 1E-5);
+
+  Diffusion<double> sim(concentrations);
+  sim.getAlphaX() = alphax;
+  sim.getAlphaY() = alphay;
+
+  // tug::Grid64 grid(concentrations);
+  // grid.setAlpha(alphax, alphay);
+
+  // tug::Boundary bc(grid);
+  auto &bc = sim.getBoundaryConditions();
+  // inner
+  bc.setInnerBoundary(2, 2, 0);
+
+  // tug::Diffusion<double> sim(grid, bc);
+  sim.setTimestep(1);
+  sim.setIterations(1);
+
+  MatrixXd input_values(concentrations);
+  sim.run();
+
+  const auto &concentrations_result = sim.getConcentrationMatrix();
+
+  EXPECT_DOUBLE_EQ(concentrations_result(2, 2), 0);
+  EXPECT_LT(concentrations_result.sum(), input_values.sum());
+
+  EXPECT_FALSE((concentrations_result.array() > 1.0).any());
+
+  EXPECT_FALSE((concentrations_result.array() < 0.0).any());
+}

test/testGrid.cpp

@@ -1,262 +0,0 @@
-#include <Eigen/Core>
-#include <tug/Grid.hpp>
-#include <vector>
-
-#include <gtest/gtest.h>
-
-using namespace Eigen;
-using namespace std;
-using namespace tug;
-
-#define GRID_TEST(x) TEST(Grid, x)
-
-GRID_TEST(InvalidConstructorParams) {
-  EXPECT_ANY_THROW(Grid64(2));
-  EXPECT_ANY_THROW(Grid64(1, 4));
-  EXPECT_ANY_THROW(Grid64(4, 1));
-}
-
-GRID_TEST(Grid64OneDimensional) {
-  int l = 12;
-  Grid64 grid(l);
-
-  {
-    EXPECT_EQ(grid.getDim(), 1);
-    EXPECT_EQ(grid.getLength(), l);
-    EXPECT_EQ(grid.getCol(), l);
-    EXPECT_EQ(grid.getRow(), 1);
-
-    EXPECT_EQ(grid.getConcentrations().rows(), 1);
-    EXPECT_EQ(grid.getConcentrations().cols(), l);
-    EXPECT_EQ(grid.getAlpha().rows(), 1);
-    EXPECT_EQ(grid.getAlpha().cols(), l);
-    EXPECT_EQ(grid.getDeltaCol(), 1);
-
-    EXPECT_ANY_THROW(grid.getAlphaX());
-    EXPECT_ANY_THROW(grid.getAlphaY());
-    EXPECT_ANY_THROW(grid.getDeltaRow());
-  }
-
-  {
-    // correct concentrations matrix
-    MatrixXd concentrations = MatrixXd::Constant(1, l, 3);
-    EXPECT_NO_THROW(grid.setConcentrations(concentrations));
-
-    // false concentrations matrix
-    MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4);
-    EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations));
-  }
-
-  {
-    // correct alpha matrix
-    MatrixXd alpha = MatrixXd::Constant(1, l, 3);
-    EXPECT_NO_THROW(grid.setAlpha(alpha));
-
-    EXPECT_ANY_THROW(grid.setAlpha(alpha, alpha));
-
-    grid.setAlpha(alpha);
-    EXPECT_EQ(grid.getAlpha(), alpha);
-    EXPECT_ANY_THROW(grid.getAlphaX());
-    EXPECT_ANY_THROW(grid.getAlphaY());
-
-    // false alpha matrix
-    MatrixXd wAlpha = MatrixXd::Constant(3, l, 2);
-    EXPECT_ANY_THROW(grid.setAlpha(wAlpha));
-  }
-
-  {
-    int d = 8;
-    // set 1D domain
-    EXPECT_NO_THROW(grid.setDomain(d));
-
-    // set 2D domain
-    EXPECT_ANY_THROW(grid.setDomain(d, d));
-
-    grid.setDomain(d);
-    EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(d) / double(l));
-    EXPECT_ANY_THROW(grid.getDeltaRow());
-
-    // set too small domain
-    d = 0;
-    EXPECT_ANY_THROW(grid.setDomain(d));
-    d = -2;
-    EXPECT_ANY_THROW(grid.setDomain(d));
-  }
-}
-
-GRID_TEST(Grid64Quadratic) {
-  int rc = 12;
-  Grid64 grid(rc, rc);
-
-  {
-    EXPECT_EQ(grid.getDim(), 2);
-    EXPECT_ANY_THROW(grid.getLength());
-    EXPECT_EQ(grid.getCol(), rc);
-    EXPECT_EQ(grid.getRow(), rc);
-
-    EXPECT_EQ(grid.getConcentrations().rows(), rc);
-    EXPECT_EQ(grid.getConcentrations().cols(), rc);
-    EXPECT_ANY_THROW(grid.getAlpha());
-
-    EXPECT_EQ(grid.getAlphaX().rows(), rc);
-    EXPECT_EQ(grid.getAlphaX().cols(), rc);
-    EXPECT_EQ(grid.getAlphaY().rows(), rc);
-    EXPECT_EQ(grid.getAlphaY().cols(), rc);
-    EXPECT_EQ(grid.getDeltaRow(), 1);
-    EXPECT_EQ(grid.getDeltaCol(), 1);
-  }
-
-  {
-    // correct concentrations matrix
-    MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2);
-    EXPECT_NO_THROW(grid.setConcentrations(concentrations));
-
-    // false concentrations matrix
-    MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1);
-    EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(rc + 3, rc, 4);
-    EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6);
-    EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations));
-  }
-
-  {
-    // correct alpha matrices
-    MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
-    MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
-    EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
-
-    EXPECT_ANY_THROW(grid.setAlpha(alphax));
-
-    grid.setAlpha(alphax, alphay);
-    EXPECT_EQ(grid.getAlphaX(), alphax);
-    EXPECT_EQ(grid.getAlphaY(), alphay);
-    EXPECT_ANY_THROW(grid.getAlpha());
-
-    // false alpha matrices
-    alphax = MatrixXd::Constant(rc + 3, rc + 1, 3);
-    EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay));
-    alphay = MatrixXd::Constant(rc + 2, rc + 1, 3);
-    EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay));
-  }
-
-  {
-    int dr = 8;
-    int dc = 9;
-
-    // set 1D domain
-    EXPECT_ANY_THROW(grid.setDomain(dr));
-
-    // set 2D domain
-    EXPECT_NO_THROW(grid.setDomain(dr, dc));
-
-    grid.setDomain(dr, dc);
-    EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(rc));
-    EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(rc));
-
-    // set too small domain
-    dr = 0;
-    EXPECT_ANY_THROW(grid.setDomain(dr, dc));
-    dr = 8;
-    dc = 0;
-    EXPECT_ANY_THROW(grid.setDomain(dr, dc));
-    dr = -2;
-    EXPECT_ANY_THROW(grid.setDomain(dr, dc));
-  }
-}
-
-GRID_TEST(Grid64NonQuadratic) {
-  int r = 12;
-  int c = 15;
-  Grid64 grid(r, c);
-
-  {
-    EXPECT_EQ(grid.getDim(), 2);
-    EXPECT_ANY_THROW(grid.getLength());
-    EXPECT_EQ(grid.getCol(), c);
-    EXPECT_EQ(grid.getRow(), r);
-
-    EXPECT_EQ(grid.getConcentrations().rows(), r);
-    EXPECT_EQ(grid.getConcentrations().cols(), c);
-    EXPECT_ANY_THROW(grid.getAlpha());
-
-    EXPECT_EQ(grid.getAlphaX().rows(), r);
-    EXPECT_EQ(grid.getAlphaX().cols(), c);
-    EXPECT_EQ(grid.getAlphaY().rows(), r);
-    EXPECT_EQ(grid.getAlphaY().cols(), c);
-    EXPECT_EQ(grid.getDeltaRow(), 1);
-    EXPECT_EQ(grid.getDeltaCol(), 1);
-  }
-
-  {
-    // correct concentrations matrix
-    MatrixXd concentrations = MatrixXd::Constant(r, c, 2);
-    EXPECT_NO_THROW(grid.setConcentrations(concentrations));
-
-    // false concentrations matrix
-    MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6);
-    EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(r + 3, c, 3);
-    EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2);
-    EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations));
-  }
-
-  {
-    // correct alpha matrices
-    MatrixXd alphax = MatrixXd::Constant(r, c, 2);
-    MatrixXd alphay = MatrixXd::Constant(r, c, 4);
-    EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
-
-    EXPECT_ANY_THROW(grid.setAlpha(alphax));
-
-    grid.setAlpha(alphax, alphay);
-    EXPECT_EQ(grid.getAlphaX(), alphax);
-    EXPECT_EQ(grid.getAlphaY(), alphay);
-    EXPECT_ANY_THROW(grid.getAlpha());
-
-    // false alpha matrices
-    alphax = MatrixXd::Constant(r + 3, c + 1, 3);
-    EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay));
-    alphay = MatrixXd::Constant(r + 2, c + 1, 5);
-    EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay));
-  }
-
-  {
-    int dr = 8;
-    int dc = 9;
-
-    // set 1D domain
-    EXPECT_ANY_THROW(grid.setDomain(dr));
-
-    // set 2D domain
-    EXPECT_NO_THROW(grid.setDomain(dr, dc));
-
-    grid.setDomain(dr, dc);
-    EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(c));
-    EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(r));
-
-    // set too small domain
-    dr = 0;
-    EXPECT_ANY_THROW(grid.setDomain(dr, dc));
-    dr = 8;
-    dc = -1;
-    EXPECT_ANY_THROW(grid.setDomain(dr, dc));
-    dr = -2;
-    EXPECT_ANY_THROW(grid.setDomain(dr, dc));
-  }
-
-  {
-    std::vector<double> concentrations(r * c);
-
-    for (int i = 0; i < r * c; i++) {
-      concentrations[i] = i;
-    }
-
-    grid.setConcentrations(concentrations.data());
-
-    EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0);
-    EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1);
-    EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c);
-  }
-}

test/testSimulation.cpp

@@ -1,182 +0,0 @@
-#include "TestUtils.hpp"
-#include <gtest/gtest.h>
-#include <stdexcept>
-#include <tug/Simulation.hpp>
-
-#include <Eigen/src/Core/Matrix.h>
-#include <string>
-
-// include the configured header file
-#include <testSimulation.hpp>
-
-#define DIFFUSION_TEST(x) TEST(Diffusion, x)
-
-using namespace Eigen;
-using namespace std;
-using namespace tug;
-
-Grid64 setupSimulation(double timestep, int iterations) {
-  int row = 11;
-  int col = 11;
-  int domain_row = 10;
-  int domain_col = 10;
-
-  // Grid
-  Grid grid = Grid64(row, col);
-  grid.setDomain(domain_row, domain_col);
-
-  MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
-  concentrations(5, 5) = 1;
-  grid.setConcentrations(concentrations);
-
-  MatrixXd alpha = MatrixXd::Constant(row, col, 1);
-  for (int i = 0; i < 5; i++) {
-    for (int j = 0; j < 6; j++) {
-      alpha(i, j) = 0.01;
-    }
-  }
-  for (int i = 0; i < 5; i++) {
-    for (int j = 6; j < 11; j++) {
-      alpha(i, j) = 0.001;
-    }
-  }
-  for (int i = 5; i < 11; i++) {
-    for (int j = 6; j < 11; j++) {
-      alpha(i, j) = 0.1;
-    }
-  }
-  grid.setAlpha(alpha, alpha);
-
-  return grid;
-}
-
-constexpr double timestep = 0.001;
-constexpr double iterations = 7000;
-
-DIFFUSION_TEST(EqualityFTCS) {
-  // set string from the header file
-  string test_path = testSimulationCSVDir;
-  MatrixXd reference = CSV2Eigen(test_path);
-  cout << "FTCS Test: " << endl;
-
-  Grid grid = setupSimulation(timestep, iterations); // Boundary
-  Boundary bc = Boundary(grid);
-
-  // Simulation
-
-  Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc);
-  // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
-  sim.setTimestep(timestep);
-  sim.setIterations(iterations);
-  sim.run();
-
-  cout << endl;
-  EXPECT_TRUE(checkSimilarity(reference, grid.getConcentrations(), 0.1));
-}
-
-DIFFUSION_TEST(EqualityBTCS) {
-  // set string from the header file
-  string test_path = testSimulationCSVDir;
-  MatrixXd reference = CSV2Eigen(test_path);
-  cout << "BTCS Test: " << endl;
-
-  Grid grid = setupSimulation(timestep, iterations); // Boundary
-  Boundary bc = Boundary(grid);
-
-  // Simulation
-  Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc);
-  // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
-  sim.setTimestep(timestep);
-  sim.setIterations(iterations);
-  sim.run();
-
-  cout << endl;
-  EXPECT_TRUE(checkSimilarityV2(reference, grid.getConcentrations(), 0.01));
-}
-
-DIFFUSION_TEST(InitializeEnvironment) {
-  int rc = 12;
-  Grid64 grid(rc, rc);
-  Boundary boundary(grid);
-
-  EXPECT_NO_THROW(Simulation sim(grid, boundary));
-}
-
-DIFFUSION_TEST(SimulationEnvironment) {
-  int rc = 12;
-  Grid64 grid(rc, rc);
-  Boundary boundary(grid);
-  Simulation<double, tug::FTCS_APPROACH> sim(grid, boundary);
-
-  EXPECT_EQ(sim.getIterations(), -1);
-
-  EXPECT_NO_THROW(sim.setIterations(2000));
-  EXPECT_EQ(sim.getIterations(), 2000);
-  EXPECT_THROW(sim.setIterations(-300), std::invalid_argument);
-
-  EXPECT_NO_THROW(sim.setTimestep(0.1));
-  EXPECT_DOUBLE_EQ(sim.getTimestep(), 0.1);
-  EXPECT_THROW(sim.setTimestep(-0.3), std::invalid_argument);
-}
-
-DIFFUSION_TEST(ClosedBoundaries) {
-
-  constexpr std::uint32_t nrows = 5;
-  constexpr std::uint32_t ncols = 5;
-
-  tug::Grid64 grid(nrows, ncols);
-
-  auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
-  auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
-  auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
-
-  grid.setConcentrations(concentrations);
-  grid.setAlpha(alphax, alphay);
-
-  tug::Boundary bc(grid);
-  bc.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1.0);
-  bc.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1.0);
-  bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0);
-  bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0);
-
-  tug::Simulation<double> sim(grid, bc);
-  sim.setTimestep(1);
-  sim.setIterations(1);
-
-  MatrixXd input_values(concentrations);
-  sim.run();
-
-  EXPECT_TRUE(checkSimilarityV2(input_values, grid.getConcentrations(), 1E-12));
-}
-
-DIFFUSION_TEST(ConstantInnerCell) {
-  constexpr std::uint32_t nrows = 5;
-  constexpr std::uint32_t ncols = 5;
-
-  tug::Grid64 grid(nrows, ncols);
-
-  auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
-  auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
-  auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
-
-  grid.setConcentrations(concentrations);
-  grid.setAlpha(alphax, alphay);
-
-  tug::Boundary bc(grid);
-  // inner
-  bc.setInnerBoundary(2, 2, 0);
-
-  tug::Simulation<double> sim(grid, bc);
-  sim.setTimestep(1);
-  sim.setIterations(1);
-
-  MatrixXd input_values(concentrations);
-  sim.run();
-
-  EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 2), 0);
-  EXPECT_LT(grid.getConcentrations().sum(), input_values.sum());
-
-  EXPECT_FALSE((grid.getConcentrations().array() > 1.0).any());
-
-  EXPECT_FALSE((grid.getConcentrations().array() < 0.0).any());
-}