Merge branch '16-include-advection-into-upstream-version' into 'main'

Draft: Resolve "Include advection into upstream version"

Closes #16

See merge request naaice/tug!37

.gitlab-ci.yml

@@ -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/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
     Diffusion

examples/Advection.cpp

@@ -0,0 +1,59 @@
+#include "tug/Boundary.hpp"
+#include <cstddef>
+#include <iostream>
+#include <tug/tug.hpp>
+
+using namespace Eigen;
+using namespace tug;
+
+int main(int argc, char *argv[]) {
+  constexpr std::size_t row = 5;
+  constexpr std::size_t col = 5;
+
+  RowMajMat<double> hydHeads = RowMajMat<double>::Constant(row, col, 1);
+  RowMajMat<double> concentrations = RowMajMat<double>::Constant(row, col, 0);
+
+  Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
+             tug::HYDRAULIC_RESOLVE::EXPLICIT>
+      velocities(hydHeads);
+
+  velocities.setDomain(5, 5);
+  velocities.setPermKX(RowMajMat<double>::Constant(row, col, 3E-7));
+  velocities.setPermKY(RowMajMat<double>::Constant(row, col, 3E-7));
+  velocities.setEpsilon(1E-8);
+
+  Advection advection = Advection(concentrations, velocities);
+
+  advection.setPorosity(RowMajMat<double>::Constant(row, col, 0.2));
+  advection.setIterations(3);
+  // 1 hour
+  advection.setTimestep(1.6666e+06);
+
+  // create boundaries
+  Boundary<double> &bcH = velocities.getBoundaryConditions();
+  bcH.setBoundarySideConstant(BC_SIDE_LEFT, 10);
+  bcH.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+  // bcH.setBoundarySideConstant(BC_SIDE_TOP, 1);
+  // bcH.setBoundarySideConstant(BC_SIDE_BOTTOM, 1);
+  // bcH.setInnerBoundary(row / 2, col / 2, 10);
+
+  Boundary<double> &bcC = advection.getBoundaryConditions();
+  bcC.setBoundarySideConstant(BC_SIDE_LEFT, 1);
+  bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+  // bcC.setInnerBoundary(row / 2, col / 2, 1);
+  // bcC.setBoundarySideConstant(BC_SIDE_LEFT, 0);
+  // bcC.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+  // bcC.setBoundarySideConstant(BC_SIDE_TOP, 0);
+  // bcC.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
+
+  advection.run();
+
+  std::cout << velocities.getConcentrationMatrix() << std::endl << std::endl;
+
+  std::cout << velocities.getVelocitiesX() << std::endl
+            << std::endl
+            << velocities.getVelocitiesY() << std::endl
+            << std::endl;
+
+  std::cout << concentrations << std::endl;
+}

examples/BTCS_2D_proto_example.cpp

@@ -1,5 +1,5 @@
 #include <Eigen/Eigen>
-#include <tug/Diffusion.hpp>
+#include <tug/tug.hpp>
 
 using namespace Eigen;
 using namespace tug;

examples/CMakeLists.txt

@@ -1,7 +1,9 @@
-add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
-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)
+# add_executable(BTCS_2D_proto_example BTCS_2D_proto_example.cpp)
+# 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)
+add_executable(Advection Advection.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)
+# 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)
+target_link_libraries(Advection tug)

examples/FTCS_2D_proto_closed_mdl.cpp

@@ -9,7 +9,7 @@
 #include <Eigen/Eigen>
 #include <cstdlib>
 #include <iostream>
-#include <tug/Diffusion.hpp>
+#include <tug/tug.hpp>
 
 using namespace Eigen;
 using namespace tug;

examples/FTCS_2D_proto_example_mdl.cpp

@@ -7,7 +7,7 @@
  */
 
 #include <Eigen/Eigen>
-#include <tug/Diffusion.hpp>
+#include <tug/tug.hpp>
 
 using namespace Eigen;
 using namespace tug;

include/tug/Advection/Advection.hpp

