Merge branch 'ml/bump-cmake-version' into 'main'

relax Eigen3 version constraint

See merge request naaice/tug!45

.gitlab-ci.yml

@@ -1,34 +1,33 @@
-image: git.gfz-potsdam.de:5000/naaice/tug:ci
+image: gcc:14
+
+before_script:
+  - apt-get update && apt-get install -y cmake ninja-build libeigen3-dev git
 
 stages:
-    - build
     - test
+    - release
     - static_analyze
-    - doc
-
-build_release:
-  stage: build
-  artifacts:
-      paths:
-          - build/test/testTug
-      expire_in: 100s
-  script:
-      - mkdir build && cd build
-      - cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON ..
-      - make -j$(nproc)
-      - cp ../test/FTCS_11_11_7000.csv test/
 
 test:
     stage: test
     script:
-        - ./build/test/testTug
+      - mkdir build && cd build
+      - cmake -DCMAKE_BUILD_TYPE=Debug -DTUG_ENABLE_TESTING=ON -G Ninja ..
+      - ninja 
+      - ctest --output-junit test_results.xml
+    artifacts:
+      when: always
+      paths:
+        - build/test_results.xml
+      reports:
+        junit: build/test_results.xml
 
-pages:
-  stage: doc
+doc:
+  stage: release
   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
@@ -40,6 +39,22 @@ pages:
   rules:
     - if: $CI_COMMIT_REF_NAME == $CI_DEFAULT_BRANCH
 
+push:
+  stage: release
+  variables:
+    GITHUB_REPOSITORY: 'git@github.com:POET-Simulator/tug.git'
+    ORIGINAL_REPO_URL: 'https://git.gfz-potsdam.de/naaice/tug.git'
+    ORIGINAL_REPO_NAME: 'tug'
+  before_script:
+    # I know that there is this file env variable in gitlab, but somehow it does not work for me (still complaining about white spaces ...)
+    # Therefore, the ssh key is stored as a base64 encoded string
+    - mkdir -p ~/.ssh && echo $GITHUB_SSH_PRIVATE_KEY | base64 -d > ~/.ssh/id_ed25519 && chmod 0600 ~/.ssh/id_ed25519
+    - ssh-keyscan github.com >> ~/.ssh/known_hosts
+  script:
+    - rm -rf $ORIGINAL_REPO_NAME.git
+    - git clone --mirror $ORIGINAL_REPO_URL "$ORIGINAL_REPO_NAME.git" && cd $ORIGINAL_REPO_NAME.git
+    - git push --mirror $GITHUB_REPOSITORY
+
 lint:
   stage: static_analyze
   before_script:

CITATION.cff

@@ -0,0 +1,43 @@
+# This CITATION.cff file was generated with cffinit.
+# Visit https://bit.ly/cffinit to generate yours today!
+
+cff-version: 1.2.0
+title: 'tug: Transport on Uniform Grids'
+message: >-
+  If you use this software, please cite it using the
+  metadata from this file.
+type: software
+authors:
+  - given-names: Marco
+    family-names: De Lucia
+    email: delucia@gfz.de
+    orcid: 'https://orcid.org/0009-0000-3058-8472'
+  - given-names: Max
+    family-names: Lübke
+    email: mluebke@uni-potsdam.de
+    orcid: 'https://orcid.org/0009-0008-9773-3038'
+  - given-names: Bettina
+    family-names: Schnor
+    email: schnor@cs.uni-potsdam.de
+    orcid: 'https://orcid.org/0000-0001-7369-8057'
+  - given-names: Hannes
+    family-names: Signer
+    email: signer@uni-potsdam.de
+    orcid: 'https://orcid.org/0009-0000-3058-8472'
+  - given-names: Philipp
+    family-names: Ungrund
+    email: ungrund@uni-potsdam.de
+    orcid: 'https://orcid.org/0009-0007-0182-4051'
+repository-code: 'https://github.com/POET-Simulator/tug'
+abstract: >-
+  tug implements different numerical approaches for
+  transport problems, notably diffusion with implicit BTCS
+  (Backward Time, Central Space) Euler and parallel 2D ADI
+  (Alternating Direction Implicit).
+keywords:
+  - diffusion
+  - advection
+  - BTCS
+  - FTCS
+  - Thomas-Algorithm
+license: GPL-2.0

CMakeLists.txt

@@ -1,40 +1,23 @@
 # debian stable (currently bullseye)
-cmake_minimum_required(VERSION 3.18)
+cmake_minimum_required(VERSION 3.20)
 
 project(
   tug
   VERSION 0.4
   LANGUAGES CXX)
 
-set(CMAKE_CXX_STANDARD 17)
-
 find_package(Eigen3 REQUIRED NO_MODULE)
 find_package(OpenMP)
 
 include(GNUInstallDirs)
 
-# find_package(easy_profiler) option(EASY_OPTION_LOG "Verbose easy_profiler" 1)
-
-# SET(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -O2 -mfma")
 option(TUG_USE_OPENMP "Compile tug with OpenMP support" ON)
 
