implementing FTCS

include/tug/Boundary.hpp

@@ -1,3 +1,6 @@
+#ifndef BOUNDARY_H_
+#define BOUNDARY_H_
+
 #include <Eigen/Core>
 #include "Grid.hpp"
 
@@ -31,7 +34,7 @@ class Boundary {
          * 
          * @return auto 
          */
-        auto getBoundaryConditionType();
+        BC_TYPE getBoundaryConditionType();
 
         /**
          * @brief Set the Boundary Condition Value object
@@ -55,3 +58,5 @@ class Boundary {
         BC_TYPE type;
         VectorXd left, right, top, bottom;
 };
+
+#endif

include/tug/Grid.hpp

@@ -32,7 +32,7 @@ class Grid {
          * 
          * @return auto 
          */
-        auto getConcentrations();
+        Matrix2d getConcentrations();
 
         /**
         * @brief Set the Alpha object
@@ -49,18 +49,34 @@ class Grid {
         */
         void setAlpha(Matrix2d alpha_x, Matrix2d alpha_y);
 
+        Matrix2d getAlphaX();
+
+        Matrix2d getAlphaY();
+
         int getDim();
 
         int getRow();
 
         int getCol();
 
+        void setDomain(int domain_col);
+
+        void setDomain(int domain_row, int domain_col);
+
+        double getDeltaCol();
+
+        double getDeltaRow();
+
 
     private:
 
         int dim;
-        int row;
         int col;
+        int row;
+        int domain_col;
+        int domain_row;
+        double delta_col;
+        double delta_row;
         Matrix2d concentrations;
         Matrix2d alpha_x;
         Matrix2d alpha_y;

src/Boundary.cpp

@@ -23,7 +23,7 @@ Boundary::Boundary(Grid &grid, BC_TYPE type) : grid(grid) {
     } 
 }
 
-auto Boundary::getBoundaryConditionType() {
+BC_TYPE Boundary::getBoundaryConditionType() {
     return this->type; 
 }
 

src/FTCS.cpp

@@ -1,5 +1,63 @@
 #include <tug/Boundary.hpp>
 
-auto FTCS(Grid &grid, Boundary &bc, double timestep) {
 
-}
+auto calc_alpha_intercell(double alpha1, double alpha2, bool useHarmonic = false) {
+    if (useHarmonic) {
+        return 2 / ((1/alpha1) + (1/alpha2));
+    } else {
+        return 0.5 * (alpha1 + alpha2);
+    }
+}
+
+auto FTCS_constant(Grid &grid, Boundary &bc, double timestep) {
+    int rowMax = grid.getRow();
+    int colMax = grid.getCol();
+    double deltaRow = grid.getDeltaRow();
+    double deltaCol = grid.getDeltaCol();
+
+    // Matrix with concentrations at time t+1
+    // TODO profiler / only use 2 matrices
+    Matrix2d concentrations_t1 = Matrix2d(rowMax, colMax);
+
+    // inner cells
+    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)) 
+                    * 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)) 
+                    * 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)
+                );
+        }
+    }
+
+    // boundary conditions
+    // left without corners / looping over rows
+
+    return 0;
+}
+
+auto FTCS_closed(Grid &grid, Boundary &bc, double timestep) {
+    return 0;
+}
+
+auto FTCS(Grid &grid, Boundary &bc, double timestep) {
+    if (bc.getBoundaryConditionType() == BC_TYPE_CONSTANT) {
+        int test = FTCS_constant(grid, bc, timestep);
+    } else if (bc.getBoundaryConditionType() == BC_TYPE_CLOSED) {
+        FTCS_closed(grid, bc, timestep);
+    }
+}

src/Grid.cpp

@@ -22,7 +22,7 @@ void Grid::setConcentrations(Matrix2d concentrations) {
     this->concentrations = concentrations;
 }
 
-auto Grid::getConcentrations() {
+Matrix2d Grid::getConcentrations() {
     return this->concentrations;
 }
 
@@ -35,6 +35,14 @@ void Grid::setAlpha(Matrix2d alpha_x, Matrix2d alpha_y) {
     this->alpha_y = alpha_y;
 }
 
+Matrix2d Grid::getAlphaX() {
+    return this->alpha_x;
+}
+
+Matrix2d Grid::getAlphaY() {
+    return this->alpha_y;
+}
+
 int Grid::getDim() {
     return dim;
 }
@@ -45,4 +53,25 @@ int Grid::getRow() {
 
 int Grid::getCol() {
     return col;
-}
+}
+
+void Grid::setDomain(int domain_col) {
+    this->domain_col = domain_col;
+    this->delta_col = this->domain_col/this->col;
+}
+
+void Grid::setDomain(int domain_row, int domain_col) {
+    this->domain_row = domain_row;
+    this->domain_col = domain_col;
+
+    this->domain_row = this->domain_row/this->row;
+    this->domain_col = this->domain_col/this->col;
+}
+
+double Grid::getDeltaCol() {
+    return this->delta_col;
+}
+
+double Grid::getDeltaRow() {
+    return this->delta_row;
+}

src/Simulation.cpp

@@ -5,6 +5,10 @@
 
 using namespace std;
 
+// auto FTCS(Grid &grid, Boundary &bc, double timestep) {
+
+// }
+
 Simulation::Simulation(Grid &grid, Boundary &bc, APPROACH approach) : grid(grid), bc(bc) {
     //probably to DEBUG assignment of grid and bc
     this->grid = grid;
@@ -41,6 +45,7 @@ auto Simulation::getIterations() {
     return this->iterations;
 }
 
+
 auto Simulation::run() {
     if (approach == FTCS_APPROACH) {
         for (int i = 0; i < iterations; i++) {
@@ -52,4 +57,4 @@ auto Simulation::run() {
             break;
         }
     }
-}
+}