@@ -0,0 +1,408 @@
+/**
+ * @file Advection.hpp
+ * @brief API of Advection class, holding information for a simulation of
+ * advection. Holds a predifined Grid object, Boundary object and Velocities
+ * object
+ */
+
+#pragma once
+
+#include "tug/Core/Matrix.hpp"
+#include <Eigen/src/Core/Array.h>
+#include <Eigen/src/Core/AssignEvaluator.h>
+#include <Eigen/src/Core/Matrix.h>
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <stdlib.h>
+#include <string>
+#include <tug/Boundary.hpp>
+
+#include <Eigen/Eigen>
+#include <tug/Core/Numeric/BTCS.hpp>
+#include <tug/Core/Numeric/FTCS.hpp>
+#include <tug/Core/TugUtils.hpp>
+#include <tug/Diffusion/Diffusion.hpp>
+
+#include <tug/Advection/Velocities.hpp>
+
+using namespace Eigen;
+namespace tug {
+template <class T, HYDRAULIC_MODE hyd_mode, HYDRAULIC_RESOLVE hyd_resolve>
+class Advection : public BaseSimulationGrid<T> {
+private:
+  T timestep{-1};
+  int numThreads{omp_get_num_procs()};
+
+  Velocities<T, hyd_mode, hyd_resolve> &velocities;
+  RowMajMat<T> porosity;
+
+public:
+  /**
+   * @brief Construct a new Advection object, used to calculate material
+   * transport. A timestep and number of iterations must be set. A transient
+   * case can be selected by initializing Steady=false. With each timestep the
+   * Velocities object will also be updated.
+   * A steady case can be selected by initializing Steady=true. The velocities
+   * object will not be updated. Velocities can be simulated to convergence
+   * beforehand. Porosity can be set, the default is 1. CSV Output is off by
+   * default.
+   *
+   * @param grid Valid grid object
+   * @param bc Valid Boundary object
+   * @param Steady Used to choose between Steady and Transient case. Either true
+   * or false
+   */
+  Advection(RowMajMat<T> &origin,
+            Velocities<T, hyd_mode, hyd_resolve> &_velocities)
+      : BaseSimulationGrid<T>(origin), velocities(_velocities) {
+    tug_assert(origin.rows() == velocities.getConcentrationMatrix().rows() &&
+                   origin.cols() == velocities.getConcentrationMatrix().cols(),
+               "Advection grid and Velocities must have the same dimensions");
+  };
+
+  Advection(T *data, std::size_t rows, std::size_t cols,
+            Velocities<T, hyd_mode, hyd_resolve> &_velocities)
+      : BaseSimulationGrid<T>(data, rows, cols), velocities(_velocities) {
+    tug_assert(rows == velocities.getConcentrationMatrix().rows() &&
+                   cols == velocities.getConcentrationMatrix().cols(),
+               "Advection grid and Velocities must have the same dimensions");
+  };
+
+  /**
+   * @brief Sets the porosity of the medium
+   *
+   * @param porosity new porosity value
+   */
+  void setPorosity(const RowMajMat<T> &porosity) {
+    tug_assert(porosity.rows() == this->rows() &&
+                   porosity.cols() == this->cols(),
+               "Porosity matrix must have the same dimensions as the grid");
+    tug_assert(porosity.minCoeff() >= 0 && porosity.maxCoeff() <= 1,
+               "Porosity must be a value between 0 and 1 (inclusive)");
+
+    this->porosity = porosity;
+  }
+
+  /**
+   * @brief Set the size of the timestep. Must be greater than zero
+   *
+   * @param timestep Size of the timestep
+   */
+  void setTimestep(T timestep) {
+    tug_assert(timestep > 0, "Timestep must be greater than zero");
+    this->timestep = timestep;
+  }
+
+  /**
+   * @brief Set the number of desired openMP Threads.
+   *
+   * @param num_threads Number of desired threads. Must have a value between
+   *                    1 and the maximum available number of processors. The
+   * maximum number of processors is set as the default case during Advection
+   * construction.
+   */
+  void setNumberThreads(int num_threads) {
+    if (num_threads > 0 && num_threads <= omp_get_num_procs()) {
+      this->numThreads = num_threads;
+    } else {
+      int maxThreadNumber = omp_get_num_procs();
+      if (num_threads > maxThreadNumber) {
+        throw std::invalid_argument(
+            "Number of threads exceeds the number of processor cores (" +
+            std::to_string(maxThreadNumber) + ").");
+      } else {
+        throw std::invalid_argument("Number of threads is less than 1.");
+      }
+    }
+  }
+
+  void run() {
+    this->setDomain(velocities.domainX(), velocities.domainY());
+
+    this->applyInnerBoundaries();
+
+    if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) {
+      velocities.run();
+    }
+
+    for (int i = 0; i < this->getIterations(); i++) {
+      if constexpr (hyd_mode == HYDRAULIC_MODE::TRANSIENT) {
+        velocities.run();
+      }
+      startAdvection();
+    }
+  }
+
+private:
+  /**
+   * @brief function calculating material transport for one timestep
+   */
+  void startAdvection() {
+    const std::size_t rows = this->rows();
+    const std::size_t cols = this->cols();
+    const T dx = this->deltaCol();
+    const T dy = this->deltaRow();
+    const T volume = this->deltaCol() * this->deltaRow();
+
+    RowMajMatMap<T> &newConcentrations = this->getConcentrationMatrix();
+
+    RowMajMat<T> veloX = this->velocities.getVelocitiesX();
+    RowMajMat<T> veloY = this->velocities.getVelocitiesY();
+    const Boundary<T> &bc = this->getBoundaryConditions();
+
+    // Calculate Courant-Levy-Frederich condition
+    const T maxFx = std::max(abs(veloX.maxCoeff()), abs(veloX.minCoeff()));
+    const T maxFy = std::max(abs(veloY.maxCoeff()), abs(veloY.minCoeff()));
+    const T maxFlux = std::max(maxFx, maxFy);
+
+    tug_assert(maxFlux != 0, "Division by zero: maxF is zero.");
+
+    // TODO: sustitute 1 by porisity
+    const RowMajMat<T> porevolume =
+        RowMajMat<T>::Constant(rows, cols, volume) * 1;
+
+    const T cfl = (porevolume / maxFlux).minCoeff();
+    const int innerSteps = (int)ceil(timestep / cfl);
+    const T innerTimestep = timestep / innerSteps;
+
+    std::cout << "CFL (adv) timestep: " << cfl << std::endl;
+    std::cout << "Inner iterations (adv): " << innerSteps << std::endl;
+    std::cout << "Inner timestep (adv): " << innerTimestep << std::endl;
+
+    // const RowMajMat<T> multiplier = porevolume * (1 / innerTimestep);
+
+    RowMajMat<T> fluxX(rows, cols + 1);
+    RowMajMat<T> fluxY(rows + 1, cols);
+
+    for (std::size_t k = 0; k < innerSteps; k++) {
+      const RowMajMat<T> oldConcentrations = newConcentrations;
+
+      // Update flux with concentrations in x-direction
+      for (std::size_t row_i = 0; row_i < fluxX.rows(); row_i++) {
+        for (std::size_t col_i = 0; col_i < fluxX.cols(); col_i++) {
+          T &currentFlux = fluxX(row_i, col_i);
+          const T &currentVelo = veloX(row_i, col_i);
+          const std::int32_t cellIndex = (currentVelo > 0) ? col_i - 1 : col_i;
+
+          if (cellIndex < 0 || cellIndex >= cols) {
+            const auto bcElement = bc.getBoundaryElement(
+                cellIndex < 0 ? BC_SIDE_LEFT : BC_SIDE_RIGHT, row_i);
+            // if we have a constant boundary, we use the value of the boundary
+            if (bcElement.getType() == BC_TYPE_CONSTANT) {
+              currentFlux = currentVelo * bcElement.getValue();
+              continue;
+            }
+
+            // otherwise it is a closed boundary
+            currentFlux = 0;
+            continue;
+          }
+          // otherwise we use the concentration of the inflow/outflow cell,
+          // multiplied by the velocity
+          currentFlux = currentVelo * oldConcentrations(row_i, cellIndex);
+        }
+      }
+
+      // Update flux with concentrations in y-direction
+      for (std::size_t row_i = 0; row_i < fluxY.rows(); row_i++) {
+        for (std::size_t col_i = 0; col_i < fluxY.cols(); col_i++) {
+          T &currentFlux = fluxY(row_i, col_i);
+          const T &currentVelo = veloY(row_i, col_i);
+          const std::int32_t cellIndex = (currentVelo > 0) ? row_i - 1 : row_i;
+
+          if (cellIndex < 0 || cellIndex >= rows) {
+            const auto bcElement = bc.getBoundaryElement(
+                cellIndex < 0 ? BC_SIDE_TOP : BC_SIDE_BOTTOM, col_i);
+            // if we have a constant boundary, we use the value of the
+            // boundary
+            if (bcElement.getType() == BC_TYPE_CONSTANT) {
+              currentFlux = currentVelo * bcElement.getValue();
+              continue;
+            }
+
+            // otherwise it is a closed boundary
+            currentFlux = 0;
+            continue;
+          }
+          // otherwise we use the concentration of the inflow/outflow cell,
+          // multiplied by the velocity
+          currentFlux = currentVelo * oldConcentrations(cellIndex, col_i);
+        }
+      }
+
+      // Update concentrations
+      for (std::size_t row_i = 0; row_i < rows; row_i++) {
+        for (std::size_t col_i = 0; col_i < cols; col_i++) {
+          const T horizontalFlux =
+              fluxX(row_i, col_i) - fluxX(row_i, col_i + 1);
+          const T verticalFlux = fluxY(row_i, col_i) - fluxY(row_i + 1, col_i);
+          newConcentrations(row_i, col_i) = oldConcentrations(row_i, col_i) +
+                                            innerTimestep /
+                                                porevolume(row_i, col_i) *
+                                                (horizontalFlux + verticalFlux);
+        }
+      }
+
+      // const Boundary<T> leftCorner = bc.getBoundaryElement(BC_SIDE_LEFT, 0);
+      // newConcentrations(0, 0) = oldConcentrations(0, 0) + ;
+      //
+      // for (std::size_t col_i = 0; col_i < cols; col_i++) {
+      //   T horizontalFlux = oldConcentrations(0, col_i);
+      //   if (col_i == 0) {
+      //     switch (bc.getBoundaryElementType(BC_SIDE_LEFT, 0)) {
+      //     case BC_TYPE_CLOSED: {
+      //       break;
+      //     }
+      //     case BC_TYPE_CONSTANT: {
+      //       newConcentrations(0, 0) =
+      //           oldConcentrations(0, 0) +
+      //           (bc.getBoundaryElementValue(BC_SIDE_LEFT, 0) * veloX(0, 0)) *
+      //               (2 * innerTimestep) / dx;
+      //       newConcentrations(rows - 1, 0) =
+      //           oldConcentrations(rows - 1, 0) +
+      //           (bc.getBoundaryElementValue(BC_SIDE_LEFT, 0) *
+      //            veloX(rows - 1, 0)) *
+      //               (2 * innerTimestep) / dx;
+      //       break;
+      //     }
+      //     }
+      //     newConcentrations(0, 0) -=
+      //         veloX(0, 1) * oldConcentrations(0, 1) * innerTimestep / dx;
+      //   }
+      // }
+      //
+      // for (std::size_t row_i = 1; row_i < rows - 1; row_i++) {
+      //   for (std::size_t col_i = 1; col_i < cols - 1; col_i++) {
+      //     const T fluxHorizontal =
+      //         veloX(row_i, col_i) * oldConcentrations(row_i, col_i - 1) -
+      //         veloX(row_i, col_i + 1) * oldConcentrations(row_i, col_i);
+      //     const T fluxVertical =
+      //         veloY(row_i, col_i) * oldConcentrations(row_i - 1, col_i) -
+      //         veloY(row_i + 1, col_i) * oldConcentrations(row_i, col_i);
+      //     newConcentrations(row_i, col_i) =
+      //         oldConcentrations(row_i, col_i) +
+      //         (fluxHorizontal / dx + fluxVertical / dy) * innerTimestep;
+      //   }
+      // }
+
+      // Calculate sum of incoming/outgoing Flow*Concentration in x-direction in
+      // each cell #pragma omp parallel for num_threads(numThreads)
+      //       for (int i = 0; i < rows; i++) {
+      //         for (int j = 0; j < cols + 1; j++) {
+      //           if (j == 0) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_LEFT, i) !=
+      //             BC_TYPE_CLOSED) {
+      //               if (veloX(i, j) > 0) {
+      //                 // outx positive -> flow from border to cell i,j
+      //                 newConcentrations(i, j) +=
+      //                     veloX(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_LEFT, i);
+      //               } else if (veloX(i, j) < 0) {
+      //                 // outx negative -> flow from i,j towards border
+      //                 newConcentrations(i, j) += veloX(i, j) *
+      //                 oldConcentrations(i, j);
+      //               }
+      //             }
+      //           } else if (j == cols) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_RIGHT, i) !=
+      //             BC_TYPE_CLOSED) {
+      //               if (veloX(i, j) > 0) {
+      //                 // outx positive-> flow from i,j-1 towards border
+      //                 newConcentrations(i, j - 1) -=
+      //                     veloX(i, j) * oldConcentrations(i, j - 1);
+      //               } else if (veloX(i, j) < 0) {
+      //                 // outx negative -> flow from border to cell i,j-1
+      //                 newConcentrations(i, j - 1) -=
+      //                     veloX(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_LEFT, i);
+      //               }
+      //             }
+      //           }
+      //           // flow between inner cells
+      //           else {
+      //             // outx positive -> flow from cell i,j-1 towards cell i,j
+      //             if (veloX(i, j) > 0) {
+      //               newConcentrations(i, j - 1) -=
+      //                   veloX(i, j) * oldConcentrations(i, j - 1);
+      //               newConcentrations(i, j) +=
+      //                   veloX(i, j) * oldConcentrations(i, j - 1);
+      //             }
+      //             // outx negative -> flow from cell i,j toward cell i,j-1
+      //             else if (veloX(i, j) < 0) {
+      //               newConcentrations(i, j - 1) +=
+      //                   veloX(i, j) * oldConcentrations(i, j);
+      //               newConcentrations(i, j) += veloX(i, j) *
+      //               oldConcentrations(i, j);
+      //             }
+      //           }
+      //         }
+      //       }
+      // // calculate sum in y-direction
+      // // parallelize outer loop over columns to ensure thread-safety, each
+      // thread only
+      // // modifies cells within its column
+      // #pragma omp parallel for num_threads(numThreads)
+      //       for (int j = 0; j < cols; j++) {
+      //         for (int i = 0; i < rows + 1; i++) {
+      //           if (i == 0) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_TOP, j) !=
+      //             BC_TYPE_CLOSED) {
+      //               if (veloY(i, j) > 0) {
+      //                 // outy positive -> flow from border to cell i,j
+      //                 newConcentrations(i, j) +=
+      //                     veloY(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_TOP, j);
+      //               } else if (veloY(i, j) < 0) {
+      //                 // outy negative -> flow from i,j towards border
+      //                 newConcentrations(i, j) += veloY(i, j) *
+      //                 oldConcentrations(i, j);
+      //               }
+      //             }
+      //           } else if (i == rows) {
+      //             if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, j) !=
+      //                 BC_TYPE_CLOSED) {
+      //               if (veloY(i, j) > 0) {
+      //                 // outy positive-> flow from i-1,j towards border
+      //                 newConcentrations(i - 1, j) -=
+      //                     veloY(i, j) * oldConcentrations(i - 1, j);
+      //               } else if (veloY(i, j) < 0) {
+      //                 // outy negative -> flow from border to cell i,j-1
+      //                 newConcentrations(i - 1, j) -=
+      //                     veloY(i, j) *
+      //                     bc.getBoundaryElementValue(BC_SIDE_BOTTOM, j);
+      //               }
+      //             }
+      //           }
+      //           // flow between inner cells
+      //           else {
+      //             // outy positive -> flow from cell i-1,j towards cell i,j
+      //             if (veloY(i, j) > 0) {
+      //               newConcentrations(i - 1, j) -=
+      //                   veloY(i, j) * oldConcentrations(i - 1, j);
+      //               newConcentrations(i, j) +=
+      //                   veloY(i, j) * oldConcentrations(i - 1, j);
+      //             }
+      //             // outy negative -> flow from cell i,j toward cell i-1,j
+      //             else if (veloY(i, j) < 0) {
+      //               newConcentrations(i - 1, j) -=
+      //                   veloY(i, j) * oldConcentrations(i, j);
+      //               newConcentrations(i, j) += veloY(i, j) *
+      //               oldConcentrations(i, j);
+      //             }
+      //           }
+      //         }
+      //       }
+
+      // for (std::size_t row_i = 0; row_i < rows; row_i++) {
+      //   for (std::size_t col_i = 0; col_i < cols; col_i++) {
+      //     newConcentrations(row_i, col_i) =
+      //         oldConcentrations(row_i, col_i) +
+      //         newConcentrations(row_i, col_i) * multiplier(row_i, col_i);
+      //   }
+      // }
+    }
+  }
+};
+} // namespace tug