-set(CMAKE_CXX_FLAGS_GENERICOPT
-    "-O3 -march=native"
-    CACHE STRING "Flags used by the C++ compiler during opt builds." FORCE)
-
-set(CMAKE_BUILD_TYPE
-    "${CMAKE_BUILD_TYPE}"
-    CACHE
-      STRING
-      "Choose the type of build, options are: None Debug Release RelWithDebInfo MinSizeRel GenericOpt."
-      FORCE)
-
 option(
   TUG_USE_UNSAFE_MATH_OPT "Use compiler options to break IEEE compliances by
     oenabling reordering of instructions when adding/multiplying of floating
     points." OFF)
 
-
 option(TUG_ENABLE_TESTING "Run tests after succesfull compilation" OFF)
 
 option(TUG_HANNESPHILIPP_EXAMPLES "Compile example applications" OFF)
@@ -90,6 +73,7 @@ install(FILES "${PROJECT_BINARY_DIR}/${PROJECT_NAME}Config.cmake"
 install(DIRECTORY ${PROJECT_SOURCE_DIR}/include/tug DESTINATION include)
 
 if(TUG_ENABLE_TESTING)
+enable_testing()
   add_subdirectory(test)
 endif()
 

README.md

@@ -0,0 +1,103 @@
+![tug boat](./doc/images/tug_logo_small.png)
+
+`tug` implements different numerical approaches for transport
+problems, notably diffusion with implicit BTCS (Backward Time, Central
+Space) Euler and parallel 2D ADI (Alternating Direction Implicit).
+
+# About
+
+This project aims to provide a library for solving transport problems -
+diffusion, advection - on uniform grids implemented in C++. The library
+is built on top of
+[Eigen](https://eigen.tuxfamily.org/index.php?title=Main_Page),
+providing easy access to its optimized data structures and linear
+equation solvers.
+
+We designed the API to be as flexible as possible. Nearly every
+built-in, framework or third-party data structure can be used to model a
+problem, as long a pointer to continuous memory can be provided. We also
+provide parallelization using [OpenMP](https://www.openmp.org/), which
+can be easily turned on/off at compile time.
+
+At the current state, both 1D and 2D diffusion problems on a regular
+grid with constant alpha for all grid cells can be solved reliably.
+
+# Requirements
+
+- C++17 compliant compiler
+- [CMake](https://cmake.org/) >= 3.20
+- [Eigen](https://eigen.tuxfamily.org/) >= 3.4.0
+
+# Getting started
+
+`tug` is designed as a framework library and it relies on
+[CMake](https://cmake.org/) for building. If you already use
+`CMake` as your build toolkit for your application, you\'re
+good to go. If you decide not to use `CMake`, you need to
+manually link your application/library to `tug`.
+
+1. Create project directory.
+
+```bash
+mkdir sample_project && cd sample_project
+```
+
+2. Clone this repository into path of choice project directory
+
+3. Add the following line into `CMakeLists.txt` file:
+
+```bash
+add_subdirectory(path_to_tug EXCLUDE_FROM_ALL)
+```
+
+4. Write application/library using `tug`\'s API, notably
+    including relevant headers (see examples).
+
+5. Link target application/library against `tug`. Do this by
+    adding into according `CMakeLists.txt` file:
+
+```bash
+target_link_libraries(<your_target> tug)
+```
+
+6. Build your application/library with `CMake`.
+
+# Usage in an application
+
+Using `tug` can be summarized into the following steps:
+
+1. Define problem dimensionality
+2. Set grid sizes for each dimension
+3. Set the timestep
+4. Define boundary conditions
+5. Run the simulation!
+
+This will run a simulation on the defined grid for one species. See the
+source code documentation of `tug` and the examples in the
+`examples/` directory for more details.
+
+# Contributing
+
+In this early stage of development every help is welcome. To do so,
+there are currently the following options:
+
+Given you have an account for GFZ\'s `gitlab` instance:
+
+1. Fork this project, create a branch and push your changes. If your
+    changes are done or you feel the need for some feedback create a
+    merge request with the destination set to the **main** branch of
+    this project.
+2. Ask for access to this repository. You most likely will get access
+    as a developer which allows you to create branches and merge
+    requests inside this repository.
+
+If you can\'t get access to this `gitlab` instance:
+
+- Download this repository and note down the SHA of the downloaded commit. Apply
+    your changes and send a mail to <mluebke@gfz-potsdam.de> or
+    <delucia@gfz-potsdam.de> with the patch/diff compared to your starting
+    point. Please split different patch types (feature, fixes, improvements ...)
+    into seperate files. Also provide us the SHA of the commit you\'ve
+    downloaded.
+
+Thank you for your contributions in advance!

README.org

@@ -1,123 +0,0 @@
-#+TITLE: TUG: a C++ framework to solve Transport on Uniform Grids
-
-[[./doc/images/tug_logo_small.png]]
-
-=tug= implements different numerical approaches for transport
-problems, notably diffusion with implicit BTCS (Backward Time, Central
-Space) Euler and parallel 2D ADI (Alternating Direction Implicit).
-
-* About
-
-This project aims to provide a library for solving transport
-problems - diffusion, advection - on uniform grids implemented in C++.
-The library is built on top of [[https://eigen.tuxfamily.org/index.php?title=Main_Page][Eigen]], providing easy access to its
-optimized data structures and linear equation solvers.
-
-We designed the API to be as flexible as possible. Nearly every
-built-in, framework or third-party data structure can be used to model
-a problem, as long a pointer to continuous memory can be provided. We
-also provide parallelization using [[https://www.openmp.org/][OpenMP]], which can be easily turned
-on/off at compile time.
-
-At the current state, both 1D and 2D diffusion problems on a regular
-grid with constant alpha for all grid cells can be solved reliably.
-
-* Getting started
-
-=tug= is designed as a framework library and it relies on [[https://cmake.org/][CMake]] for
-building. If you already use =CMake= as your build toolkit for your
-application, you're good to go. If you decide not to use =CMake=, you
-need to manually link your application/library to =tug=.
-
-1. Create project directory.
-
-   #+BEGIN_SRC
-   $ mkdir sample_project && cd sample_project
-   #+END_SRC
-
-2. Clone this repository into path of choice project directory
-   - with =ssh=:
-
-   #+BEGIN_SRC
-   $ git clone git@git.gfz-potsdam.de:sec34/tug.git
-   #+END_SRC
-
-   - with =https=: 
-   #+BEGIN_SRC
-   $ git clone https://git.gfz-potsdam.de/sec34/tug.git
-   #+END_SRC
-
-3. Add the following line into =CMakeLists.txt= file:
-
-   #+BEGIN_SRC
-   add_subdirectory(path_to_tug EXCLUDE_FROM_ALL)
-   #+END_SRC
-
-4. Write application/library using =tug='s API, notably including
-   relevant headers (see examples).
-
-5. Link target application/library against =tug=. Do this by adding
-   into according =CMakeLists.txt= file:
-
-   #+BEGIN_SRC
-   target_link_libraries(<your_target> tug)
-   #+END_SRC
-
-6. Build your application/library with =CMake=.
-
-
-* Usage in an application
-
-Using =tug= can be summarized into the following steps:
-
-1. Define problem dimensionality
-2. Set grid sizes for each dimension
-3. Set the timestep 
-4. Define boundary conditions
-5. Run the simulation!
-
-This will run a simulation on the defined grid for one species. See
-the source code documentation of =tug= and the examples in the
-=examples/= directory for more details.
-
-* Roadmap
-
-- [X] 1D diffusion using BTCS
-- [X] 2D diffusion with ADI
-- [ ] 3D diffusion (?)
-- [X] R-API (see [[https://git.gfz-potsdam.de/sec34/rcppbtcs][RcppBTCS]])
-- [-] Python-API (?)
-- [X] Testing
-
-* Contributing
-** *PLEASE NOTE*
-
-Starting from version v0.2 we would like to use more meaningful commit
-messages. An overview of good practices and conventions can be found
-[[https://www.conventionalcommits.org/en/v1.0.0/][here]].
-
-** Workflow
-
-In this early stage of development every help is welcome. To do so,
-there are currently the following options:
-
-Given you have an account for GFZ's =gitlab= instance:
-
-1. Fork this project, create a branch and push your changes. If your
-   changes are done or you feel the need for some feedback create a
-   merge request with the destination set to the *main* branch of this
-   project.
-2. Ask for access to this repository. You most likely will get access
-   as a developer which allows you to create branches and merge
-   requests inside this repository.
-
-If you can't get access to this =gitlab= instance:
-
-3. Download this repository and note down the SHA of the downloaded
-   commit. Apply your changes and send a mail to
-   [[mailto:mluebke@gfz-potsdam.de][mluebke@gfz-potsdam.de]] or [[mailto:delucia@gfz-potsdam.de][delucia@gfz-potsdam.de]] with the
-   patch/diff compared to your starting point. Please split different
-   patch types (feature, fixes, improvements ...) into seperate files.
-   Also provide us the SHA of the commit you've downloaded.
-
-Thank you for your contributions in advance!

docs_sphinx/Diffusion.rst

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

docs_sphinx/Grid.rst

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

docs_sphinx/conf.py

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

docs_sphinx/user.rst

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

examples/BTCS_1D_proto_example.cpp

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

examples/BTCS_2D_proto_example.cpp

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

examples/CMakeLists.txt

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

examples/CRNI_2D_proto_example.cpp

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

examples/FTCS_1D_proto_example.cpp

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

examples/FTCS_2D_proto_closed_mdl.cpp

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

examples/FTCS_2D_proto_example.cpp

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

examples/FTCS_2D_proto_example_mdl.cpp

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

examples/_deps/doctest-src

@@ -1 +0,0 @@
-Subproject commit b7c21ec5ceeadb4951b00396fc1e4642dd347e5f

examples/profiling_openmp.cpp

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

examples/profiling_speedup.cpp

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

examples/reference-FTCS_2D_closed.cpp

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

include/tug/Boundary.hpp

@@ -7,10 +7,15 @@
 #ifndef BOUNDARY_H_
 #define BOUNDARY_H_
 
-#include "Grid.hpp"
+#include "tug/Core/TugUtils.hpp"
+
+#include <Eigen/Dense>
 #include <cstddef>
 #include <cstdint>
-#include <exception>
+#include <map>
+#include <stdexcept>
+#include <utility>
+#include <vector>
 
 namespace tug {
 
@@ -41,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.
@@ -107,14 +112,20 @@ private:
 template <class T> class Boundary {
 public:
   /**
-   * @brief Creates a boundary object based on the passed grid object and
-   *        initializes the boundaries as closed.
+   * @brief Creates a boundary object for a 1D grid
    *
-   * @param grid Grid object on the basis of which the simulation takes place
-   *             and from which the dimensions (in 2D case) are taken.
+   * @param length Length of the grid
    */
-  Boundary(const Grid<T> &grid)
-      : dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) {
+  Boundary(std::uint32_t length) : Boundary(1, length) {};
+
+  /**
+   * @brief Creates a boundary object for a 2D grid
+   *
+   * @param rows Number of rows of the grid
+   * @param cols Number of columns of the grid
+   */
+  Boundary(std::uint32_t rows, std::uint32_t cols)
+      : 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
@@ -133,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.
@@ -142,13 +182,11 @@ public:
    * @param side Side to be set to closed, e.g. BC_SIDE_LEFT.
    */
   void setBoundarySideClosed(BC_SIDE side) {
-    if (this->dim == 1) {
-      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
-        throw std::invalid_argument(
-            "For the one-dimensional case, only the BC_SIDE_LEFT and "
-            "BC_SIDE_RIGHT borders exist.");
-      }
-    }
+    tug_assert((this->dim > 1) ||
+                   ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
+               "For the "
+               "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
+               "borders exist.");
 
     const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
     const int n = is_vertical ? this->rows : this->cols;
@@ -167,13 +205,11 @@ public:
    * page.
    */
   void setBoundarySideConstant(BC_SIDE side, double value) {
-    if (this->dim == 1) {
-      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
-        throw std::invalid_argument(
-            "For the one-dimensional case, only the BC_SIDE_LEFT and "
-            "BC_SIDE_RIGHT borders exist.");
-      }
-    }
+    tug_assert((this->dim > 1) ||
+                   ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
+               "For the "
+               "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
+               "borders exist.");
 
     const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT;
     const int n = is_vertical ? this->rows : this->cols;
@@ -193,10 +229,9 @@ public:
    */
   void setBoundaryElemenClosed(BC_SIDE side, int index) {
     // tests whether the index really points to an element of the boundary side.
-    if ((boundaries[side].size() < index) || index < 0) {
-      throw std::invalid_argument(
-          "Index is selected either too large or too small.");
-    }
+    tug_assert(boundaries[side].size() > index && index >= 0,
+               "Index is selected either too large or too small.");
+
     this->boundaries[side][index].setType(BC_TYPE_CLOSED);
   }
 
@@ -214,10 +249,8 @@ public:
    */
   void setBoundaryElementConstant(BC_SIDE side, int index, double value) {
     // tests whether the index really points to an element of the boundary side.
-    if ((boundaries[side].size() < index) || index < 0) {
-      throw std::invalid_argument(
-          "Index is selected either too large or too small.");
-    }
+    tug_assert(boundaries[side].size() > index && index >= 0,
+               "Index is selected either too large or too small.");
     this->boundaries[side][index].setType(BC_TYPE_CONSTANT);
     this->boundaries[side][index].setValue(value);
   }
@@ -232,13 +265,11 @@ public:
    * BoundaryElement<T> objects.
    */
   const std::vector<BoundaryElement<T>> &getBoundarySide(BC_SIDE side) const {
-    if (this->dim == 1) {
-      if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) {
-        throw std::invalid_argument(
-            "For the one-dimensional trap, only the BC_SIDE_LEFT and "
-            "BC_SIDE_RIGHT borders exist.");
-      }
-    }
+    tug_assert((this->dim > 1) ||
+                   ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)),
+               "For the "
+               "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT "
+               "borders exist.");
     return this->boundaries[side];
   }
 
@@ -276,10 +307,8 @@ public:
    * object.
    */
   BoundaryElement<T> getBoundaryElement(BC_SIDE side, int index) const {
-    if ((boundaries[side].size() < index) || index < 0) {
-      throw std::invalid_argument(
-          "Index is selected either too large or too small.");
-    }
+    tug_assert(boundaries[side].size() > index && index >= 0,
+               "Index is selected either too large or too small.");
     return this->boundaries[side][index];
   }
 
@@ -294,10 +323,8 @@ public:
    * @return Boundary Type of the corresponding boundary condition.
    */
   BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const {
-    if ((boundaries[side].size() < index) || index < 0) {
-      throw std::invalid_argument(
-          "Index is selected either too large or too small.");
-    }
+    tug_assert(boundaries[side].size() > index && index >= 0,
+               "Index is selected either too large or too small.");
     return this->boundaries[side][index].getType();
   }
 
@@ -314,17 +341,187 @@ public:
    * object.
    */
   T getBoundaryElementValue(BC_SIDE side, int index) const {
-    if ((boundaries[side].size() < index) || index < 0) {
-      throw std::invalid_argument(
-          "Index is selected either too large or too small.");
-    }
-    if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) {
-      throw std::invalid_argument(
-          "A value can only be output if it is a constant boundary condition.");
-    }
+    tug_assert(boundaries[side].size() > index && index >= 0,
+               "Index is selected either too large or too small.");
+    tug_assert(
+        boundaries[side][index].getType() == BC_TYPE_CONSTANT,
+        "A value can only be output if it is a constant boundary condition.");
+
     return this->boundaries[side][index].getValue();
   }
 
+  /**
+   * @brief Serializes the boundary conditions into a vector of BoundaryElement
+   * objects.
+   *
+   * @return Vector with BoundaryElement objects.
+   */
+  std::vector<BoundaryElement<T>> serialize() const {
+    std::vector<BoundaryElement<T>> serialized;
+    for (std::size_t side = 0; side < boundaries.size(); side++) {
+      for (std::size_t index = 0; index < boundaries[side].size(); index++) {
+        serialized.push_back(boundaries[side][index]);
+      }
+    }
+    return serialized;
+  }
+
+  /**
+   * @brief Deserializes the boundary conditions from a vector of
+   * BoundaryElement objects.
+   *
+   * @param serialized Vector with BoundaryElement objects.
+   */
+  void deserialize(const std::vector<BoundaryElement<T>> &serialized) {
+    std::size_t index = 0;
+    for (std::size_t side = 0; side < boundaries.size(); side++) {
+      for (std::size_t i = 0; i < boundaries[side].size(); i++) {
+        boundaries[side][i] = serialized[index];
+        index++;
+      }
+    }
+  }
+
+  /**
+   *
+   * @param index Index of the inner constant boundary condition
+   * @param value Value of the inner constant boundary condition
+   */
+  void setInnerBoundary(std::uint32_t index, T value) {
+    tug_assert(this->dim == 1, "This function is only available for 1D grids.");
+    tug_assert(index < this->cols, "Index is out of bounds.");
+
+    this->inner_boundary[std::make_pair(0, index)] = value;
+  }
+
+  /**
+   * @brief Set inner constant boundary condition in 2D case.
+   *
+   * @param row Row index of the inner constant boundary condition
+   * @param col Column index of the inner constant boundary condition
+   * @param value Value of the inner constant boundary condition
+   */
+  void setInnerBoundary(std::uint32_t row, std::uint32_t col, T value) {
+    tug_assert(this->dim == 2, "This function is only available for 2D grids.");
+    tug_assert(row < this->rows && col < this->cols, "Index is out of bounds.");
+
+    this->inner_boundary[std::make_pair(row, col)] = value;
+  }
+
+  /**
+   * @brief Get inner constant boundary condition in 1D case.
+   *
+   * @param index Index of the inner constant boundary condition
+   * @return std::pair<bool, T> Pair of boolean (whether constant boundary was
+   * set or not) and value of the inner constant boundary condition
+   */
+  std::pair<bool, T> getInnerBoundary(std::uint32_t index) const {
+    tug_assert(this->dim == 1, "This function is only available for 1D grids.");
+    tug_assert(index < this->cols, "Index is out of bounds.");
+
+    auto it = this->inner_boundary.find(std::make_pair(0, index));
+    if (it == this->inner_boundary.end()) {
+      return std::make_pair(false, -1);
+    }
+    return std::make_pair(true, it->second);
+  }
+
+  /**
+   * @brief Get inner constant boundary condition in 2D case.
+   *
+   * @param row Row index of the inner constant boundary condition
+   * @param col Column index of the inner constant boundary condition
+   * @return std::pair<bool, T> Pair of boolean (whether constant boundary was
+   * set or not) and value of the inner constant boundary condition
+   */
+  std::pair<bool, T> getInnerBoundary(std::uint32_t row,
+                                      std::uint32_t col) const {
+    tug_assert(this->dim == 2, "This function is only available for 2D grids.");
+    tug_assert(row < this->rows && col < this->cols, "Index is out of bounds.");
+
+    auto it = this->inner_boundary.find(std::make_pair(row, col));
+    if (it == this->inner_boundary.end()) {
+      return std::make_pair(false, -1);
+    }
+    return std::make_pair(true, it->second);
+  }
+
+  /**
+   * @brief Get inner constant boundary conditions of a row as a vector. Can be
+   * used for 1D grids (row == 0) or 2D grids.
+   *
+   * @param row Index of the row for which the inner boundary conditions are to
+   * be returned.
+   * @return std::vector<std::pair<bool, T>> Vector of pairs of boolean (whether
+   * constant boundary was set or not) and value of the inner constant boundary
+   * condition
+   */
+  std::vector<std::pair<bool, T>> getInnerBoundaryRow(std::uint32_t row) const {
+    tug_assert(row < this->rows, "Index is out of bounds.");
+
+    if (inner_boundary.empty()) {
+      return std::vector<std::pair<bool, T>>(this->cols,
+                                             std::make_pair(false, -1));
+    }
+
+    std::vector<std::pair<bool, T>> row_values;
+    for (std::uint32_t col = 0; col < this->cols; col++) {
+      row_values.push_back(getInnerBoundary(row, col));
+    }
+
+    return row_values;
+  }
+
+  /**
+   * @brief Get inner constant boundary conditions of a column as a vector. Can
+   * only be used for 2D grids.
+   *
+   * @param col Index of the column for which the inner boundary conditions are
+   * to be returned.
+   * @return std::vector<std::pair<bool, T>> Vector of pairs of boolean (whether
+   * constant boundary was set or not) and value of the inner constant boundary
+   * condition
+   */
+  std::vector<std::pair<bool, T>> getInnerBoundaryCol(std::uint32_t col) const {
+    tug_assert(this->dim == 2, "This function is only available for 2D grids.");
+    tug_assert(col < this->cols, "Index is out of bounds.");
+
+    if (inner_boundary.empty()) {
+      return std::vector<std::pair<bool, T>>(this->rows,
+                                             std::make_pair(false, -1));
+    }
+
+    std::vector<std::pair<bool, T>> col_values;
+    for (std::uint32_t row = 0; row < this->rows; row++) {
+      col_values.push_back(getInnerBoundary(row, col));
+    }
+
+    return col_values;
+  }
+
+  /**
+   * @brief Get inner constant boundary conditions as a map. Can be read by
+   * setInnerBoundaries.
+   *
+   * @return Map of inner constant boundary conditions
+   */
+  std::map<std::pair<std::uint32_t, std::uint32_t>, T>
+  getInnerBoundaries() const {
+    return this->inner_boundary;
+  }
+
+  /**
+   * @brief Set inner constant boundary conditions as a map. Can be obtained by
+   * getInnerBoundaries.
+   *
+   * @param inner_boundary Map of inner constant boundary conditions
+   */
+  void
+  setInnerBoundaries(const std::map<std::pair<std::uint32_t, std::uint32_t>, T>
+                         &inner_boundary) {
+    this->inner_boundary = inner_boundary;
+  }
+
 private:
   const std::uint8_t dim;
   const std::uint32_t cols;
@@ -332,6 +529,9 @@ private:
 
   std::vector<std::vector<BoundaryElement<T>>>
       boundaries; // Vector with Boundary Element information
+
+  // Inner boundary conditions for 1D and 2D grids identified by (row,col)
+  std::map<std::pair<std::uint32_t, std::uint32_t>, T> inner_boundary;
 };
 } // namespace tug
 #endif // BOUNDARY_H_

include/tug/Core/BaseSimulation.hpp

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

include/tug/Core/FTCS.hpp

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

include/tug/Core/Matrix.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <Eigen/Core>
+
+namespace tug {
+/**
+ * @brief Alias template for a row-major matrix using Eigen library.
+ *
+ * This alias template defines a type `RowMajMat` which represents a row-major
+ * matrix using the Eigen library. It is a template that takes a type `T` as its
+ * template parameter. The matrix is dynamically sized with `Eigen::Dynamic` for
+ * both rows and columns. The matrix is stored in row-major order.
+ *
+ * @tparam T The type of the matrix elements.
+ */
+template <typename T>
+using RowMajMat =
+    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+template <typename T> using RowMajMatMap = Eigen::Map<RowMajMat<T>>;
+} // namespace tug

include/tug/Core/Numeric/BTCS.hpp

@@ -10,12 +10,15 @@
 #ifndef BTCS_H_
 #define BTCS_H_
 
-#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>
@@ -25,6 +28,18 @@
 
 namespace tug {
 
+// optimization to remove Eigen sparse matrix
+template <class T> class Diagonals {
+public:
+  Diagonals() : left(), center(), right() {};
+  Diagonals(std::size_t size) : left(size), center(size), right(size) {};
+
+public:
+  std::vector<T> left;
+  std::vector<T> center;
+  std::vector<T> right;
+};
+
 // calculates coefficient for boundary in constant case
 template <class T>
 constexpr std::pair<T, T> calcBoundaryCoeffConstant(T alpha_center,
@@ -49,76 +64,84 @@ constexpr std::pair<T, T> calcBoundaryCoeffClosed(T alpha_center, T alpha_side,
 
 // creates coefficient matrix for next time step from alphas in x-direction
 template <class T>
-static Eigen::SparseMatrix<T>
-createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
+static Diagonals<T>
+createCoeffMatrix(const RowMajMat<T> &alpha,
                   const std::vector<BoundaryElement<T>> &bcLeft,
-                  const std::vector<BoundaryElement<T>> &bcRight, int numCols,
+                  const std::vector<BoundaryElement<T>> &bcRight,
+                  const std::vector<std::pair<bool, T>> &inner_bc, int numCols,
                   int rowIndex, T sx) {
 
   // square matrix of column^2 dimension for the coefficients
-  Eigen::SparseMatrix<T> cm(numCols, numCols);
-  cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
+  Diagonals<T> cm(numCols);
 
   // left column
-
-  switch (bcLeft[rowIndex].getType()) {
-  case BC_TYPE_CONSTANT: {
-    auto [centerCoeffTop, rightCoeffTop] =
-        calcBoundaryCoeffConstant(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
-    cm.insert(0, 0) = centerCoeffTop;
-    cm.insert(0, 1) = rightCoeffTop;
-    break;
-  }
-  case BC_TYPE_CLOSED: {
-    auto [centerCoeffTop, rightCoeffTop] =
-        calcBoundaryCoeffClosed(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
-    cm.insert(0, 0) = centerCoeffTop;
-    cm.insert(0, 1) = rightCoeffTop;
-    break;
-  }
-  default: {
-    throw_invalid_argument(
-        "Undefined Boundary Condition Type somewhere on Left or Top!");
-  }
+  if (inner_bc[0].first) {
+    cm.center[0] = 1;
+  } else {
+    switch (bcLeft[rowIndex].getType()) {
+    case BC_TYPE_CONSTANT: {
+      auto [centerCoeffTop, rightCoeffTop] =
+          calcBoundaryCoeffConstant(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
+      cm.center[0] = centerCoeffTop;
+      cm.right[0] = rightCoeffTop;
+      break;
+    }
+    case BC_TYPE_CLOSED: {
+      auto [centerCoeffTop, rightCoeffTop] =
+          calcBoundaryCoeffClosed(alpha(rowIndex, 0), alpha(rowIndex, 1), sx);
+      cm.center[0] = centerCoeffTop;
+      cm.right[0] = rightCoeffTop;
+      break;
+    }
+    default: {
+      throw_invalid_argument(
+          "Undefined Boundary Condition Type somewhere on Left or Top!");
+    }
+    }
   }
 
   // inner columns
   int n = numCols - 1;
   for (int i = 1; i < n; i++) {
-    cm.insert(i, i - 1) =
+    if (inner_bc[i].first) {
+      cm.center[i] = 1;
+      continue;
+    }
+    cm.left[i] =
         -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i));
-    cm.insert(i, i) =
+    cm.center[i] =
         1 +
         sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) +
               calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)));
-    cm.insert(i, i + 1) =
+    cm.right[i] =
         -sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1));
   }
 
   // right column
-
-  switch (bcRight[rowIndex].getType()) {
-  case BC_TYPE_CONSTANT: {
-    auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant(
-        alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
-    cm.insert(n, n - 1) = leftCoeffBottom;
-    cm.insert(n, n) = centerCoeffBottom;
-    break;
+  if (inner_bc[n].first) {
+    cm.center[n] = 1;
+  } else {
+    switch (bcRight[rowIndex].getType()) {
+    case BC_TYPE_CONSTANT: {
+      auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant(
+          alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
+      cm.left[n] = leftCoeffBottom;
+      cm.center[n] = centerCoeffBottom;
+      break;
+    }
+    case BC_TYPE_CLOSED: {
+      auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffClosed(
+          alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
+      cm.left[n] = leftCoeffBottom;
+      cm.center[n] = centerCoeffBottom;
+      break;
+    }
+    default: {
+      throw_invalid_argument(
+          "Undefined Boundary Condition Type somewhere on Right or Bottom!");
+    }
+    }
   }
-  case BC_TYPE_CLOSED: {
-    auto [centerCoeffBottom, leftCoeffBottom] =
-        calcBoundaryCoeffClosed(alpha(rowIndex, n), alpha(rowIndex, n - 1), sx);
-    cm.insert(n, n - 1) = leftCoeffBottom;
-    cm.insert(n, n) = centerCoeffBottom;
-    break;
-  }
-  default: {
-    throw_invalid_argument(
-        "Undefined Boundary Condition Type somewhere on Right or Bottom!");
-  }
-  }
-
-  cm.makeCompressed(); // important for Eigen solver
 
   return cm;
 }
@@ -146,15 +169,15 @@ 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 Eigen::MatrixX<T> &concentrations,
-                     const Eigen::MatrixX<T> &alphaX,
-                     const Eigen::MatrixX<T> &alphaY,
+createSolutionVector(const EigenType &concentrations,
+                     const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY,
                      const std::vector<BoundaryElement<T>> &bcLeft,
                      const std::vector<BoundaryElement<T>> &bcRight,
                      const std::vector<BoundaryElement<T>> &bcTop,
                      const std::vector<BoundaryElement<T>> &bcBottom,
+                     const std::vector<std::pair<bool, T>> &inner_bc,
                      int length, int rowIndex, T sx, T sy) {
 
   Eigen::VectorX<T> sv(length);
@@ -163,6 +186,10 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
   // inner rows
   if (rowIndex > 0 && rowIndex < numRows - 1) {
     for (int i = 0; i < length; i++) {
+      if (inner_bc[i].first) {
+        sv(i) = inner_bc[i].second;
+        continue;
+      }
       sv(i) =
           sy *
               calcAlphaIntercell(alphaY(rowIndex, i), alphaY(rowIndex + 1, i)) *
@@ -181,6 +208,10 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
   // first row
   else if (rowIndex == 0) {
     for (int i = 0; i < length; i++) {
+      if (inner_bc[i].first) {
+        sv(i) = inner_bc[i].second;
+        continue;
+      }
       switch (bcTop[i].getType()) {
       case BC_TYPE_CONSTANT: {
         sv(i) = calcExplicitConcentrationsBoundaryConstant(
@@ -204,6 +235,10 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
   // last row
   else if (rowIndex == numRows - 1) {
     for (int i = 0; i < length; i++) {
+      if (inner_bc[i].first) {
+        sv(i) = inner_bc[i].second;
+        continue;
+      }
       switch (bcBottom[i].getType()) {
       case BC_TYPE_CONSTANT: {
         sv(i) = calcExplicitConcentrationsBoundaryConstant(
@@ -226,13 +261,14 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
 
   // first column -> additional fixed concentration change from perpendicular
   // dimension in constant bc case
-  if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
+  if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT && !inner_bc[0].first) {
     sv(0) += 2 * sx * alphaX(rowIndex, 0) * bcLeft[rowIndex].getValue();
   }
 
   // last column -> additional fixed concentration change from perpendicular
   // dimension in constant bc case
-  if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
+  if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT &&
+      !inner_bc[length - 1].first) {
     sv(length - 1) +=
         2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue();
   }
@@ -244,12 +280,28 @@ createSolutionVector(const Eigen::MatrixX<T> &concentrations,
 // b to the solution vector
 // use of EigenLU solver
 template <class T>
-static Eigen::VectorX<T> EigenLUAlgorithm(Eigen::SparseMatrix<T> &A,
+static Eigen::VectorX<T> EigenLUAlgorithm(Diagonals<T> &A,
                                           Eigen::VectorX<T> &b) {
 
+  // convert A to Eigen sparse matrix
+  size_t dim_A = A.center.size();
+  Eigen::SparseMatrix<T> A_sparse(dim_A, dim_A);
+
+  A_sparse.insert(0, 0) = A.center[0];
+  A_sparse.insert(0, 1) = A.right[0];
+
+  for (size_t i = 1; i < dim_A - 1; i++) {
+    A_sparse.insert(i, i - 1) = A.left[i];
+    A_sparse.insert(i, i) = A.center[i];
+    A_sparse.insert(i, i + 1) = A.right[i];
+  }
+
+  A_sparse.insert(dim_A - 1, dim_A - 2) = A.left[dim_A - 1];
+  A_sparse.insert(dim_A - 1, dim_A - 1) = A.center[dim_A - 1];
+
   Eigen::SparseLU<Eigen::SparseMatrix<T>> solver;
-  solver.analyzePattern(A);
-  solver.factorize(A);
+  solver.analyzePattern(A_sparse);
+  solver.factorize(A_sparse);
 
   return solver.solve(b);
 }
@@ -258,27 +310,11 @@ static Eigen::VectorX<T> EigenLUAlgorithm(Eigen::SparseMatrix<T> &A,
 // b to the solution vector
 // implementation of Thomas Algorithm
 template <class T>
-static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
+static Eigen::VectorX<T> ThomasAlgorithm(Diagonals<T> &A,
                                          Eigen::VectorX<T> &b) {
   Eigen::Index n = b.size();
-
-  Eigen::VectorX<T> a_diag(n);
-  Eigen::VectorX<T> b_diag(n);
-  Eigen::VectorX<T> c_diag(n);
   Eigen::VectorX<T> x_vec = b;
 
-  // Fill diagonals vectors
-  b_diag[0] = A.coeff(0, 0);
-  c_diag[0] = A.coeff(0, 1);
-
-  for (Eigen::Index i = 1; i < n - 1; i++) {
-    a_diag[i] = A.coeff(i, i - 1);
-    b_diag[i] = A.coeff(i, i);
-    c_diag[i] = A.coeff(i, i + 1);
-  }
-  a_diag[n - 1] = A.coeff(n - 1, n - 2);
-  b_diag[n - 1] = A.coeff(n - 1, n - 1);
-
   // HACK: write CSV to file
 #ifdef WRITE_THOMAS_CSV
 #include <fstream>
@@ -304,20 +340,20 @@ static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
 
   // start solving - c_diag and x_vec are overwritten
   n--;
-  c_diag[0] /= b_diag[0];
-  x_vec[0] /= b_diag[0];
+  A.right[0] /= A.center[0];
+  x_vec[0] /= A.center[0];
 
   for (Eigen::Index i = 1; i < n; i++) {
-    c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1];
-    x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) /
-               (b_diag[i] - a_diag[i] * c_diag[i - 1]);
+    A.right[i] /= A.center[i] - A.left[i] * A.right[i - 1];
+    x_vec[i] = (x_vec[i] - A.left[i] * x_vec[i - 1]) /
+               (A.center[i] - A.left[i] * A.right[i - 1]);
   }
 
-  x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) /
-             (b_diag[n] - a_diag[n] * c_diag[n - 1]);
+  x_vec[n] = (x_vec[n] - A.left[n] * x_vec[n - 1]) /
+             (A.center[n] - A.left[n] * A.right[n - 1]);
 
   for (Eigen::Index i = n; i-- > 0;) {
-    x_vec[i] -= c_diag[i] * x_vec[i + 1];
+    x_vec[i] -= A.right[i] * x_vec[i + 1];
   }
 
   return x_vec;
@@ -325,130 +361,122 @@ 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,
-                    Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
+static void BTCS_1D(SimulationInput<T> &input,
+                    Eigen::VectorX<T> (*solverFunc)(Diagonals<T> &A,
                                                     Eigen::VectorX<T> &b)) {
-  int length = grid.getLength();
-  T sx = timestep / (grid.getDelta() * grid.getDelta());
+  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;
+  Diagonals<T> A;
   Eigen::VectorX<T> b(length);
 
-  const auto &alpha = grid.getAlpha();
+  const auto &alpha = input.alphaX;
+
+  const auto &bc = input.boundaries;
 
   const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
   const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
-  Eigen::MatrixX<T> concentrations = grid.getConcentrations();
+  const auto inner_bc = bc.getInnerBoundaryRow(0);
+
+  RowMajMatMap<T> &concentrations = input.concentrations;
   int rowIndex = 0;
-  A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex,
+  A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex,
                         sx); // this is exactly same as in 2D
   for (int i = 0; i < length; i++) {
     b(i) = concentrations(0, i);
   }
-  if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) {
+  if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT &&
+      !inner_bc[0].first) {
     b(0) += 2 * sx * alpha(0, 0) * bcLeft[0].getValue();
   }
-  if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) {
+  if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT &&
+      !inner_bc[length - 1].first) {
     b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue();
   }
 
-  concentrations_t1 = solverFunc(A, b);
-
-  for (int j = 0; j < length; j++) {
-    concentrations(0, j) = concentrations_t1(j);
-  }
-
-  grid.setConcentrations(concentrations);
+  concentrations = solverFunc(A, b);
 }
 
 // BTCS solution for 2D grid
 template <class T>
-static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
-                    Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
+static void BTCS_2D(SimulationInput<T> &input,
+                    Eigen::VectorX<T> (*solverFunc)(Diagonals<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);
 
-  Eigen::MatrixX<T> concentrations_t1 =
-      Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
-  Eigen::VectorX<T> row_t1(colMax);
+  RowMajMat<T> concentrations_t1(rowMax, colMax);
 
-  Eigen::SparseMatrix<T> A;
+  Diagonals<T> A;
   Eigen::VectorX<T> b;
 
-  auto alphaX = grid.getAlphaX();
-  auto 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);
 
-  Eigen::MatrixX<T> concentrations = grid.getConcentrations();
+  RowMajMatMap<T> &concentrations = input.concentrations;
 
-#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
+#pragma omp parallel for num_threads(numThreads) private(A, b)
   for (int i = 0; i < rowMax; i++) {
+    auto inner_bc = bc.getInnerBoundaryRow(i);
 
-    A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
+    A = createCoeffMatrix(alphaX, bcLeft, bcRight, inner_bc, colMax, i, sx);
     b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
-                             bcTop, bcBottom, colMax, i, sx, sy);
+                             bcTop, bcBottom, inner_bc, colMax, i, sx, sy);
 
-    row_t1 = solverFunc(A, b);
-
-    concentrations_t1.row(i) = row_t1;
+    concentrations_t1.row(i) = solverFunc(A, b);
   }
 
   concentrations_t1.transposeInPlace();
-  concentrations.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, row_t1)
+#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, rowMax, i, sy);
-    b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
-                             bcLeft, bcRight, 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);
 
-    row_t1 = solverFunc(A, b);
-
-    concentrations.row(i) = row_t1;
+    concentrations.col(i) = solverFunc(A, b);
   }
-
-  concentrations.transposeInPlace();
-
-  grid.setConcentrations(concentrations);
 }
 
 // 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, 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

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

include/tug/Core/Numeric/SimulationInput.hpp

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

include/tug/Core/TugUtils.hpp

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

include/tug/Diffusion.hpp

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

include/tug/Grid.hpp

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

naaice/BTCS_2D_NAAICE.cpp

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

naaice/NAAICE_dble_vs_float.cpp

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

test/CMakeLists.txt

@@ -1,18 +1,24 @@
-find_library(DOCTEST_LIB doctest)
+include(FetchContent)
 
-if(NOT DOCTEST_LIB)
-  include(FetchContent)
+FetchContent_Declare( 
+  googletest 
+  GIT_REPOSITORY https://github.com/google/googletest.git
+  GIT_TAG v1.15.2
+)
 
-  FetchContent_Declare(
-    DocTest
-    GIT_REPOSITORY https://github.com/doctest/doctest.git
-    GIT_TAG v2.4.9)
+FetchContent_MakeAvailable(googletest)
 
-  FetchContent_MakeAvailable(DocTest)
-endif()
 
-add_executable(testTug setup.cpp testSimulation.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp)
-target_link_libraries(testTug doctest tug)
+add_executable(testTug 
+setup.cpp
+testDiffusion.cpp 
+testFTCS.cpp 
+testBoundary.cpp
+)
+target_link_libraries(testTug tug GTest::gtest)
+  
+include(GoogleTest)
+gtest_discover_tests(testTug)
 
 # get relative path of the CSV file
 get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH)
@@ -20,8 +26,3 @@ get_filename_component(testSimulationCSV "FTCS_11_11_7000.csv" REALPATH)
 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")
-
-add_custom_target(
-  check
-  COMMAND $<TARGET_FILE:testTug>
-  DEPENDS testTug)

test/TestUtils.hpp

@@ -1,12 +1,15 @@
+#include "tug/Core/Matrix.hpp"
 #include <Eigen/Core>
 #include <Eigen/Dense>
 #include <fstream>
-#include <ios>
-#include <iostream>
 #include <sstream>
 #include <stdexcept>
 
-inline Eigen::MatrixXd CSV2Eigen(std::string file2Convert) {
+#include <gtest/gtest.h>
+
+#define TUG_TEST(x) TEST(Tug, x)
+
+inline tug::RowMajMat<double> CSV2Eigen(std::string file2Convert) {
 
   std::vector<double> matrixEntries;
 
@@ -29,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/setup.cpp

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

test/testBoundary.cpp

@@ -1,77 +1,102 @@
-#include <doctest/doctest.h>
-#include <iostream>
-#include <stdio.h>
-#include <string>
+#include "gtest/gtest.h"
+#include <stdexcept>
 #include <tug/Boundary.hpp>
-#include <typeinfo>
+#include <utility>
+#include <vector>
 
 using namespace std;
 using namespace tug;
 
-TEST_CASE("BoundaryElement") {
+#include <gtest/gtest.h>
 
-  SUBCASE("Closed case") {
-    BoundaryElement boundaryElementClosed = BoundaryElement<double>();
-    CHECK_NOTHROW(BoundaryElement<double>());
-    CHECK_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED);
-    CHECK_EQ(boundaryElementClosed.getValue(), -1);
-    CHECK_THROWS(boundaryElementClosed.setValue(0.2));
-  }
+#define BOUNDARY_TEST(x) TEST(Boundary, x)
 
-  SUBCASE("Constant case") {
-    BoundaryElement boundaryElementConstant = BoundaryElement(0.1);
-    CHECK_NOTHROW(BoundaryElement(0.1));
-    CHECK_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT);
-    CHECK_EQ(boundaryElementConstant.getValue(), 0.1);
-    CHECK_NOTHROW(boundaryElementConstant.setValue(0.2));
-    CHECK_EQ(boundaryElementConstant.getValue(), 0.2);
-  }
+BOUNDARY_TEST(Element) {
+
+  BoundaryElement boundaryElementClosed = BoundaryElement<double>();
+  EXPECT_NO_THROW(BoundaryElement<double>());
+  EXPECT_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED);
+  EXPECT_DOUBLE_EQ(boundaryElementClosed.getValue(), -1);
+  EXPECT_THROW(boundaryElementClosed.setValue(0.2), std::invalid_argument);
+
+  BoundaryElement boundaryElementConstant = BoundaryElement(0.1);
+  EXPECT_NO_THROW(BoundaryElement(0.1));
+  EXPECT_EQ(boundaryElementConstant.getType(), BC_TYPE_CONSTANT);
+  EXPECT_DOUBLE_EQ(boundaryElementConstant.getValue(), 0.1);
+  EXPECT_NO_THROW(boundaryElementConstant.setValue(0.2));
+  EXPECT_DOUBLE_EQ(boundaryElementConstant.getValue(), 0.2);
 }
 
-TEST_CASE("Boundary Class") {
-  Grid grid1D = Grid64(10);
-  Grid grid2D = Grid64(10, 12);
-  Boundary boundary1D = Boundary(grid1D);
-  Boundary boundary2D = Boundary(grid2D);
+BOUNDARY_TEST(Class) {
+  Boundary<double> boundary1D(10);
+  Boundary<double> boundary2D(10, 12);
   vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));
 
-  SUBCASE("Boundaries 1D case") {
-    CHECK_NOTHROW(Boundary boundary(grid1D));
-    CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1);
-    CHECK_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1);
-    CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
-             BC_TYPE_CLOSED);
-    CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_TOP));
-    CHECK_THROWS(boundary1D.getBoundarySide(BC_SIDE_BOTTOM));
-    CHECK_NOTHROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT));
-    CHECK_THROWS(boundary1D.setBoundarySideClosed(BC_SIDE_TOP));
-    CHECK_NOTHROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
-    CHECK_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
-    CHECK_THROWS(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2));
-    CHECK_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
-             BC_TYPE_CONSTANT);
-    CHECK_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
-             boundary1DVector[0].getType());
+  constexpr double inner_condition_value = -5;
+  constexpr std::pair<bool, double> innerBoundary =
+      std::make_pair(true, inner_condition_value);
+
+  std::vector<std::pair<bool, double>> row_ibc(12, std::make_pair(false, -1));
+  row_ibc[1] = innerBoundary;
+
+  std::vector<std::pair<bool, double>> col_ibc(10, std::make_pair(false, -1));
+  col_ibc[0] = innerBoundary;
+
+  {
+    EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_LEFT).size(), 1);
+    EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1);
+    EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
+              BC_TYPE_CLOSED);
+    EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_TOP),
+                 ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*");
+    EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_BOTTOM),
+                 ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*");
+    EXPECT_NO_THROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT));
+    EXPECT_DEATH(boundary1D.setBoundarySideClosed(BC_SIDE_TOP),
+                 ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*");
+    EXPECT_NO_THROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
+    EXPECT_DOUBLE_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
+    EXPECT_DEATH(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2),
+                 ".*Index is selected either too large or too small.*");
+    EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0),
+              BC_TYPE_CONSTANT);
+    EXPECT_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
+              boundary1DVector[0].getType());
+
+    EXPECT_NO_THROW(boundary1D.setInnerBoundary(0, inner_condition_value));
+    EXPECT_DEATH(boundary1D.setInnerBoundary(0, 0, inner_condition_value),
+                 ".*only available for 2D grids.*");
+    EXPECT_EQ(boundary1D.getInnerBoundary(0), innerBoundary);
+    EXPECT_FALSE(boundary1D.getInnerBoundary(1).first);
   }
 
