implemented the 1D FTCS case and created an example prototype

examples/CMakeLists.txt

@@ -2,8 +2,10 @@ add_executable(first_example first_example.cpp)
 add_executable(second_example second_example.cpp)
 add_executable(boundary_example1D boundary_example1D.cpp)
 add_executable(FTCS_2D_proto_example FTCS_2D_proto_example.cpp)
+add_executable(FTCS_1D_proto_example FTCS_1D_proto_example.cpp)
 
 target_link_libraries(first_example tug)
 target_link_libraries(second_example tug)
 target_link_libraries(boundary_example1D tug)
-target_link_libraries(FTCS_2D_proto_example tug)
+target_link_libraries(FTCS_2D_proto_example tug)
+target_link_libraries(FTCS_1D_proto_example tug)

examples/FTCS_1D_proto_example.cpp

@@ -0,0 +1,44 @@
+#include <tug/Simulation.hpp>
+
+int main(int argc, char *argv[]) {
+    // **************
+    // **** GRID ****
+    // **************
+
+    // create a linear grid with 20 cells
+    int cells = 20;
+    Grid grid = Grid(cells);
+
+    MatrixXd concentrations = MatrixXd::Constant(1,20,20);
+    grid.setConcentrations(concentrations);
+
+
+    // ******************
+    // **** BOUNDARY ****
+    // ******************
+
+    // create a boundary with constant values
+    Boundary bc = Boundary(grid, BC_TYPE_CONSTANT);
+
+
+    // ************************
+    // **** SIMULATION ENV ****
+    // ************************
+
+    // set up a simulation environment
+    Simulation simulation = Simulation(grid, bc, FTCS_APPROACH); // 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_example.cpp

@@ -25,7 +25,7 @@ int main(int argc, char *argv[]) {
     // (optional) set the concentrations, e.g.:
     // MatrixXd concentrations = MatrixXd::Constant(20,20,1000); // #row,#col,value
     // grid.setConcentrations(concentrations);
-    MatrixXd concentrations = MatrixXd::Constant(20,20,0);
+    MatrixXd concentrations = MatrixXd::Constant(20,20,20);
     // concentrations(0,0) = 2000;
     grid.setConcentrations(concentrations);
 
@@ -45,12 +45,12 @@ int main(int argc, char *argv[]) {
     // (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);
+    // VectorXd bc_zero_values = VectorXd::Constant(20,0);
-    bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_LEFT, bc_zero_values);
-    bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_RIGHT, bc_zero_values);
-    VectorXd bc_front_values = VectorXd::Constant(20,2000);
+    // VectorXd bc_front_values = VectorXd::Constant(20,2000);
-    bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_TOP, bc_front_values);
-    bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
+    // bc.setBoundaryConditionValue(BC_SIDE_BOTTOM, bc_zero_values);
 
 
     // ************************
@@ -64,10 +64,10 @@ int main(int argc, char *argv[]) {
     simulation.setTimestep(0.1); // timestep
 
     // (optional) set the number of iterations
-    simulation.setIterations(1000);
+    simulation.setIterations(100);
 
     // (optional) set kind of output [CSV_OUTPUT_OFF (default), CSV_OUTPUT_ON, CSV_OUTPUT_VERBOSE]
-    simulation.setOutputCSV(CSV_OUTPUT_VERBOSE);
+    simulation.setOutputCSV(CSV_OUTPUT_OFF);
     
     // **** RUN SIMULATION ****
 

src/FTCS.cpp

@@ -13,59 +13,156 @@ double calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = fal
 }
 
 // IN PROGRESS