include/tug/Advection/Velocities.hpp

@@ -0,0 +1,376 @@
+/**
+ * @file Velocities.hpp
+ * @brief API of Velocities class, holding information for a simulation of
+ * Hydraulic Charge and Darcy-Velocities. Holds a predifined Grid object and
+ * Boundary object.
+ *
+ */
+
+#pragma once
+
+#include "tug/Core/Numeric/SimulationInput.hpp"
+#include <algorithm>
+#include <cstddef>
+#include <iostream>
+#include <stdexcept>
+#include <stdlib.h>
+#include <string>
+
+#include <tug/Boundary.hpp>
+#include <tug/Core/BaseSimulation.hpp>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Core/Numeric/BTCS.hpp>
+#include <tug/Core/Numeric/FTCS.hpp>
+#include <tug/Core/TugUtils.hpp>
+#include <tug/Diffusion/Diffusion.hpp>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_num_procs() 1
+#endif
+
+using namespace Eigen;
+namespace tug {
+
+enum class HYDRAULIC_MODE { TRANSIENT, STEADY_STATE };
+enum class HYDRAULIC_RESOLVE { EXPLICIT, IMPLICIT };
+
+template <class T, HYDRAULIC_MODE hyd_mode, HYDRAULIC_RESOLVE hyd_resolve>
+class Velocities : public BaseSimulationGrid<T> {
+private:
+  int innerIterations{1};
+  T timestep{-1};
+  T epsilon{1E-5};
+  int numThreads{omp_get_num_procs()};
+  bool steady{false};
+
+  RowMajMat<T> velocitiesX;
+  RowMajMat<T> velocitiesY;
+
+  RowMajMat<T> permKX;
+  RowMajMat<T> permKY;
+
+public:
+  Velocities(RowMajMat<T> &origin)
+      : BaseSimulationGrid<T>(origin),
+        velocitiesX(origin.rows(), origin.cols() + 1),
+        velocitiesY(origin.rows() + 1, origin.cols()),
+        permKX(origin.rows(), origin.cols()),
+        permKY(origin.rows(), origin.cols()) {};
+
+  Velocities(T *data, std::size_t rows, std::size_t cols)
+      : BaseSimulationGrid<T>(data, rows, cols), velocitiesX(rows, cols + 1),
+        velocitiesY(rows + 1, cols), permKX(rows, cols), permKY(rows, cols) {};
+
+  // Velocities(T *data, std::size_t length)
+  //     : BaseSimulationGrid<T>(data, 1, length), velocitiesX(1, length + 1),
+  //       alphaX(1, length) {};
+
+  /**
+   * @brief Set the epsilon value, the relativ error allowed for convergence
+   *
+   * @param epsilon the new epsilon value
+   */
+  void setEpsilon(T epsilon) {
+    tug_assert(0 <= epsilon && epsilon < 1,
+               "Relative Error epsilon must be between 0 and 1");
+
+    this->epsilon = epsilon;
+  }
+
+  /**
+   * @brief Get the alphaX matrix.
+   *
+   * @return RowMajMat<T>& Reference to the alphaX matrix.
+   */
+  RowMajMat<T> &getPermKX() { return permKX; }
+
+  /**
+   * @brief Get the alphaY matrix.
+   *
+   * @return RowMajMat<T>& Reference to the alphaY matrix.
+   */
+  RowMajMat<T> &getPermKY() {
+    tug_assert(
+        this->getDim(),
+        "Grid is not two dimensional, there is no domain in y-direction!");
+
+    return permKY;
+  }
+
+  /**
+   * @brief Set the alphaX matrix.
+   *
+   * @param alphaX The new alphaX matrix.
+   */
+  void setPermKX(const RowMajMat<T> &alphaX) { this->permKX = alphaX; }
+
+  /**
+   * @brief Set the alphaY matrix.
+   *
+   * @param alphaY The new alphaY matrix.
+   */
+  void setPermKY(const RowMajMat<T> &alphaY) {
+    tug_assert(
+        this->getDim(),
+        "Grid is not two dimensional, there is no domain in y-direction!");
+
+    this->permKY = alphaY;
+  }
+
+  bool isSteady() const { return steady; }
+
+  /**
+   * @brief Set the timestep per iteration
+   *
+   * @param timestep timestep per iteration
+   */
+  void setTimestep(T timestep) override {
+    if (timestep <= 0) {
+      throw std::invalid_argument("Timestep must be greater than zero");
+    }
+    this->timestep = timestep;
+    const T deltaColSquare = this->deltaCol() * this->deltaCol();
+    const T deltaRowSquare = this->deltaRow() * this->deltaRow();
+    const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
+
+    const T maxK = std::max(this->permKX.maxCoeff(), this->permKY.maxCoeff());
+
+    T cfl = minDeltaSquare / (4 * maxK);
+    if (timestep > cfl) {
+      this->innerIterations = (int)ceil(timestep / cfl);
+      this->timestep = timestep / (double)innerIterations;
+      std::cerr << "Warning :: Timestep was adjusted, because of stability "
+                   "conditions. Time duration was approximately preserved by "
+                   "adjusting internal number of iterations."
+                << std::endl;
+      std::cout << "FTCS" << "_" << "2D" << " :: Required "
+                << this->innerIterations
+                << " inner iterations with dt=" << this->timestep << std::endl;
+    } else {
+      this->innerIterations = 1;
+    }
+  };
+
+  /**
+   * @brief Set the number of desired openMP Threads.
+   *
+   * @param num_threads Number of desired threads. Must have a value between
+   *                    1 and the maximum available number of processors. The
+   * maximum number of processors is set as the default case during Velocities
+   * construction.
+   */
+  void setNumberThreads(int num_threads) {
+    if (num_threads > 0 && num_threads <= omp_get_num_procs()) {
+      this->numThreads = num_threads;
+    } else {
+      int maxThreadNumber = omp_get_num_procs();
+      throw std::invalid_argument(
+          "Number of threads exceeds the number of processor cores (" +
+          std::to_string(maxThreadNumber) + ") or is less than 1.");
+    }
+  }
+
+  /**
+   * @brief Getter function for outx, the matrix containing velocities in
+   * x-Direction; returns a reference to outx
+   *
+   * */
+  const RowMajMat<T> &getVelocitiesX() const { return this->velocitiesX; }
+
+  /**
+   * @brief Getter function for outy, the matrix containing velocities in
+   * y-Direction; return a reference to outy
+   */
+  const RowMajMat<T> &getVelocitiesY() const { return this->velocitiesY; }
+
+  /**
+   * @brief Simulation of hydraulic charge either until convergence,
+   * or for a number of selected timesteps. Calculation of Darcy-velocities.
+   */
+  void run() override {
+    // if iterations < 1 calculate hydraulic charge until steady state is
+    // reached
+
+    this->applyInnerBoundaries();
+
+    SimulationInput<T> input = {.concentrations =
+                                    this->getConcentrationMatrix(),
+                                .alphaX = this->getPermKX(),
+                                .alphaY = this->getPermKY(),
+                                .boundaries = this->getBoundaryConditions(),
+                                .dim = this->getDim(),
+                                .timestep = this->timestep,
+                                .rowMax = this->rows(),
+                                .colMax = this->cols(),
+                                .deltaRow = this->deltaRow(),
+                                .deltaCol = this->deltaCol()};
+
+    if constexpr (hyd_mode == HYDRAULIC_MODE::STEADY_STATE) {
+      if (steady) {
+        return;
+      }
+
+      const T deltaColSquare = this->deltaCol() * this->deltaCol();
+      const T deltaRowSquare = this->deltaRow() * this->deltaRow();
+      const T minDeltaSquare = std::min(deltaColSquare, deltaRowSquare);
+
+      const T maxK = std::max(this->permKX.maxCoeff(), this->permKY.maxCoeff());
+      // Calculate largest possible timestep, depending on K and gridsize
+
+      std::cout << "CFL (hydHead) timestep: " << minDeltaSquare / (4 * maxK)
+                << std::endl;
+
+      setTimestep(minDeltaSquare / (4 * maxK));
+
+      input.timestep = this->timestep;
+
+      RowMajMat<T> oldConcentrations;
+      do {
+        oldConcentrations = this->getConcentrationMatrix();
+        (void)calculate_hydraulic_flow(input);
+      } while (!checkConvergance(oldConcentrations));
+
+      steady = true;
+
+    } else {
+      if (timestep == -1) {
+        throw_invalid_argument("Timestep is not set");
+      }
+      for (int i = 0; i < innerIterations; i++) {
+        (void)calculate_hydraulic_flow(input);
+      }
+    }
+
+    (void)computeFluidVelocities();
+  };
+
+private:
+  /**
+   * @brief Calculate the new hydraulic charge using FTCS
+   */
+  void calculate_hydraulic_flow(SimulationInput<T> &sim_in) {
+    if constexpr (hyd_resolve == HYDRAULIC_RESOLVE::EXPLICIT) {
+      FTCS_2D(sim_in, numThreads);
+    } else {
+      BTCS_2D(sim_in, ThomasAlgorithm, numThreads);
+    }
+  };
+
+  /**
+   * @brief checks if the matrix of Hydraulic Heads has converged
+   *
+   * @return bool true if for all corresponding cells of the matrices,
+   * containing old and new Charge values, the relative error is below the
+   * selected Epsilon
+   */
+  bool checkConvergance(const RowMajMat<T> &oldHeads) {
+    const auto abs_err = (oldHeads - this->getConcentrationMatrix()).cwiseAbs();
+    const auto rel_err = abs_err.cwiseQuotient(this->getConcentrationMatrix());
+
+    return rel_err.maxCoeff() < epsilon;
+  }
+
+  /**
+   * @brief Update the matrices containing Darcy velocities in x and y
+   * directions
+   */
+  void computeFluidVelocities() {
+    const std::size_t rows = this->rows();
+    const std::size_t cols = this->cols();
+    const T dx = this->deltaCol();
+    const T dy = this->deltaRow();
+    const RowMajMat<T> &hydraulicCharges = this->getConcentrationMatrix();
+
+    const RowMajMat<T> &permKX = this->permKX;
+    const RowMajMat<T> &permKY = this->permKY;
+
+    const Boundary<T> &bc = this->getBoundaryConditions();
+
+    // calculate velocities in x-direction
+    for (std::size_t i_rows = 0; i_rows < rows; i_rows++) {
+      const auto bc_left = bc.getBoundaryElement(BC_SIDE_LEFT, i_rows);
+      switch (bc_left.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesX(i_rows, 0) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesX(i_rows, 0) =
+            -this->permKX(i_rows, 0) *
+            (hydraulicCharges(i_rows, 0) - bc_left.getValue()) / (dx / 2);
+        break;
+      }
+      }
+
+      const auto bc_right = bc.getBoundaryElement(BC_SIDE_RIGHT, i_rows);
+      switch (bc_right.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesX(i_rows, cols) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesX(i_rows, cols) =
+            -permKX(i_rows, cols - 1) *
+            (bc_right.getValue() - hydraulicCharges(i_rows, cols - 1)) /
+            (dx / 2);
+        break;
+      }
+      }
+    }
+
+// main loop for calculating velocities in x-direction for inner cells
+#pragma omp parallel for num_threads(numThreads)
+    for (int i = 0; i < rows; i++) {
+      for (int j = 1; j < cols; j++) {
+        velocitiesX(i, j) =
+            -permKX(i, j - 1) *
+            (hydraulicCharges(i, j) - hydraulicCharges(i, j - 1)) / dx;
+      }
+    }
+
+    // calculate velocities in y-direction
+    for (std::size_t i_cols = 0; i_cols < cols; i_cols++) {
+      const auto bc_top = bc.getBoundaryElement(BC_SIDE_TOP, i_cols);
+      switch (bc_top.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesY(0, i_cols) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesY(0, i_cols) =
+            -permKY(0, i_cols) *
+            (hydraulicCharges(0, i_cols) - bc_top.getValue()) / (dy / 2);
+        break;
+      }
+      }
+
+      const auto bc_bottom = bc.getBoundaryElement(BC_SIDE_BOTTOM, i_cols);
+      switch (bc_bottom.getType()) {
+      case BC_TYPE_CLOSED: {
+        velocitiesY(rows, i_cols) = 0;
+        break;
+      }
+      case BC_TYPE_CONSTANT: {
+        velocitiesY(rows, i_cols) =
+            -permKY(rows - 1, i_cols) *
+            (bc_bottom.getValue() - hydraulicCharges(rows - 1, i_cols)) /
+            (dy / 2);
+        break;
+      }
+      }
+    }
+
+// main loop for calculating velocities in y-direction for inner cells
+#pragma omp parallel for num_threads(numThreads)
+    for (int i = 1; i < rows; i++) {
+      for (int j = 0; j < cols; j++) {
+        velocitiesY(i, j) =
+            -permKY(i, j) *
+            (hydraulicCharges(i, j) - hydraulicCharges(i - 1, j)) / dy;
+      }
+    }
+  };
+};
+} // namespace tug