-  SUBCASE("Boundaries 2D case") {
-    CHECK_NOTHROW(Boundary boundary(grid1D));
-    CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_LEFT).size(), 10);
-    CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_RIGHT).size(), 10);
-    CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_TOP).size(), 12);
-    CHECK_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12);
-    CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
-             BC_TYPE_CLOSED);
-    CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_TOP));
-    CHECK_NOTHROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM));
-    CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT));
-    CHECK_NOTHROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP));
-    CHECK_NOTHROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
-    CHECK_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
-    CHECK_THROWS(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12));
-    CHECK_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
-             BC_TYPE_CONSTANT);
-    CHECK_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
-             boundary1DVector[0].getType());
+  {
+    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);
+    EXPECT_EQ(boundary2D.getBoundarySide(BC_SIDE_BOTTOM).size(), 12);
+    EXPECT_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
+              BC_TYPE_CLOSED);
+    EXPECT_NO_THROW(boundary2D.getBoundarySide(BC_SIDE_TOP));
+    EXPECT_NO_THROW(boundary2D.getBoundarySide(BC_SIDE_BOTTOM));
+    EXPECT_NO_THROW(boundary2D.setBoundarySideClosed(BC_SIDE_LEFT));
+    EXPECT_NO_THROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP));
+    EXPECT_NO_THROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0));
+    EXPECT_DOUBLE_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0);
+    EXPECT_DEATH(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12),
+                 ".*too large or too small.*");
+    EXPECT_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0),
+              BC_TYPE_CONSTANT);
+    EXPECT_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(),
+              boundary1DVector[0].getType());
+
+    EXPECT_DEATH(boundary2D.setInnerBoundary(0, inner_condition_value),
+                 ".* 1D .*");
+    EXPECT_NO_THROW(boundary2D.setInnerBoundary(0, 1, inner_condition_value));
+    EXPECT_EQ(boundary2D.getInnerBoundary(0, 1), innerBoundary);
+    EXPECT_FALSE(boundary2D.getInnerBoundary(0, 2).first);
+
+    EXPECT_EQ(boundary2D.getInnerBoundaryRow(0), row_ibc);
+    EXPECT_EQ(boundary2D.getInnerBoundaryCol(1), col_ibc);
   }
 }