-double verticalChange(Grid grid, int row, int col) {
+double calcHorizontalChange(Grid grid, int row, int col) {
+
+    double result = 
+        calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+            * grid.getConcentrations()(row,col+1)
+        - (
+            calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
+            + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            )
+            * grid.getConcentrations()(row,col)
+        + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            * grid.getConcentrations()(row,col-1);
+
+    return result;
+}
+
+// IN PROGRESS
+double calcVerticalChange(Grid grid, int row, int col) {
     
     double result =    
-                calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+        calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
-                * grid.getConcentrations()(row+1,col)
+            * grid.getConcentrations()(row+1,col)
-                - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+        - (
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
+            calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
-                * grid.getConcentrations()(row,col)
+            + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
+            )
-                * grid.getConcentrations()(row-1,col);
+            * grid.getConcentrations()(row,col)
+        + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
+            * grid.getConcentrations()(row-1,col);
 
     return result;
 }
 
 // IN PROGRESS
-double horizontalChange(Grid grid, int row, int col) {
+double calcHorizontalChangeLeftBoundary(Grid grid, Boundary bc, int row, int col) {
+
+    double result = 
+        calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+            * grid.getConcentrations()(row,col+1)
+        - (
+            calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+            + 2 * grid.getAlphaX()(row,col)
+            ) 
+            * grid.getConcentrations()(row,col) 
+        + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row);
+
+    return result;
+}
+
+// IN PROGRESS
+double calcHorizontalChangeRightBoundary(Grid grid, Boundary bc, int row, int col) {
+
+    double result = 
+        2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+        - (
+            calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            + 2 * grid.getAlphaX()(row,col)
+            ) 
+            * grid.getConcentrations()(row,col)
+        + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            * grid.getConcentrations()(row,col-1);
+
+    return result;
+}
+
+// IN PROGRESS
+double calcVerticalChangeTopBoundary(Grid grid, Boundary bc, int row, int col) {
     
     double result = 
-                calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+        calc_alpha_intercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) 
-                * grid.getConcentrations()(row,col+1)
+            * grid.getConcentrations()(row+1,col)
-                - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
+        - (
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)))
+            calc_alpha_intercell(grid.getAlphaY()(row+1, col), grid.getAlphaY()(row, col)) 
-                * grid.getConcentrations()(row,col)
+            + 2 * grid.getAlphaY()(row, col)
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+            ) 
-                * grid.getConcentrations()(row,col-1);
+            * grid.getConcentrations()(row, col)
+        + 2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col);
 
     return result;
 }
 
 // IN PROGRESS
-MatrixXd FTCS_1D_constant(Grid grid, Boundary bc, double timestep) {
+double calcVerticalChangeBottomBoundary(Grid grid, Boundary bc, int row, int col) {
+
+    double result = 
+        2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+        - (
+            calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) 
+            + 2 * grid.getAlphaY()(row, col)
+            ) 
+            * grid.getConcentrations()(row, col)
+        + calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) 
+            * grid.getConcentrations()(row-1,col);
+
+    return result;
+}
+
+
+// IN PROGRESS
+MatrixXd FTCS_1D(Grid grid, Boundary bc, double timestep) {
     int colMax = grid.getCol();
     double deltaCol = grid.getDeltaCol();
 
     MatrixXd concentrations_t1 = MatrixXd::Constant(1, colMax, 0);
 
+    // only one row in 1D case
     int row = 0;
-    // for (int col = 1; col < colMax-1; col++) {
+
-    //     concentrations_t1 = grid.getConcentrations()(row,col)
+    // inner cells
-    //         + horizontal_term();
+    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
+    int col = 0;
+    concentrations_t1(row, col) = grid.getConcentrations()(row,col)
+            + timestep / (deltaCol*deltaCol) 
+                * (
+                    calcHorizontalChangeLeftBoundary(grid, bc, row, col)
+                )
+            ;
+
+
+    // right boundary
+    col = colMax-1;
+    concentrations_t1(row,col) = grid.getConcentrations()(row,col)
+            + timestep / (deltaCol*deltaCol) 
+                * (
+                    calcHorizontalChangeRightBoundary(grid, bc, row, col)
+                )
+            ;
+
+    return concentrations_t1;
 }
 
 // IN PROGRESS
-MatrixXd FTCS_2D_constant(Grid grid, Boundary bc, double timestep) {
+// MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
-    int rowMax = grid.getRow();
+//     int rowMax = grid.getRow();
-    int colMax = grid.getCol();
+//     int colMax = grid.getCol();
-    double deltaRow = grid.getDeltaRow();
+//     double deltaRow = grid.getDeltaRow();
-    double deltaCol = grid.getDeltaCol();
+//     double deltaCol = grid.getDeltaCol();
 
-}
+// }
 
-MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
+MatrixXd FTCS_2D(Grid grid, Boundary bc, double timestep) {
     int rowMax = grid.getRow();
     int colMax = grid.getCol();
     double deltaRow = grid.getDeltaRow();
@@ -76,29 +173,17 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
     MatrixXd concentrations_t1 = MatrixXd::Constant(rowMax, colMax, 0);
 
     // inner cells
-    // (should have 7 calls to current concentration)
+    // these do not depend on the boundary condition type
     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) 
                     * (
-                    calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                        calcVerticalChange(grid, row, col)
-                    * grid.getConcentrations()(row+1,col)
-                    - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
-                    + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
-                    * grid.getConcentrations()(row,col)
-                    + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
-                    * grid.getConcentrations()(row-1,col)
                     )
                 + timestep / (deltaCol*deltaCol) 
                     * (
-                    calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                        calcHorizontalChange(grid, row, col)
-                    * grid.getConcentrations()(row,col+1)
-                    - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
-                    + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col)))
-                    * grid.getConcentrations()(row,col)
-                    + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                    * grid.getConcentrations()(row,col-1)
                     )
                 ;
         }