include/tug/Boundary.hpp

@@ -7,12 +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>
 
@@ -45,7 +46,7 @@ public:
    *        BC_TYPE_CLOSED, where the value takes -1 and does not hold any
    *        physical meaning.
    */
-  BoundaryElement(){};
+  BoundaryElement() {};
 
   /**
    * @brief Construct a new Boundary Element object for the constant case.
@@ -115,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
@@ -124,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)){};
-
-  /**
-   * @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()) {
+      : dim(rows == 1 ? 1 : 2), cols(cols), rows(rows) {
     if (this->dim == 1) {
       this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(
           2); // in 1D only left and right boundary
@@ -153,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.

include/tug/Core/BaseSimulation.hpp

@@ -1,6 +1,10 @@
 #pragma once
 
-#include <stdexcept>
+#include "tug/Boundary.hpp"
+#include <cstddef>
+#include <cstdint>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Core/TugUtils.hpp>
 
 namespace tug {
 
@@ -8,11 +12,11 @@ namespace tug {
  * @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
+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 */
 };
 
@@ -20,30 +24,278 @@ enum CSV_OUTPUT {
  * @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 */
+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 TIME_MEASURE {
-  TIME_MEASURE_OFF, /*!< do not print any time measures */
-  TIME_MEASURE_ON   /*!< print time measure after last iteration */
+enum class TIME_MEASURE {
+  OFF, /*!< do not print any time measures */
+  ON   /*!< print time measure after last iteration */
 };
 
-class BaseSimulation {
-protected:
-  CSV_OUTPUT csv_output{CSV_OUTPUT_OFF};
-  CONSOLE_OUTPUT console_output{CONSOLE_OUTPUT_OFF};
-  TIME_MEASURE time_measure{TIME_MEASURE_OFF};
+/**
+ * @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;
+
+protected:
+  void applyInnerBoundaries() {
+    const auto &inner_bc = boundaryConditions.getInnerBoundaries();
+    if (inner_bc.empty()) {
+      return;
+    }
+
+    for (const auto &[rowcol, value] : inner_bc) {
+      concentrationMatrix(rowcol.first, rowcol.second) = value;
+    }
+  }
 
 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.
    *
@@ -58,13 +310,16 @@ public:
    *                     - 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!");
-    }
+  void setOutputCSV(CSV_OUTPUT csv_output) { this->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
@@ -80,14 +335,20 @@ public:
    * matrices to console
    */
   void setOutputConsole(CONSOLE_OUTPUT console_output) {
-    if (console_output < CONSOLE_OUTPUT_OFF &&
-        console_output > CONSOLE_OUTPUT_VERBOSE) {
-      throw std::invalid_argument("Invalid console output option given!");
-    }
-
     this->console_output = console_output;
   }
 
+  /**
+   * @brief 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.
    *
@@ -98,13 +359,16 @@ public:
    * console
    */
   void setTimeMeasure(TIME_MEASURE time_measure) {
-    if (time_measure < TIME_MEASURE_OFF && time_measure > TIME_MEASURE_ON) {
-      throw std::invalid_argument("Invalid time measure option given!");
-    }
-
     this->time_measure = time_measure;
   }
 
+  /**
+   * @brief 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.
@@ -112,10 +376,9 @@ public:
    * @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.");
-    }
+    tug_assert(iterations > 0,
+               "Number of iterations must be greater than zero.");
+
     this->iterations = iterations;
   }
 
@@ -131,5 +394,7 @@ public:
    *        parameters.
    */
   virtual void run() = 0;
+
+  virtual void setTimestep(T timestep) = 0;
 };
-} // namespace tug
+} // namespace tug

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.getCol();
-  T sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol());
+  const std::size_t &length = input.colMax;
+  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.getAlphaX();
+  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();
-  int colMax = grid.getCol();
-  T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
-  T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
+  const std::size_t &rowMax = input.rowMax;
+  const std::size_t &colMax = input.colMax;
+  const T sx = input.timestep / (2 * input.deltaCol * input.deltaCol);
+  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();
-  RowMajMat<T> alphaY = grid.getAlphaY();
+  const RowMajMat<T> &alphaX = input.alphaX;
+  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();
-  alphaY.transposeInPlace();
+  const RowMajMat<T> alphaX_t = alphaX.transpose();
+  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);
-    b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
-                             bcLeft, bcRight, inner_bc, rowMax, i, sy, sx);
+    A = createCoeffMatrix(alphaY_t, bcTop, bcBottom, inner_bc, rowMax, i, sy);
+    b = createSolutionVector(concentrations_t1, alphaY_t, alphaX_t, bcTop,
+                             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>
-void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
-  if (grid.getDim() == 1) {
-    BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
-  } else if (grid.getDim() == 2) {
-    BTCS_2D(grid, bc, timestep, EigenLUAlgorithm, numThreads);
+template <class T> void BTCS_LU(SimulationInput<T> &input, int numThreads) {
+  tug_assert(input.dim <= 2,
+             "Error: Only 1- and 2-dimensional grids are defined!");
+
+  if (input.dim == 1) {
+    BTCS_1D(input, EigenLUAlgorithm);
   } else {
-    throw_invalid_argument(
-        "Error: Only 1- and 2-dimensional grids are defined!");
+    BTCS_2D(input.dim, EigenLUAlgorithm, numThreads);
   }
 }
 
 // entry point for Thomas algorithm solver; differentiate 1D and 2D grid
-template <class T>
-void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
-  if (grid.getDim() == 1) {
-    BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
-  } else if (grid.getDim() == 2) {
-    BTCS_2D(grid, bc, timestep, ThomasAlgorithm, numThreads);
+template <class T> void BTCS_Thomas(SimulationInput<T> &input, int numThreads) {
+  tug_assert(input.dim <= 2,
+             "Error: Only 1- and 2-dimensional grids are defined!");
+
+  if (input.dim == 1) {
+    BTCS_1D(input, ThomasAlgorithm);
   } else {
-    throw_invalid_argument(
-        "Error: Only 1- and 2-dimensional grids are defined!");
+    BTCS_2D(input, ThomasAlgorithm, numThreads);
   }
 }
 } // namespace tug

include/tug/Core/Numeric/FTCS.hpp

@@ -8,11 +8,13 @@
 #ifndef FTCS_H_
 #define FTCS_H_
 
-#include "../TugUtils.hpp"
 #include "tug/Core/Matrix.hpp"
-
+#include <cstddef>
 #include <cstring>
 #include <tug/Boundary.hpp>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Core/Numeric/SimulationInput.hpp>
+#include <tug/Core/TugUtils.hpp>
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -22,384 +24,223 @@
 
 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) {
+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 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);
+  return alpha_center_left * conc_left -
+         (alpha_center_left + alpha_center_right) * conc_c +
+         alpha_center_right * conc_right;
 }
 