test/testDiffusion.cpp

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

test/testFTCS.cpp

@@ -1,20 +1,18 @@
+#include <gtest/gtest.h>
 #include <tug/Core/TugUtils.hpp>
 
-#include <doctest/doctest.h>
-#include <limits>
+#include <gtest/gtest.h>
 
-TEST_CASE("Maths") {
-  SUBCASE("mean between two alphas") {
-    double alpha1 = 10;
-    double alpha2 = 20;
-    double average = 15;
-    double harmonicMean =
-        double(2) / ((double(1) / alpha1) + (double(1) / alpha2));
+TEST(FTCS, calcAlphaIntercell) {
+  double alpha1 = 10;
+  double alpha2 = 20;
+  double average = 15;
+  double harmonicMean =
+      double(2) / ((double(1) / alpha1) + (double(1) / alpha2));
 
-    // double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) -
-    // harmonicMean); CHECK(difference <
-    // std::numeric_limits<double>::epsilon());
-    CHECK_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean);
-    CHECK_EQ(calcAlphaIntercell(alpha1, alpha2, false), average);
-  }
+  // double difference = std::fabs(calcAlphaIntercell(alpha1, alpha2) -
+  // harmonicMean); CHECK(difference <
+  // std::numeric_limits<double>::epsilon());
+  EXPECT_DOUBLE_EQ(calcAlphaIntercell(alpha1, alpha2), harmonicMean);
+  EXPECT_DOUBLE_EQ(calcAlphaIntercell(alpha1, alpha2, false), average);
 }