@@ -106,184 +191,119 @@ MatrixXd FTCS_constant(Grid grid, Boundary bc, double timestep) {
 
     // boundary conditions
     // left without corners / looping over rows
-    // (should have 6 calls to current concentration)
     int col = 0;
     for (int row = 1; row < rowMax-1; row++) {
         concentrations_t1(row, col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
-                * (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                * (
-                * grid.getConcentrations()(row,col+1)
+                    calcHorizontalChangeLeftBoundary(grid, bc, row, col)
-                - (calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) 
+                )
-                + 2 * grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col) 
-                + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row))
             + timestep / (deltaRow*deltaRow)
-                * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                * (
-                * grid.getConcentrations()(row+1,col)
+                    calcVerticalChange(grid, row, col)
-                - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+                )
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
+            ;
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
-                * grid.getConcentrations()(row-1,col));
     }
 
     // right without corners / looping over columns
-    // (should have 6 calls to current concentration)
     col = colMax-1;
     for (int row = 1; row < rowMax-1; row++) {
         concentrations_t1(row,col) = grid.getConcentrations()(row,col)
             + timestep / (deltaCol*deltaCol) 
-                * (2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                * (
-                - (calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
+                    calcHorizontalChangeRightBoundary(grid, bc, row, col)
-                + 2 * grid.getAlphaX()(row,col)) 
+                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                * grid.getConcentrations()(row,col-1))
             + timestep / (deltaRow*deltaRow)
-                * (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
+                * (
-                * grid.getConcentrations()(row+1,col)
+                    calcVerticalChange(grid, row, col)
-                - (calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) 
+                )
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col)))
+            ;
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row-1,col), grid.getAlphaY()(row,col))
-                * grid.getConcentrations()(row-1,col));
     }
 
 
     // top without corners / looping over cols
-    // (should have 6 calls to current concentration)
     int row = 0;
     for (int col=1; col<colMax-1;col++){
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
-            + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) * grid.getConcentrations()(1,col)
+            + timestep / (deltaRow*deltaRow) 
-                - (calc_alpha_intercell(grid.getAlphaY()(1, col), grid.getAlphaY()(0, col)) + 2 * grid.getAlphaY()(0, col)) * grid.getConcentrations()(0, col)
+                * (
-                + 2 * grid.getAlphaY()(0, col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col))
+                    calcVerticalChangeTopBoundary(grid, bc, row, col)
-            + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col+1)
+                )
-                - (calc_alpha_intercell(grid.getAlphaX()(0, col+1), grid.getAlphaX()(0, col)) + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col))) * grid.getConcentrations()(0, col)
+            + timestep / (deltaCol*deltaCol) 
-                + calc_alpha_intercell(grid.getAlphaX()(0, col-1), grid.getAlphaX()(0, col)) * grid.getConcentrations()(0, col-1));
+                * (
+                    calcHorizontalChange(grid, row, col)
+                )
+            ;
     }
 
     // bottom without corners / looping over cols
-    // (should have 6 calls to current concentration)
     row = rowMax-1;
     for(int col=1; col<colMax-1;col++){
         concentrations_t1(row, col) = grid.getConcentrations()(row, col)
-            + timestep/(grid.getDeltaRow()*grid.getDeltaRow()) * (2 * grid.getAlphaY()(row, col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+            + timestep / (deltaRow*deltaRow) 
-                - (calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) + 2 * grid.getAlphaY()(row, col)) * grid.getConcentrations()(row, col)
+                * (
-                + calc_alpha_intercell(grid.getAlphaY()(row, col), grid.getAlphaY()(row-1, col)) * grid.getConcentrations()(row-1,col))
+                    calcVerticalChangeBottomBoundary(grid, bc, row, col)
-            + timestep/(grid.getDeltaCol()*grid.getDeltaCol()) * (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) * grid.getConcentrations()(row, col+1)
+                )
-                - (calc_alpha_intercell(grid.getAlphaX()(row, col+1), grid.getAlphaX()(row, col)) + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col))) * grid.getConcentrations()(row, col)
+            + timestep / (deltaCol*deltaCol) 
-                + calc_alpha_intercell(grid.getAlphaX()(row, col-1), grid.getAlphaX()(row, col)) * grid.getConcentrations()(row, col-1));
+                * (
+                    calcHorizontalChange(grid, row, col)
+                )
+            ;
     }
 
     // corner top left