-// calculates vertical change on one cell independent of boundary type
 template <class T>
-static inline T calcVerticalChange(Grid<T> &grid, int &row, int &col) {
+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();
 
-  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!");
+  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;
   }
-}
-
-// 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!");
+  case BC_TYPE_CLOSED: {
+    return (alpha_center_neighbor * (conc_neighbor - conc_c));
   }
-}
-
-// 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!");
-  }
+  tug_assert(false, "Undefined Boundary Condition Type!");
+  return 0;
 }
 
 // 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();
+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;
 
-  RowMajMat<T> &concentrations_grid = grid.getConcentrations();
+  RowMajMatMap<T> &concentrations_grid = input.concentrations;
+
+  const auto &alphaX = input.alphaX;
+  const auto &bc = input.boundaries;
   // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(1, colMax, 0);
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
 
   // only one row in 1D case -> row constant at index 0
   int row = 0;
 
+  const auto inner_bc = bc.getInnerBoundaryRow(0);
+
   // inner cells
   // independent of boundary condition type
   for (int col = 1; col < colMax - 1; col++) {
-    concentrations_t1(row, col) = concentrations_grid(row, col) +
-                                  timestep / (deltaCol * deltaCol) *
-                                      (calcHorizontalChange(grid, row, col));
+    if (inner_bc[col].first) {
+      continue;
+    }
+    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;
-  concentrations_t1(row, col) =
-      concentrations_grid(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeLeftBoundary(grid, bc, row, col));
+  if (!inner_bc[col].first) {
+    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
   col = colMax - 1;
-  concentrations_t1(row, col) =
-      concentrations_grid(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeRightBoundary(grid, bc, row, col));
+  if (!inner_bc[col].first) {
+    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
-  std::memcpy(concentrations_grid.data(), concentrations_t1.data(),
-              colMax * sizeof(T));
+  concentrations_grid = 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();
+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;
 
-  RowMajMat<T> &concentrations_grid = grid.getConcentrations();
+  RowMajMatMap<T> &concentrations_grid = input.concentrations;
+  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);
 
   // matrix for concentrations at time t+1
-  RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
+  RowMajMat<T> concentrations_t1 = concentrations_grid;
 
-  // 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) = concentrations_grid(row, col) +
-                                    timestep / (deltaRow * deltaRow) *
-                                        (calcVerticalChange(grid, row, col)) +
-                                    timestep / (deltaCol * deltaCol) *
-                                        (calcHorizontalChange(grid, row, col));
+  for (std::size_t row_i = 0; row_i < rowMax; row_i++) {
+    for (std::size_t col_i = 0; col_i < colMax; col_i++) {
+      if (bc.getInnerBoundary(row_i, col_i).first) {
+        continue;
+      }
+      // 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;
     }
   }
 
-  // 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) =
-        concentrations_grid(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) =
-        concentrations_grid(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) =
-        concentrations_grid(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) =
-        concentrations_grid(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) =
-      concentrations_grid(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) =
-      concentrations_grid(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) =
-      concentrations_grid(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) =
-      concentrations_grid(row, col) +
-      timestep / (deltaCol * deltaCol) *
-          (calcHorizontalChangeRightBoundary(grid, bc, row, col)) +
-      timestep / (deltaRow * deltaRow) *
-          (calcVerticalChangeBottomBoundary(grid, bc, row, col));
-
   // overwrite obsolete concentrations
-  std::memcpy(concentrations_grid.data(), concentrations_t1.data(),
-              rowMax * colMax * sizeof(T));
-  // }
+  concentrations_grid = 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);
+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 {
-    throw_invalid_argument(
-        "Error: Only 1- and 2-dimensional grids are defined!");
+    FTCS_2D(input, numThreads);
   }
 }
 } // namespace tug

include/tug/Core/Numeric/SimulationInput.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <cstdint>
+#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;
+  T timestep;
+  const std::size_t rowMax;
+  const std::size_t colMax;
+  const T deltaRow;
+  const T deltaCol;
+};
+} // namespace tug

include/tug/Diffusion/Diffusion.hpp

@@ -8,8 +8,7 @@
 
 #pragma once
 
-#include "Boundary.hpp"
-#include "Grid.hpp"
+#include "tug/Core/Matrix.hpp"
 #include <algorithm>
 #include <filesystem>
 #include <fstream>
@@ -18,10 +17,9 @@
 #include <string>
 #include <vector>
 