test/testGrid.cpp

@@ -1,253 +0,0 @@
-#include <Eigen/Core>
-#include <doctest/doctest.h>
-#include <tug/Grid.hpp>
-
-using namespace Eigen;
-using namespace std;
-using namespace tug;
-
-TEST_CASE("1D Grid, too small length") {
-  int l = 2;
-  CHECK_THROWS(Grid64(l));
-}
-
-TEST_CASE("2D Grid64, too small side") {
-  int r = 2;
-  int c = 4;
-  CHECK_THROWS(Grid64(r, c));
-
-  r = 4;
-  c = 2;
-  CHECK_THROWS(Grid64(r, c));
-}
-
-TEST_CASE("1D Grid64") {
-  int l = 12;
-  Grid64 grid(l);
-
-  SUBCASE("correct construction") {
-    CHECK_EQ(grid.getDim(), 1);
-    CHECK_EQ(grid.getLength(), l);
-    CHECK_EQ(grid.getCol(), l);
-    CHECK_EQ(grid.getRow(), 1);
-
-    CHECK_EQ(grid.getConcentrations().rows(), 1);
-    CHECK_EQ(grid.getConcentrations().cols(), l);
-    CHECK_EQ(grid.getAlpha().rows(), 1);
-    CHECK_EQ(grid.getAlpha().cols(), l);
-    CHECK_EQ(grid.getDeltaCol(), 1);
-
-    CHECK_THROWS(grid.getAlphaX());
-    CHECK_THROWS(grid.getAlphaY());
-    CHECK_THROWS(grid.getDeltaRow());
-  }
-
-  SUBCASE("setting concentrations") {
-    // correct concentrations matrix
-    MatrixXd concentrations = MatrixXd::Constant(1, l, 3);
-    CHECK_NOTHROW(grid.setConcentrations(concentrations));
-
-    // false concentrations matrix
-    MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4);
-    CHECK_THROWS(grid.setConcentrations(wConcentrations));
-  }
-
-  SUBCASE("setting alpha") {
-    // correct alpha matrix
-    MatrixXd alpha = MatrixXd::Constant(1, l, 3);
-    CHECK_NOTHROW(grid.setAlpha(alpha));
-
-    CHECK_THROWS(grid.setAlpha(alpha, alpha));
-
-    grid.setAlpha(alpha);
-    CHECK_EQ(grid.getAlpha(), alpha);
-    CHECK_THROWS(grid.getAlphaX());
-    CHECK_THROWS(grid.getAlphaY());
-
-    // false alpha matrix
-    MatrixXd wAlpha = MatrixXd::Constant(3, l, 2);
-    CHECK_THROWS(grid.setAlpha(wAlpha));
-  }
-
-  SUBCASE("setting domain") {
-    int d = 8;
-    // set 1D domain
-    CHECK_NOTHROW(grid.setDomain(d));
-
-    // set 2D domain
-    CHECK_THROWS(grid.setDomain(d, d));
-
-    grid.setDomain(d);
-    CHECK_EQ(grid.getDeltaCol(), double(d) / double(l));
-    CHECK_THROWS(grid.getDeltaRow());
-
-    // set too small domain
-    d = 0;
-    CHECK_THROWS(grid.setDomain(d));
-    d = -2;
-    CHECK_THROWS(grid.setDomain(d));
-  }
-}
-
-TEST_CASE("2D Grid64 quadratic") {
-  int rc = 12;
-  Grid64 grid(rc, rc);
-
-  SUBCASE("correct construction") {
-    CHECK_EQ(grid.getDim(), 2);
-    CHECK_THROWS(grid.getLength());
-    CHECK_EQ(grid.getCol(), rc);
-    CHECK_EQ(grid.getRow(), rc);
-
-    CHECK_EQ(grid.getConcentrations().rows(), rc);
-    CHECK_EQ(grid.getConcentrations().cols(), rc);
-    CHECK_THROWS(grid.getAlpha());
-
-    CHECK_EQ(grid.getAlphaX().rows(), rc);
-    CHECK_EQ(grid.getAlphaX().cols(), rc);
-    CHECK_EQ(grid.getAlphaY().rows(), rc);
-    CHECK_EQ(grid.getAlphaY().cols(), rc);
-    CHECK_EQ(grid.getDeltaRow(), 1);
-    CHECK_EQ(grid.getDeltaCol(), 1);
-  }
-
-  SUBCASE("setting concentrations") {
-    // correct concentrations matrix
-    MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2);
-    CHECK_NOTHROW(grid.setConcentrations(concentrations));
-
-    // false concentrations matrix
-    MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1);
-    CHECK_THROWS(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(rc + 3, rc, 4);
-    CHECK_THROWS(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6);
-    CHECK_THROWS(grid.setConcentrations(wConcentrations));
-  }
-
-  SUBCASE("setting alphas") {
-    // correct alpha matrices
-    MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
-    MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
-    CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
-
-    CHECK_THROWS(grid.setAlpha(alphax));
-
-    grid.setAlpha(alphax, alphay);
-    CHECK_EQ(grid.getAlphaX(), alphax);
-    CHECK_EQ(grid.getAlphaY(), alphay);
-    CHECK_THROWS(grid.getAlpha());
-
-    // false alpha matrices
-    alphax = MatrixXd::Constant(rc + 3, rc + 1, 3);
-    CHECK_THROWS(grid.setAlpha(alphax, alphay));
-    alphay = MatrixXd::Constant(rc + 2, rc + 1, 3);
-    CHECK_THROWS(grid.setAlpha(alphax, alphay));
-  }
-
-  SUBCASE("setting domain") {
-    int dr = 8;
-    int dc = 9;
-
-    // set 1D domain
-    CHECK_THROWS(grid.setDomain(dr));
-
-    // set 2D domain
-    CHECK_NOTHROW(grid.setDomain(dr, dc));
-
-    grid.setDomain(dr, dc);
-    CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc));
-    CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc));
-
-    // set too small domain
-    dr = 0;
-    CHECK_THROWS(grid.setDomain(dr, dc));
-    dr = 8;
-    dc = 0;
-    CHECK_THROWS(grid.setDomain(dr, dc));
-    dr = -2;
-    CHECK_THROWS(grid.setDomain(dr, dc));
-  }
-}
-
-TEST_CASE("2D Grid64 non-quadratic") {
-  int r = 12;
-  int c = 15;
-  Grid64 grid(r, c);
-
-  SUBCASE("correct construction") {
-    CHECK_EQ(grid.getDim(), 2);
-    CHECK_THROWS(grid.getLength());
-    CHECK_EQ(grid.getCol(), c);
-    CHECK_EQ(grid.getRow(), r);
-
-    CHECK_EQ(grid.getConcentrations().rows(), r);
-    CHECK_EQ(grid.getConcentrations().cols(), c);
-    CHECK_THROWS(grid.getAlpha());
-
-    CHECK_EQ(grid.getAlphaX().rows(), r);
-    CHECK_EQ(grid.getAlphaX().cols(), c);
-    CHECK_EQ(grid.getAlphaY().rows(), r);
-    CHECK_EQ(grid.getAlphaY().cols(), c);
-    CHECK_EQ(grid.getDeltaRow(), 1);
-    CHECK_EQ(grid.getDeltaCol(), 1);
-  }
-
-  SUBCASE("setting concentrations") {
-    // correct concentrations matrix
-    MatrixXd concentrations = MatrixXd::Constant(r, c, 2);
-    CHECK_NOTHROW(grid.setConcentrations(concentrations));
-
-    // false concentrations matrix
-    MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6);
-    CHECK_THROWS(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(r + 3, c, 3);
-    CHECK_THROWS(grid.setConcentrations(wConcentrations));
-    wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2);
-    CHECK_THROWS(grid.setConcentrations(wConcentrations));
-  }
-
-  SUBCASE("setting alphas") {
-    // correct alpha matrices
-    MatrixXd alphax = MatrixXd::Constant(r, c, 2);
-    MatrixXd alphay = MatrixXd::Constant(r, c, 4);
-    CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
-
-    CHECK_THROWS(grid.setAlpha(alphax));
-
-    grid.setAlpha(alphax, alphay);
-    CHECK_EQ(grid.getAlphaX(), alphax);
-    CHECK_EQ(grid.getAlphaY(), alphay);
-    CHECK_THROWS(grid.getAlpha());
-
-    // false alpha matrices
-    alphax = MatrixXd::Constant(r + 3, c + 1, 3);
-    CHECK_THROWS(grid.setAlpha(alphax, alphay));
-    alphay = MatrixXd::Constant(r + 2, c + 1, 5);
-    CHECK_THROWS(grid.setAlpha(alphax, alphay));
-  }
-
-  SUBCASE("setting domain") {
-    int dr = 8;
-    int dc = 9;
-
-    // set 1D domain
-    CHECK_THROWS(grid.setDomain(dr));
-
-    // set 2D domain
-    CHECK_NOTHROW(grid.setDomain(dr, dc));
-
-    grid.setDomain(dr, dc);
-    CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
-    CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
-
-    // set too small domain
-    dr = 0;
-    CHECK_THROWS(grid.setDomain(dr, dc));
-    dr = 8;
-    dc = -1;
-    CHECK_THROWS(grid.setDomain(dr, dc));
-    dr = -2;
-    CHECK_THROWS(grid.setDomain(dr, dc));
-  }
-}