-    // (should have 5 calls to current concentration)
     row = 0;
     col = 0;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
         + timestep/(deltaCol*deltaCol)
             * (
-                calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col+1)
+                calcHorizontalChangeLeftBoundary(grid, bc, row, col)
-                - ( 
-                    calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) * grid.getConcentrations()(row+1,col)
+                calcVerticalChangeTopBoundary(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col)
             )
         ;
 
     // corner top right
-    // (should have 5 calls to current concentration)
     row = 0;
     col = colMax-1;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
         + timestep/(deltaCol*deltaCol)
             * (
-                2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                calcHorizontalChangeRightBoundary(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                * grid.getConcentrations()(row,col-1)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col)) * grid.getConcentrations()(row+1,col)
+                calcVerticalChangeTopBoundary(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row+1,col), grid.getAlphaY()(row,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_TOP)(col)
             )
         ;
 
     // corner bottom left
-    // (should have 5 calls to current concentration)
     row = rowMax-1;
     col = 0;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col)
         + timestep/(deltaCol*deltaCol)
             * (
-                calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col)) * grid.getConcentrations()(row,col+1)
+                calcHorizontalChangeLeftBoundary(grid, bc, row, col)
-                - ( 
-                    calc_alpha_intercell(grid.getAlphaX()(row,col+1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + 2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_LEFT)(row)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+                calcVerticalChangeBottomBoundary(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                * grid.getConcentrations()(row-1,col)
             )
         ;
 
     // corner bottom right
-    // (should have 5 calls to current concentration)
     row = rowMax-1;
     col = colMax-1;
     concentrations_t1(row,col) = grid.getConcentrations()(row,col) 
         + timestep/(deltaCol*deltaCol)
             * (
-                2 * grid.getAlphaX()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_RIGHT)(row)
+                calcHorizontalChangeRightBoundary(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                    + 2 * grid.getAlphaX()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaX()(row,col-1), grid.getAlphaX()(row,col))
-                * grid.getConcentrations()(row,col-1)
             )
         + timestep/(deltaRow*deltaRow)
             * (
-                2 * grid.getAlphaY()(row,col) * bc.getBoundaryConditionValue(BC_SIDE_BOTTOM)(col)
+                calcVerticalChangeBottomBoundary(grid, bc, row, col)
-                - (
-                    calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                    + 2 * grid.getAlphaY()(row,col)
-                )
-                * grid.getConcentrations()(row,col)
-                + calc_alpha_intercell(grid.getAlphaY()(row,col), grid.getAlphaY()(row-1,col))
-                * grid.getConcentrations()(row-1,col)
             )
         ;
 
@@ -297,13 +317,29 @@ MatrixXd FTCS_closed(Grid grid, Boundary bc, double timestep) {
 }
 
 MatrixXd FTCS(Grid grid, Boundary bc, double timestep) {
-    switch (bc.getBoundaryConditionType()) {
+    // inner cells 
-        case BC_TYPE_CONSTANT:
+    // TODO only the boundary cells are different in constant and closed case
-            return FTCS_constant(grid, bc, timestep);
+
-        case BC_TYPE_CLOSED:
+    // if 1D:
-            return FTCS_closed(grid, bc, timestep);
+    //      do inner cells 
-        default:
+    //      do left boundary according to bc type
-            // TODO handle
+    //      do right boundary according to bc type
-            return MatrixXd();
+    // if 2D:
+    //      do inner cells
+    //      do left boundaries according to bc type
+    //      do right boundaries according to bc type
+    //      ...
+
+    if (grid.getDim() == 1) {
+        return FTCS_1D(grid, bc, timestep);
+    } else {
+        return FTCS_2D(grid, bc, timestep);
     }
+
+    // checking the boundary condition type first does not work
+    // if the boundary condition types change dynamically for a grid
+    // meaning:
+    // - boundary condition type needs to be checked for every single boundary cell
+    //      -> this check is last in order
+    // - 
 }

src/Simulation.cpp

@@ -84,7 +84,11 @@ void Simulation::run() {
         }
         printConcentrationsConsole();
         if (csv_output >= CSV_OUTPUT_ON) {
-            printConcentrationsCSV("test", true);
+            bool append = false;
+            if (csv_output == CSV_OUTPUT_VERBOSE) {
+                append = true;
+            }
+            printConcentrationsCSV("test", append);
         }
     } else if (approach == BTCS_APPROACH) {
         for (int i = 0; i < iterations; i++) {