-#include "Core/Numeric/BTCS.hpp"
-#include "Core/Numeric/FTCS.hpp"
-#include "Core/TugUtils.hpp"
-#include "tug/Core/BaseSimulation.hpp"
+#include <tug/Core/BaseSimulation.hpp>
+#include <tug/Core/Numeric/BTCS.hpp>
+#include <tug/Core/Numeric/FTCS.hpp>
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -63,31 +61,95 @@ enum SOLVER {
  */
 template <class T, APPROACH approach = BTCS_APPROACH,
           SOLVER solver = THOMAS_ALGORITHM_SOLVER>
-class Diffusion : public BaseSimulation {
+class Diffusion : public BaseSimulationGrid<T> {
 private:
   T timestep{-1};
   int innerIterations{1};
   int numThreads{omp_get_num_procs()};
 
-  Grid<T> &grid;
-  Boundary<T> &bc;
+  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
-   * 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.
+   * @brief Construct a new Diffusion object from a given Eigen matrix
    *
-   * @param grid Valid grid object
-   * @param bc Valid boundary condition object
-   * @param approach Approach to solving the problem. Either FTCS or BTCS.
    */
-  Diffusion(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
+  Diffusion(RowMajMat<T> &origin) : BaseSimulationGrid<T>(origin) {
+    init_alpha();
+  }
+
+  /**
+   * @brief Construct a new 2D Diffusion object from a given data pointer and
+   * the dimensions.
+   *
+   */
+  Diffusion(T *data, int rows, int cols)
+      : BaseSimulationGrid<T>(data, rows, cols) {
+    init_alpha();
+  }
+
+  /**
+   * @brief Construct a new 1D Diffusion object from a given data pointer and
+   * the length.
+   *
+   */
+  Diffusion(T *data, std::size_t length) : BaseSimulationGrid<T>(data, length) {
+    init_alpha();
+  }
+
+  /**
+   * @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;
+  }
 
   /**
    * @brief Setting the time step for each iteration step. Time step must be
@@ -95,31 +157,31 @@ public:
    *
    * @param timestep Valid timestep greater than zero.
    */
-  void setTimestep(T timestep) {
-      tug_assert(timestep > 0, "Timestep has to be greater than zero.");
+  void setTimestep(T timestep) override {
+    tug_assert(timestep > 0, "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.getDeltaCol();
-        const T maxAlpha = grid.getAlphaX().maxCoeff();
+        const T deltaSquare = this->deltaCol();
+        const T maxAlpha = this->alphaX.maxCoeff();
 
         // Courant-Friedrichs-Lewy condition
         cfl = deltaSquare / (4 * maxAlpha);
-      } else if (grid.getDim() == 2) {
-        const T deltaColSquare = grid.getDeltaCol() * grid.getDeltaCol();
+      } else if (this->getDim() == 2) {
+        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
@@ -183,8 +245,8 @@ public:
    * @brief Outputs the current concentrations of the grid on the console.
    *
    */
-  inline void printConcentrationsConsole() const {
-    std::cout << grid.getConcentrations() << std::endl;
+  void printConcentrationsConsole() const {
+    std::cout << this->getConcentrationMatrix() << std::endl;
     std::cout << std::endl;
   }
 
@@ -204,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 col = std::to_string(grid.getCol());
-    std::string numIterations = std::to_string(iterations);
+    std::string row = std::to_string(this->rows());
+    std::string col = std::to_string(this->cols());
+    std::string numIterations = std::to_string(this->getIterations());
 
     std::string filename =
         approachString + "_" + row + "_" + col + "_" + numIterations + ".csv";
@@ -225,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)
@@ -260,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();
   }
@@ -269,30 +333,44 @@ public:
    * @brief Method starts the simulation process with the previously set
    *        parameters.
    */
-  void run() {
-    tug_assert(this->timestep > 0, "Timestep is not set!");
-    tug_assert(this->iterations > 0, "Number of iterations are not set!");
+  void run() override {
+    tug_assert(this->getTimestep() > 0, "Timestep is not set!");
+    tug_assert(this->getIterations() > 0, "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();
 
+    this->applyInnerBoundaries();
+
+    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()};
+
     if constexpr (approach == FTCS_APPROACH) { // FTCS case
-      for (int i = 0; i < iterations * innerIterations; i++) {
-        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 *
@@ -305,29 +383,28 @@ public:
     } else if constexpr (approach == BTCS_APPROACH) { // BTCS case
 
       if constexpr (solver == EIGEN_LU_SOLVER) {
-        for (int i = 0; i < iterations; i++) {
-          if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
+        for (int i = 0; i < this->getIterations(); i++) {
+          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++) {
-          if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
+        for (int i = 0; i < this->getIterations(); i++) {
+          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
 
@@ -339,22 +416,22 @@ public:
       RowMajMat<T> concentrations;
       RowMajMat<T> concentrationsFTCS;
       RowMajMat<T> concentrationsResult;
-      for (int i = 0; i < iterations * innerIterations; i++) {
-        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);
         }
 
-        concentrations = grid.getConcentrations();
+        concentrations = this->getConcentrationMatrix();
         FTCS(this->grid, this->bc, this->timestep, this->numThreads);
-        concentrationsFTCS = grid.getConcentrations();
-        grid.setConcentrations(concentrations);
-        BTCS_Thomas(this->grid, this->bc, this->timestep, this->numThreads);
-        concentrationsResult =
-            beta * concentrationsFTCS + (1 - beta) * grid.getConcentrations();
-        grid.setConcentrations(concentrationsResult);
+        concentrationsFTCS = this->getConcentrationMatrix();
+        this->getConcentrationMatrix() = concentrations;
+        BTCS_Thomas(sim_input, this->numThreads);
+        concentrationsResult = beta * concentrationsFTCS +
+                               (1 - beta) * this->getConcentrationMatrix();
+        this->getConcentrationMatrix() = concentrationsResult;
       }
     }
 
@@ -362,15 +439,15 @@ 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;
     }

include/tug/Grid.hpp

@@ -1,342 +0,0 @@
-#pragma once
-
-/**
- * @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 "tug/Core/TugUtils.hpp"
-#include <Eigen/Core>
-#include <Eigen/Sparse>
-#include <Eigen/src/Core/Matrix.h>
-#include <Eigen/src/Core/util/Constants.h>
-#include <cstddef>
-
-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 Construct a new Grid object.
-   *
-   * Constructs a new Grid object with given concentrations, defined by an
-   * Eigen::Matrix. The domain length is per default the same as the length. The
-   * alpha coefficients are set to 1. The dimensions of the grid are determined
-   * by the given matrix, which can also be an Eigen::Vector for a 1D-Grid.
-   *
-   * @param concentrations An Eigen3 MatrixX<T> holding the concentrations.
-   */
-  Grid(const RowMajMat<T> &concentrations) {
-    if (concentrations.rows() == 1) {
-      this->dim = 1;
-      this->domainCol = static_cast<T>(concentrations.cols());
-      this->deltaCol = static_cast<T>(this->domainCol) /
-                       static_cast<T>(concentrations.cols()); // -> 1
-
-      this->concentrations = concentrations;
-      return;
-    }
-
-    if (concentrations.cols() == 1) {
-      this->dim = 1;
-      this->domainCol = static_cast<T>(concentrations.rows());
-      this->deltaCol = static_cast<T>(this->domainCol) /
-                       static_cast<T>(concentrations.rows()); // -> 1
-
-      this->concentrations = concentrations.transpose();
-      return;
-    }
-
-    this->dim = 2;
-    this->domainRow = static_cast<T>(concentrations.rows());
-    this->domainCol = static_cast<T>(concentrations.cols());
-    this->deltaRow = static_cast<T>(this->domainRow) /
-                     static_cast<T>(concentrations.rows()); // -> 1
-    this->deltaCol = static_cast<T>(this->domainCol) /
-                     static_cast<T>(concentrations.cols()); // -> 1
-
-    this->concentrations = concentrations;
-    // this->alphaX = RowMajMat<T>::Constant(concentrations.rows(),
-    //                                       concentrations.cols(),
-    //                                       MAT_INIT_VAL);
-    // this->alphaY = RowMajMat<T>::Constant(concentrations.rows(),
-    //                                       concentrations.cols(),
-    //                                       MAT_INIT_VAL);
-  }
-
-  /**
-   * @brief Construct a new Grid object.
-   *
-   * Constructs a new 1D Grid object with given concentrations, defined by a
-   * pointer to consecutive memory and the length of the array. The domain
-   * length is per default the same as the count of grid cells (length of
-   * array). The memory region is mapped internally, changes will affect the
-   * original array and the memory shall not be freed. There is no check for
-   * correct memory size!
-   *
-   * @param concentrations Pointer to consecutive memory holding concentrations.
-   * @param length Length of the array/the 1D grid.
-   */
-  Grid(T *concentrations, std::size_t length) : dim(1) {
-    this->domainCol = static_cast<T>(length); // -> 1
-    this->deltaCol =
-        static_cast<T>(this->domainCol) / static_cast<T>(length); // -> 1
-
-    this->concentrations = RowMajMatMap<T>(concentrations, 1, length);
-  }
-
-  /**
-   * @brief Construct a new Grid object.
-   *
-   * Constructs a new 2D Grid object with given concentrations, defined by a
-   * pointer to consecutive memory and the number of rows and columns. The
-   * domain size is per default the same as the number of rows and columns. The
-   * memory region is mapped internally, changes will affect the original array
-   * and the memory shall not be freed. There is no check for correct memory
-   * size!
-   *
-   * @param concentrations Pointer to consecutive memory holding concentrations.
-   * @param row Number of rows.
-   * @param col Number of columns.
-   */
-  Grid(T *concentrations, std::size_t row, std::size_t col) : dim(2) {
-    this->domainRow = static_cast<T>(row); // -> 1
-    this->domainCol = static_cast<T>(col); // -> 1
-    this->deltaCol =
-        static_cast<T>(this->domainCol) / static_cast<T>(col); // -> 1
-    this->deltaRow =
-        static_cast<T>(this->domainRow) / static_cast<T>(row); // -> 1
-
-    this->concentrations = RowMajMatMap<T>(concentrations, row, col);
-  }
-
-  /**
-   * @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; }
-
-  void initAlpha() {
-    this->alphaX = RowMajMat<T>::Constant(
-        this->concentrations.rows(), this->concentrations.cols(), MAT_INIT_VAL);
-    if (dim > 1) {
-
-      this->alphaY =
-          RowMajMat<T>::Constant(this->concentrations.rows(),
-                                 this->concentrations.cols(), MAT_INIT_VAL);
-    }
-  }
-
-  /**
-   * @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) {
-    tug_assert(dim == 1,
-               "Grid is not one dimensional, use 2D setter function!");
-
-    tug_assert(
-        alpha.rows() == 1 && alpha.cols() == this->concentrations.cols(),
-        "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) {
-    tug_assert(dim == 1,
-               "Grid is not one dimensional, 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) {
-    tug_assert(dim == 2,
-               "Grid is not two dimensional, use 1D setter function!");
-
-    tug_assert(alphaX.rows() == this->concentrations.rows(),
-               "Alpha in x-direction "
-               "has wrong number of rows!");
-    tug_assert(alphaX.cols() == this->concentrations.cols(),
-               "Alpha in x-direction has wrong number of columns!");
-
-    tug_assert(alphaY.rows() == this->concentrations.rows(),
-               "Alpha in y-direction "
-               "has wrong number of rows!");
-
-    tug_assert(alphaY.cols() == this->concentrations.cols(),
-               "Alpha in y-direction has wrong number of columns!");
-
-    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) {
-    tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!");
-
-    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 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 {
-    tug_assert(this->alphaX.size() > 0, "AlphaX is empty!");
-    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 {
-    tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!");
-    tug_assert(this->alphaY.size() > 0, "AlphaY is empty!");
-
-    return this->alphaY;
-  }
-
-  /**
-   * @brief Gets the dimensions of the grid.
-   *
-   * @return Dimensions, either 1 or 2.
-   */
-  int getDim() const { return this->dim; }
-
-  /**
-   * @brief Gets the number of rows of the grid.
-   *
-   * @return Number of rows.
-   */
-  int getRow() const { return this->concentrations.rows(); }
-
-  /**
-   * @brief Gets the number of columns of the grid.
-   *
-   * @return Number of columns.
-   */
-  int getCol() const { return this->concentrations.cols(); }
-
-  /**
-   * @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) {
-    tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!");
-    tug_assert(domainLength > 0, "Given domain length is not positive!");
-
-    this->domainCol = domainLength;
-    this->deltaCol =
-        double(this->domainCol) / double(this->concentrations.cols());
-  }
-
-  /**
-   * @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) {
-    tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!");
-    tug_assert(domainCol > 0,
-               "Given domain size in x-direction is not positive!");
-    tug_assert(domainRow > 0,
-               "Given domain size in y-direction is not positive!");
-
-    this->domainRow = domainRow;
-    this->domainCol = domainCol;
-    this->deltaRow =
-        double(this->domainRow) / double(this->concentrations.rows());
-    this->deltaCol =
-        double(this->domainCol) / double(this->concentrations.cols());
-  }
-
-  /**
-   * @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 {
-    tug_assert(dim == 2, "Grid is not two dimensional, there is no delta in "
-                         "y-direction!");
-
-    return this->deltaRow;
-  }
-
-private:
-  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;     // 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

include/tug/tug.hpp

@@ -0,0 +1,7 @@
+#pragma once
+
+#include <tug/Advection/Advection.hpp>
+#include <tug/Advection/Velocities.hpp>
+#include <tug/Boundary.hpp>
+#include <tug/Core/Matrix.hpp>
+#include <tug/Diffusion/Diffusion.hpp>

test/CMakeLists.txt

@@ -12,7 +12,8 @@ FetchContent_MakeAvailable(googletest)
 add_executable(testTug 
 setup.cpp
 testDiffusion.cpp 
-testGrid.cpp 
+testVelocities.cpp
+testAdvection.cpp
 testFTCS.cpp 
 testBoundary.cpp
 )
@@ -26,4 +27,4 @@ get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH)
 # set relative path in header file
 configure_file(testSimulation.hpp.in testSimulation.hpp)
 # include test directory with generated header file from above
-target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src")
+target_include_directories(testTug PUBLIC "${CMAKE_CURRENT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}/src")

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<
-      Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
-      matrixEntries.data(), matrixRowNumber,
-      matrixEntries.size() / matrixRowNumber);
+  return tug::RowMajMatMap<double>(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,
-                              double maxDiff) {
+inline bool checkSimilarityV2(tug::RowMajMat<double> &a,
+                              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/testAdvection.cpp

@@ -0,0 +1,69 @@
+#include <gtest/gtest.h>
+#include <tug/tug.hpp>
+
+#define ADVECTION_TEST(x) TEST(Advection, x)
+
+ADVECTION_TEST(LeftToRight) {
+  constexpr std::size_t rows = 21;
+  constexpr std::size_t cols = 21;
+
+  constexpr double K = 1E-2;
+  constexpr double timestep = 5039.05;
+  constexpr std::size_t iterations = 21;
+  constexpr double porosity = 0.2;
+
+  constexpr double epsilon = 1E-13;
+
+  tug::RowMajMat<double> hydHeads =
+      tug::RowMajMat<double>::Constant(rows, cols, 1);
+  tug::RowMajMat<double> concentrations =
+      tug::RowMajMat<double>::Constant(rows, cols, 0);
+
+  tug::RowMajMat<double> permK =
+      tug::RowMajMat<double>::Constant(rows, cols, K);
+
+  tug::Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
+                  tug::HYDRAULIC_RESOLVE::IMPLICIT>
+      velocities(hydHeads);
+  velocities.setDomain(100, 100);
+  velocities.setPermKX(permK);
+  velocities.setPermKY(permK);
+  velocities.setEpsilon(1E-8);
+
+  tug::Advection advection(concentrations, velocities);
+
+  advection.setPorosity(tug::RowMajMat<double>::Constant(rows, cols, porosity));
+  advection.setIterations(iterations);
+  advection.setTimestep(timestep);
+
+  tug::Boundary<double> &bcH = velocities.getBoundaryConditions();
+  bcH.setBoundarySideConstant(tug::BC_SIDE_LEFT, 10);
+  bcH.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1);
+
+  tug::Boundary<double> &bcC = advection.getBoundaryConditions();
+  bcC.setBoundarySideConstant(tug::BC_SIDE_LEFT, 0.1);
+  bcC.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1);
+
+  advection.run();
+
+  // check if the concentration is transported from left to right
+  for (std::size_t i_rows = 0; i_rows < rows; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols < cols - 1; i_cols++) {
+      if (i_cols == 0) {
+        EXPECT_LE(concentrations(i_rows, i_cols), 10);
+      } else {
+        EXPECT_GE(concentrations(i_rows, i_cols),
+                  concentrations(i_rows, i_cols + 1));
+      }
+    }
+  }
+
+  // the values should also be equal from top to bottom
+  for (std::size_t i_cols = 0; i_cols < cols; i_cols++) {
+    const double &ref = concentrations(0, i_cols);
+    for (std::size_t i_rows = 1; i_rows < rows; i_rows++) {
+      // check if the values are equal within the epsilon range
+      EXPECT_NEAR(ref, concentrations(i_rows, i_cols), epsilon);
+    }
+  }
+}

test/testBoundary.cpp

@@ -28,12 +28,8 @@ BOUNDARY_TEST(Element) {
 }
 
 BOUNDARY_TEST(Class) {
-  Eigen::VectorXd conc(10);
-  Grid grid1D = Grid64(conc);
-  Eigen::MatrixXd conc2D(10, 12);
-  Grid grid2D = Grid64(conc2D);
-  Boundary boundary1D = Boundary(grid1D);
-  Boundary boundary2D = Boundary(grid2D);
+  Boundary<double> boundary1D(10);
+  Boundary<double> boundary2D(10, 12);
   vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));
 
   constexpr double inner_condition_value = -5;
@@ -47,7 +43,6 @@ 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),
@@ -76,7 +71,6 @@ BOUNDARY_TEST(Class) {
   }
 
   {
-    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);

test/testDiffusion.cpp

@@ -1,8 +1,9 @@
 #include "TestUtils.hpp"
+#include "tug/Core/Matrix.hpp"
 #include "gtest/gtest.h"
 #include <gtest/gtest.h>
 #include <stdexcept>
-#include <tug/Diffusion.hpp>
+#include <tug/Diffusion/Diffusion.hpp>
 
 #include <Eigen/src/Core/Matrix.h>
 #include <string>
@@ -16,18 +17,27 @@ using namespace Eigen;
 using namespace std;
 using namespace tug;
 
-Grid64 setupSimulation(double timestep, int iterations) {
-  int row = 11;
-  int col = 11;
+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
-  MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
+  // RowMajMat<double> concentrations = MatrixXd::Constant(row, col, 0);
   concentrations(5, 5) = 1;
 
-  Grid grid = Grid64(concentrations);
-  grid.setDomain(domain_row, domain_col);
+  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++) {
@@ -45,9 +55,10 @@ Grid64 setupSimulation(double timestep, int iterations) {
       alpha(i, j) = 0.1;
     }
   }
-  grid.setAlpha(alpha, alpha);
+  diffusiongrid.setAlphaX(alpha);
+  diffusiongrid.setAlphaY(alpha);
 
-  return grid;
+  return diffusiongrid;
 }
 
 constexpr double timestep = 0.001;
@@ -56,127 +67,194 @@ constexpr double iterations = 7000;
 DIFFUSION_TEST(EqualityFTCS) {
   // set string from the header file
   string test_path = testSimulationCSVDir;
-  MatrixXd reference = CSV2Eigen(test_path);
+  RowMajMat<double> reference = CSV2Eigen(test_path);
   cout << "FTCS Test: " << endl;
 
-  Grid grid = setupSimulation(timestep, iterations); // Boundary
-  Boundary bc = Boundary(grid);
+  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);
+  // Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
   // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
-  sim.setTimestep(timestep);
-  sim.setIterations(iterations);
+  // sim.setTimestep(timestep);
+  // sim.setIterations(iterations);
   sim.run();
 
   cout << endl;
-  EXPECT_TRUE(checkSimilarity(reference, grid.getConcentrations(), 0.1));
+  EXPECT_TRUE(checkSimilarity(reference, sim.getConcentrationMatrix(), 0.1));
 }
 
 DIFFUSION_TEST(EqualityBTCS) {
   // set string from the header file
   string test_path = testSimulationCSVDir;
-  MatrixXd reference = CSV2Eigen(test_path);
+  RowMajMat<double> reference = CSV2Eigen(test_path);
   cout << "BTCS Test: " << endl;
 
-  Grid grid = setupSimulation(timestep, iterations); // Boundary
-  Boundary bc = Boundary(grid);
+  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);
+  // Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
   // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
-  sim.setTimestep(timestep);
-  sim.setIterations(iterations);
+  // sim.setTimestep(timestep);
+  // sim.setIterations(iterations);
   sim.run();
 
   cout << endl;
-  EXPECT_TRUE(checkSimilarityV2(reference, grid.getConcentrations(), 0.01));
+  EXPECT_TRUE(checkSimilarityV2(reference, sim.getConcentrationMatrix(), 0.01));
 }
 
 DIFFUSION_TEST(InitializeEnvironment) {
   int rc = 12;
-  Eigen::MatrixXd concentrations(rc, rc);
-  Grid64 grid(concentrations);
-  Boundary boundary(grid);
+  RowMajMat<double> concentrations(rc, rc);
+  // Grid64 grid(concentrations);
+  // Boundary boundary(grid);
 
-  EXPECT_NO_THROW(Diffusion sim(grid, boundary));
+  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);
+// 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_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.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.*");
-}
+//   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;
 
-  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);
+  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);
 
-  tug::Grid64 grid(concentrations);
+  Diffusion<double> sim(concentrations);
+  sim.getAlphaX() = alphax;
+  sim.getAlphaY() = alphay;
 
-  grid.setAlpha(alphax, alphay);
+  // tug::Grid64 grid(concentrations);
 
-  tug::Boundary bc(grid);
+  // 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);
+  // tug::Diffusion<double> sim(grid, bc);
   sim.setTimestep(1);
   sim.setIterations(1);
 
-  MatrixXd input_values(concentrations);
+  RowMajMat<double> input_values(concentrations);
   sim.run();
 
-  EXPECT_TRUE(checkSimilarityV2(input_values, grid.getConcentrations(), 1E-12));
+  EXPECT_TRUE(
+      checkSimilarityV2(input_values, sim.getConcentrationMatrix(), 1E-12));
 }
 
 DIFFUSION_TEST(ConstantInnerCell) {
   constexpr std::uint32_t nrows = 5;
   constexpr std::uint32_t ncols = 5;
 
-  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);
+  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);
 
-  tug::Grid64 grid(concentrations);
-  grid.setAlpha(alphax, alphay);
+  Diffusion<double> sim(concentrations);
+  sim.getAlphaX() = alphax;
+  sim.getAlphaY() = alphay;
 
-  tug::Boundary bc(grid);
+  // 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);
+  // tug::Diffusion<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());
+  const auto &concentrations_result = sim.getConcentrationMatrix();
 
-  EXPECT_FALSE((grid.getConcentrations().array() > 1.0).any());
+  EXPECT_DOUBLE_EQ(concentrations_result(2, 2), 0);
+  EXPECT_LT(concentrations_result.sum(), input_values.sum());
 
-  EXPECT_FALSE((grid.getConcentrations().array() < 0.0).any());
+  EXPECT_FALSE((concentrations_result.array() > 1.0).any());
+
+  EXPECT_FALSE((concentrations_result.array() < 0.0).any());
+}
+
+DIFFUSION_TEST(Symmetry) {
+  // Arrange
+  constexpr std::size_t rows = 25;
+  constexpr std::size_t cols = 25;
+
+  constexpr std::size_t center_row = rows / 2;
+  constexpr std::size_t center_col = cols / 2;
+
+  tug::RowMajMat<double> concentrations =
+      tug::RowMajMat<double>::Constant(rows, cols, 1);
+
+  tug::RowMajMat<double> alpha =
+      tug::RowMajMat<double>::Constant(rows, cols, 1E-2);
+
+  tug::Diffusion<double, tug::FTCS_APPROACH> sim(concentrations);
+
+  sim.setDomain(100, 100);
+  sim.setAlphaX(alpha);
+  sim.setAlphaY(alpha);
+  // choose a high number of iterations, which lead to small changes in ULP
+  // between symmetric cells
+  sim.setIterations(10000);
+  sim.setTimestep(10);
+
+  tug::Boundary<double> &bcH = sim.getBoundaryConditions();
+
+  bcH.setInnerBoundary(center_row, center_col, 10);
+
+  sim.run();
+
+  // check symmetry
+  for (std::size_t i_rows = 0; i_rows <= center_row; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols <= center_col; i_cols++) {
+      if (i_rows == center_row && i_cols == center_col) {
+        continue;
+      }
+      // to avoid floating point errors, we check with ASSERT_DOUBLE_EQ with a
+      // precision of ULP(4), see https://stackoverflow.com/a/4149599
+      ASSERT_DOUBLE_EQ(concentrations(i_rows, i_cols),
+                       concentrations(rows - i_rows - 1, cols - i_cols - 1));
+    }
+  }
 }

test/testGrid.cpp

@@ -1,264 +0,0 @@
-#include "gtest/gtest.h"
-#include <Eigen/Core>
-#include <tug/Grid.hpp>
-
-#include <gtest/gtest.h>
-
-using namespace Eigen;
-using namespace std;
-using namespace tug;
-
-#define GRID_TEST(x) TEST(Grid, x)
-
-GRID_TEST(Grid64OneDimensional) {
-  int l = 12;
-  Eigen::VectorXd conc(l);
-  Grid64 grid(conc);
-  grid.initAlpha();
-
-  {
-    EXPECT_EQ(grid.getDim(), 1);
-    EXPECT_EQ(grid.getCol(), 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.getAlphaX().rows(), 1);
-    EXPECT_EQ(grid.getAlphaX().cols(), l);
-    EXPECT_EQ(grid.getDeltaCol(), 1);
-
-    EXPECT_DEATH(grid.getAlphaY(), ".* no alphaY!.*");
-    EXPECT_DEATH(grid.getDeltaRow(), ".* not two dimensional, .*");
-  }
-
-  {
-    // correct alpha matrix
-    MatrixXd alpha = MatrixXd::Constant(1, l, 3);
-    EXPECT_NO_THROW(grid.setAlpha(alpha));
-
-    EXPECT_DEATH(grid.setAlpha(alpha, alpha), ".* is not two dimensional, .*");
-
-    grid.setAlpha(alpha);
-    EXPECT_EQ(grid.getAlphaX(), alpha);
-    EXPECT_NO_THROW(grid.getAlphaX());
-    EXPECT_DEATH(grid.getAlphaY(), ".* no alphaY!.*");
-
-    // false alpha matrix
-    MatrixXd wAlpha = MatrixXd::Constant(3, l, 2);
-    EXPECT_DEATH(grid.setAlpha(wAlpha), ".* mismatch with Grid dimensions!.*");
-  }
-
-  {
-    int d = 8;
-    // set 1D domain
-    EXPECT_NO_THROW(grid.setDomain(d));
-
-    // set 2D domain
-    EXPECT_DEATH(grid.setDomain(d, d), ".* not two dimensional, .*");
-
-    grid.setDomain(d);
-    EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(d) / double(l));
-    EXPECT_DEATH(grid.getDeltaRow(), ".* not two dimensional, .*");
-
-    // set too small domain
-    EXPECT_DEATH(grid.setDomain(-2), "Given domain length .*");
-  }
-}
-
-GRID_TEST(Grid64Quadratic) {
-  int rc = 12;
-  Eigen::MatrixXd conc(rc, rc);
-  Grid64 grid(conc);
-  grid.initAlpha();
-
-  {
-    EXPECT_EQ(grid.getDim(), 2);
-    EXPECT_EQ(grid.getCol(), rc);
-    EXPECT_EQ(grid.getRow(), rc);
-
-    EXPECT_EQ(grid.getConcentrations().rows(), rc);
-    EXPECT_EQ(grid.getConcentrations().cols(), rc);
-
-    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 alpha matrices
-    MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
-    MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
-    EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
-
-    EXPECT_DEATH(grid.setAlpha(alphax), ".* 2D setter function!.*");
-
-    grid.setAlpha(alphax, alphay);
-    EXPECT_EQ(grid.getAlphaX(), alphax);
-    EXPECT_EQ(grid.getAlphaY(), alphay);
-
-    // false alpha matrices
-    alphax = MatrixXd::Constant(rc + 3, rc + 1, 3);
-    EXPECT_DEATH(grid.setAlpha(alphax, alphay),
-                 ".*has wrong number of rows!.*");
-    alphay = MatrixXd::Constant(rc + 2, rc + 1, 3);
-    EXPECT_DEATH(grid.setAlpha(alphax, alphay),
-                 ".*has wrong number of rows!.*");
-  }
-
-  {
-    int dr = 8;
-    int dc = 9;
-
-    // set 1D domain
-    EXPECT_DEATH(grid.setDomain(dr), ".* 2D domain setter!.*");
-
-    // 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_DEATH(grid.setDomain(dr, dc), ".* not positive!.*");
-    dr = 8;
-    dc = 0;
-    EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*");
-    dr = -2;
-    EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*");
-  }
-}
-
-GRID_TEST(Grid64NonQuadratic) {
-  int r = 12;
-  int c = 15;
-  Eigen::MatrixXd conc(r, c);
-  Grid64 grid(conc);
-  grid.initAlpha();
-
-  {
-    EXPECT_EQ(grid.getDim(), 2);
-    EXPECT_EQ(grid.getCol(), c);
-    EXPECT_EQ(grid.getRow(), r);
-
-    EXPECT_EQ(grid.getConcentrations().rows(), r);
-    EXPECT_EQ(grid.getConcentrations().cols(), c);
-
-    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 alpha matrices
-    MatrixXd alphax = MatrixXd::Constant(r, c, 2);
-    MatrixXd alphay = MatrixXd::Constant(r, c, 4);
-    EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
-
-    grid.setAlpha(alphax, alphay);
-    EXPECT_EQ(grid.getAlphaX(), alphax);
-    EXPECT_EQ(grid.getAlphaY(), alphay);
-  }
-
-  {
-    int dr = 8;
-    int dc = 9;
-
-    // set 2D domain
-    EXPECT_NO_THROW(grid.setDomain(dr, dc));
-
-    grid.setDomain(dr, dc);
-    EXPECT_EQ(grid.getDeltaCol(), double(dc) / double(c));
-    EXPECT_EQ(grid.getDeltaRow(), double(dr) / double(r));
-  }
-
-  {
-    int r = 4;
-    int c = 5;
-    std::vector<double> concentrations(r * c);
-
-    for (int i = 0; i < r * c; i++) {
-      concentrations[i] = i;
-    }
-    Grid64 grid(concentrations.data(), r, c);
-    grid.initAlpha();
-
-    {
-      EXPECT_EQ(grid.getDim(), 2);
-      EXPECT_EQ(grid.getCol(), c);
-      EXPECT_EQ(grid.getRow(), r);
-
-      EXPECT_EQ(grid.getConcentrations().rows(), r);
-      EXPECT_EQ(grid.getConcentrations().cols(), c);
-
-      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 alpha matrices
-      MatrixXd alphax = MatrixXd::Constant(r, c, 2);
-      MatrixXd alphay = MatrixXd::Constant(r, c, 4);
-      EXPECT_NO_THROW(grid.setAlpha(alphax, alphay));
-
-      EXPECT_DEATH(grid.setAlpha(alphax), ".* 2D setter function!.*");
-
-      grid.setAlpha(alphax, alphay);
-      EXPECT_EQ(grid.getAlphaX(), alphax);
-      EXPECT_EQ(grid.getAlphaY(), alphay);
-
-      // false alpha matrices
-      alphax = MatrixXd::Constant(r + 3, c + 1, 3);
-      EXPECT_DEATH(grid.setAlpha(alphax, alphay),
-                   ".*has wrong number of rows!.*");
-      alphay = MatrixXd::Constant(r + 2, c + 1, 5);
-      EXPECT_DEATH(grid.setAlpha(alphax, alphay),
-                   ".*has wrong number of rows!.*");
-
-      {
-        int dr = 8;
-        int dc = 9;
-
-        // set 1D domain
-        EXPECT_DEATH(grid.setDomain(dr), ".* 2D domain setter!.*");
-
-        // 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));
-      }
-
-      {
-        auto &concentrations = grid.getConcentrations();
-
-        for (int i = 0; i < r; i++) {
-          for (int j = 0; j < c; j++) {
-            concentrations(i, j) = i * c + j;
-          }
-        }
-
-        EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0);
-        EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1);
-        EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c);
-
-        EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0);
-        EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1);
-        EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c);
-        EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 1), 2 * c + 1);
-      }
-    }
-  }
-}

test/testVelocities.cpp

@@ -0,0 +1,78 @@
+#include "tug/Boundary.hpp"
+#include "tug/Core/Matrix.hpp"
+#include <cstddef>
+#include <tug/Advection/Velocities.hpp>
+
+#include <gtest/gtest.h>
+
+#define VELOCITIES_TEST(x) TEST(Velocities, x)
+
+VELOCITIES_TEST(SteadyStateCenter) {
+  // Arrange
+  constexpr std::size_t rows = 25;
+  constexpr std::size_t cols = 25;
+
+  constexpr std::size_t center_row = rows / 2;
+  constexpr std::size_t center_col = cols / 2;
+
+  constexpr double K = 1E-2;
+
+  tug::RowMajMat<double> hydHeads =
+      tug::RowMajMat<double>::Constant(rows, cols, 1);
+
+  tug::RowMajMat<double> permK =
+      tug::RowMajMat<double>::Constant(rows, cols, K);
+
+  tug::Velocities<double, tug::HYDRAULIC_MODE::STEADY_STATE,
+                  tug::HYDRAULIC_RESOLVE::EXPLICIT>
+      velo(hydHeads);
+
+  velo.setDomain(100, 100);
+  velo.setPermKX(permK);
+  velo.setPermKY(permK);
+
+  tug::Boundary<double> &bcH = velo.getBoundaryConditions();
+
+  bcH.setInnerBoundary(center_row, center_col, 10);
+
+  velo.run();
+
+  const auto &velocitiesX = velo.getVelocitiesX();
+  const auto &velocitiesY = velo.getVelocitiesY();
+
+  // Assert
+
+  // check velocities in x-direction
+  for (std::size_t i_rows = 0; i_rows < rows; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols < cols + 1; i_cols++) {
+      if (i_rows <= center_row && i_cols <= center_col) {
+        EXPECT_LE(velocitiesX(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols > center_col) {
+        EXPECT_GE(velocitiesX(i_rows, i_cols), 0);
+      } else if (i_rows <= center_row && i_cols > center_col) {
+        EXPECT_GE(velocitiesX(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols <= center_col) {
+        EXPECT_LE(velocitiesX(i_rows, i_cols), 0);
+      } else {
+        FAIL() << "Uncovered case";
+      }
+    }
+  }
+
+  // check velocities in y-direction
+  for (std::size_t i_rows = 0; i_rows < rows + 1; i_rows++) {
+    for (std::size_t i_cols = 0; i_cols < cols; i_cols++) {
+      if (i_rows <= center_row && i_cols <= center_col) {
+        EXPECT_LE(velocitiesY(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols > center_col) {
+        EXPECT_GE(velocitiesY(i_rows, i_cols), 0);
+      } else if (i_rows <= center_row && i_cols > center_col) {
+        EXPECT_LE(velocitiesY(i_rows, i_cols), 0);
+      } else if (i_rows > center_row && i_cols <= center_col) {
+        EXPECT_GE(velocitiesY(i_rows, i_cols), 0);
+      } else {
+        FAIL() << "Uncovered case";
+      }
+    }
+  }
+}