test/testSimulation.cpp

@@ -1,153 +0,0 @@
-#include "TestUtils.hpp"
-
-#include <Eigen/src/Core/Matrix.h>
-#include <doctest/doctest.h>
-#include <stdio.h>
-#include <string>
-#include <tug/Simulation.hpp>
-
-// include the configured header file
-#include <testSimulation.hpp>
-
-using namespace Eigen;
-using namespace std;
-using namespace tug;
-
-Grid64 setupSimulation(double timestep, int iterations) {
-  int row = 11;
-  int col = 11;
-  int domain_row = 10;
-  int domain_col = 10;
-
-  // Grid
-  Grid grid = Grid64(row, col);
-  grid.setDomain(domain_row, domain_col);
-
-  MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
-  concentrations(5, 5) = 1;
-  grid.setConcentrations(concentrations);
-
-  MatrixXd alpha = MatrixXd::Constant(row, col, 1);
-  for (int i = 0; i < 5; i++) {
-    for (int j = 0; j < 6; j++) {
-      alpha(i, j) = 0.01;
-    }
-  }
-  for (int i = 0; i < 5; i++) {
-    for (int j = 6; j < 11; j++) {
-      alpha(i, j) = 0.001;
-    }
-  }
-  for (int i = 5; i < 11; i++) {
-    for (int j = 6; j < 11; j++) {
-      alpha(i, j) = 0.1;
-    }
-  }
-  grid.setAlpha(alpha, alpha);
-
-  return grid;
-}
-
-constexpr double timestep = 0.001;
-constexpr double iterations = 7000;
-
-TEST_CASE("equality to reference matrix with FTCS") {
-  // set string from the header file
-  string test_path = testSimulationCSVDir;
-  MatrixXd reference = CSV2Eigen(test_path);
-  cout << "FTCS Test: " << endl;
-
-  Grid grid = setupSimulation(timestep, iterations); // Boundary
-  Boundary bc = Boundary(grid);
-
-  // Simulation
-
-  Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc);
-  // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
-  sim.setTimestep(timestep);
-  sim.setIterations(iterations);
-  sim.run();
-
-  cout << endl;
-  CHECK(checkSimilarity(reference, grid.getConcentrations(), 0.1) == true);
-}
-
-TEST_CASE("equality to reference matrix with BTCS") {
-  // set string from the header file
-  string test_path = testSimulationCSVDir;
-  MatrixXd reference = CSV2Eigen(test_path);
-  cout << "BTCS Test: " << endl;
-
-  Grid grid = setupSimulation(timestep, iterations); // Boundary
-  Boundary bc = Boundary(grid);
-
-  // Simulation
-  Simulation sim = Simulation<double, tug::FTCS_APPROACH>(grid, bc);
-  // sim.setOutputConsole(CONSOLE_OUTPUT_ON);
-  sim.setTimestep(timestep);
-  sim.setIterations(iterations);
-  sim.run();
-
-  cout << endl;
-  CHECK(checkSimilarityV2(reference, grid.getConcentrations(), 0.01) == true);
-}
-
-TEST_CASE("Initialize environment") {
-  int rc = 12;
-  Grid64 grid(rc, rc);
-  Boundary boundary(grid);
-
-  CHECK_NOTHROW(Simulation sim(grid, boundary));
-}
-
-TEST_CASE("Simulation environment") {
-  int rc = 12;
-  Grid64 grid(rc, rc);
-  Boundary boundary(grid);
-  Simulation<double, tug::FTCS_APPROACH> sim(grid, boundary);
-
-  SUBCASE("default paremeters") { CHECK_EQ(sim.getIterations(), -1); }
-
-  SUBCASE("set iterations") {
-    CHECK_NOTHROW(sim.setIterations(2000));
-    CHECK_EQ(sim.getIterations(), 2000);
-    CHECK_THROWS(sim.setIterations(-300));
-  }
-
-  SUBCASE("set timestep") {
-    CHECK_NOTHROW(sim.setTimestep(0.1));
-    CHECK_EQ(sim.getTimestep(), 0.1);
-    CHECK_THROWS(sim.setTimestep(-0.3));
-  }
-}
-
-TEST_CASE("Closed Boundaries - no change expected") {
-
-  constexpr std::uint32_t nrows = 5;
-  constexpr std::uint32_t ncols = 5;
-
-  tug::Grid64 grid(nrows, ncols);
-
-  auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0);
-  auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
-  auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
-
-  grid.setConcentrations(concentrations);
-  grid.setAlpha(alphax, alphay);
-
-  tug::Boundary bc(grid);
-  bc.setBoundarySideConstant(tug::BC_SIDE_LEFT, 1.0);
-  bc.setBoundarySideConstant(tug::BC_SIDE_RIGHT, 1.0);
-  bc.setBoundarySideConstant(tug::BC_SIDE_TOP, 1.0);
-  bc.setBoundarySideConstant(tug::BC_SIDE_BOTTOM, 1.0);
-
-  tug::Simulation<double> sim(grid, bc);
-  sim.setTimestep(1);
-  sim.setIterations(1);
-
-  MatrixXd input_values(concentrations);
-  sim.run();
-
-  CHECK(checkSimilarityV2(input_values, grid.getConcentrations(), 1E-12) ==
-        true);
-}

utils/ci.Dockerfile

@@ -1,4 +1,14 @@
 FROM alpine
-    MAINTAINER Max Luebke <mluebke@uni-potsdam.de>
+
+LABEL maintainer="Max Luebke <mluebke@uni-potsdam.de"
 
 RUN apk add --no-cache build-base openmp cmake git eigen-dev
+
+RUN git clone https://github.com/doctest/doctest.git /doctest \
+    && cd /doctest \
+    && mkdir build \
+    && cd build \
+    && cmake .. \
+    && make install \
+    && cd / \
+    && rm -rf /doctest