implementation of closed 1D and 2D cases

examples/BTCS_2D_proto_example.cpp

@@ -35,10 +35,14 @@ int main(int argc, char *argv[]) {
 
     // create a boundary with constant values
     Boundary bc = Boundary(grid);
-    bc.setBoundarySideConstant(BC_SIDE_LEFT, 0);
+    bc.setBoundarySideClosed(BC_SIDE_LEFT);
-    bc.setBoundarySideConstant(BC_SIDE_RIGHT, 0);
+    bc.setBoundarySideClosed(BC_SIDE_RIGHT);
-    bc.setBoundarySideConstant(BC_SIDE_TOP, 0);
+    bc.setBoundarySideClosed(BC_SIDE_TOP);
-    bc.setBoundarySideConstant(BC_SIDE_BOTTOM, 0);
+    bc.setBoundarySideClosed(BC_SIDE_BOTTOM);
+    // 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.:

src/BTCSv2.cpp

@@ -4,16 +4,71 @@
 using namespace Eigen;
 
 
-static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int numCols, int rowIndex, double sx) {
+static tuple<double, double> calcLeftBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, double &sx) {
+    double centerCoeff;
+    double rightCoeff;
+
+    centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1))
+                                + 2 * alpha(rowIndex,0));
+    rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
+
+    return {centerCoeff, rightCoeff};
+}
+
+
+static tuple<double, double> calcLeftBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, double &sx) {
+    double centerCoeff;
+    double rightCoeff;
+
+    centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
+    rightCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
+
+    return {centerCoeff, rightCoeff};
+}
+
+
+static tuple<double, double> calcRightBoundaryCoeffConstant(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
+    double leftCoeff;
+    double centerCoeff;
+
+    leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
+    centerCoeff = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n))
+                                + 2 * alpha(rowIndex,n));
+
+    return {leftCoeff, centerCoeff};
+}
+
+
+static tuple<double, double> calcRightBoundaryCoeffClosed(MatrixXd &alpha, int &rowIndex, int &n, double &sx) {
+    double leftCoeff;
+    double centerCoeff;
+
+    leftCoeff = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
+    centerCoeff = 1 + sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
+
+    return {leftCoeff, centerCoeff};
+}
+
+
+static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, vector<BoundaryElement> bcLeft, vector<BoundaryElement> bcRight, int numCols, int rowIndex, double sx) {
 
     // square matrix of column^2 dimension for the coefficients
     SparseMatrix<double> cm(numCols, numCols);
     cm.reserve(VectorXi::Constant(numCols, 3));
 
     // left column
-    cm.insert(0,0) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1))
+    BC_TYPE type = bcLeft[rowIndex].getType();
-                                + 2 * alpha(rowIndex,0));
+    if (type == BC_TYPE_CONSTANT) {
-    cm.insert(0,1) = -sx * calcAlphaIntercell(alpha(rowIndex,0), alpha(rowIndex,1));
+        auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffConstant(alpha, rowIndex, sx);
+        cm.insert(0,0) = centerCoeffTop;
+        cm.insert(0,1) = rightCoeffTop;
+    } else if (type == BC_TYPE_CLOSED) { 
+        auto [centerCoeffTop, rightCoeffTop] = calcLeftBoundaryCoeffClosed(alpha, rowIndex, sx);
+        cm.insert(0,0) = centerCoeffTop;
+        cm.insert(0,1) = rightCoeffTop;
+    } else {
+        throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!");
+    }
 
     // inner columns
     int n = numCols-1;
@@ -28,22 +83,101 @@ static SparseMatrix<double> createCoeffMatrix(MatrixXd &alpha, int numCols, int
     }
 
     // right column
-    cm.insert(n,n-1) = -sx * calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n));
+    type = bcRight[rowIndex].getType();
-    cm.insert(n,n) = 1 + sx * (calcAlphaIntercell(alpha(rowIndex,n-1), alpha(rowIndex,n))
+    if (type == BC_TYPE_CONSTANT) {
-                                + 2 * alpha(rowIndex,n));
+        auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffConstant(alpha, rowIndex, n, sx);
+        cm.insert(n,n-1) = leftCoeffBottom;
+        cm.insert(n,n) = centerCoeffBottom;
+    } else if (type == BC_TYPE_CLOSED) { 
+        auto [leftCoeffBottom, centerCoeffBottom] = calcRightBoundaryCoeffClosed(alpha, rowIndex, n, sx);
+        cm.insert(n,n-1) = leftCoeffBottom;
+        cm.insert(n,n) = centerCoeffBottom;
+    } else {
+        throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!");
+    }
 
-    cm.makeCompressed();
+    cm.makeCompressed(); // important for Eigen solver
 
     return cm;
 }
 
 
+static double calcExplicitConcentrationsTopBoundaryConstant(MatrixXd &concentrations, 
+                            MatrixXd &alpha, vector<BoundaryElement> &bcTop, int &rowIndex, int &i, double &sy) {
+    double c;
+
+    c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
+                            * concentrations(rowIndex,i)
+                        + (
+                            1 - sy * (
+                                calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
+                                + 2 * alpha(rowIndex,i)
+                            )
+                        ) * concentrations(rowIndex,i)
+                        + sy * alpha(rowIndex,i) * bcTop[i].getValue();
+
+    return c;
+}
+
+
+static double calcExplicitConcentrationsTopBoundaryClosed(MatrixXd &concentrations, 
+                            MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
+    double c;
+
+    c = sy * calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
+                            * concentrations(rowIndex,i)
+                        + (
+                            1 - sy * (
+                                calcAlphaIntercell(alpha(rowIndex,i), alpha(rowIndex+1,i))
+                            )
+                        ) * concentrations(rowIndex,i);
+
+    return c;
+}
+
+
+static double calcExplicitConcentrationsBottomBoundaryConstant(MatrixXd &concentrations, 
+                            MatrixXd &alpha, vector<BoundaryElement> &bcBottom, int &rowIndex, int &i, double &sy) {
+    double c;
+
+    c = sy * alpha(rowIndex,i) * bcBottom[i].getValue()
+                        + (
+                            1 - sy * (
+                                2 * alpha(rowIndex,i)
+                                + calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
+                            )
+                        ) * concentrations(rowIndex,i)
+                        + sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) 
+                            * concentrations(rowIndex-1,i);
+
+    return c;
+}
+
+
+static double calcExplicitConcentrationsBottomBoundaryClosed(MatrixXd &concentrations, 
+                            MatrixXd &alpha, int &rowIndex, int &i, double &sy) {
+    double c;
+
+    c = (
+            1 - sy * (
+                + calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i))
+            )
+        ) * concentrations(rowIndex,i)
+        + sy * calcAlphaIntercell(alpha(rowIndex-1,i), alpha(rowIndex,i)) 
+            * concentrations(rowIndex-1,i);
+
+    return c;
+}
+
+
 static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, MatrixXd &alphaY, 
-                                        VectorXd &bcLeft, VectorXd &bcRight, VectorXd &bcTop, 
+                                        vector<BoundaryElement> &bcLeft, vector<BoundaryElement> &bcRight, 
-                                        VectorXd &bcBottom, int length, int rowIndex, double sx, double sy) {
+                                        vector<BoundaryElement> &bcTop, vector<BoundaryElement> &bcBottom, 
+                                        int length, int rowIndex, double sx, double sy) {
 
     VectorXd sv(length);
     int numRows = concentrations.rows();
+    BC_TYPE type;
 
     // inner rows
     if (rowIndex > 0 && rowIndex < numRows-1) {
@@ -65,40 +199,40 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX,
     // first row
     if (rowIndex == 0) {
         for (int i = 0; i < length; i++) {
-            sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
+            type = bcTop[i].getType();
-                            * concentrations(rowIndex,i)
+            if (type == BC_TYPE_CONSTANT) {
-                        + (
+                sv(i) = calcExplicitConcentrationsTopBoundaryConstant(concentrations, alphaY, bcTop, rowIndex, i, sy);
-                            1 - sy * (
+            } else if (type == BC_TYPE_CLOSED) {
-                                calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i))
+                sv(i) = calcExplicitConcentrationsTopBoundaryClosed(concentrations, alphaY, rowIndex, i, sy);
-                                + 2 * alphaY(rowIndex,i)
+            } else {
-                            )
+                throw_invalid_argument("Undefined Boundary Condition Type somewhere on Left or Top!");
-                        ) * concentrations(rowIndex,i)
+            }
-                        + sy * alphaY(rowIndex,i) * bcTop(i)
-                    ;
         }
     }
 
     // last row
     if (rowIndex == numRows-1) {
         for (int i = 0; i < length; i++) {
-            sv(i) = sy * alphaY(rowIndex,i) * bcBottom(i)
+            type = bcBottom[i].getType();
-                        + (
+            if (type == BC_TYPE_CONSTANT) {
-                            1 - sy * (
+                sv(i) = calcExplicitConcentrationsBottomBoundaryConstant(concentrations, alphaY, bcBottom, rowIndex, i, sy);
-                                2 * alphaY(rowIndex,i)
+            } else if (type == BC_TYPE_CLOSED) {
-                                + calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i))
+                sv(i) = calcExplicitConcentrationsBottomBoundaryClosed(concentrations, alphaY, rowIndex, i, sy);
-                            )
+            } else {
-                        ) * concentrations(rowIndex,i)
+                throw_invalid_argument("Undefined Boundary Condition Type somewhere on Right or Bottom!");
-                        + sy * calcAlphaIntercell(alphaY(rowIndex-1,i), alphaY(rowIndex,i)) 
+            }
-                            * concentrations(rowIndex-1,i)
-                    ;
         }
     }
 
-    // first column -> additional fixed concentration change from perpendicular dimension
+    // first column -> additional fixed concentration change from perpendicular dimension in constant bc case
-    sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft(rowIndex);
+    if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT) {
+        sv(0) += 2 * sx * alphaX(rowIndex,0) * bcLeft[rowIndex].getValue();
+    }
 
-    // last column -> additional fixed concentration change from perpendicular dimension
+    // last column -> additional fixed concentration change from perpendicular dimension in constant bc case
-    sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight(rowIndex);
+    if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
+        sv(length-1) += 2 * sx * alphaX(rowIndex,length-1) * bcRight[rowIndex].getValue();
+    }
 
     return sv;
 }
@@ -124,16 +258,20 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double &timestep) {
     VectorXd b(length);
 
     MatrixXd alpha = grid.getAlpha();
-    VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT);
+    vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-    VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT);
+    vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
     MatrixXd concentrations = grid.getConcentrations();
-    A = createCoeffMatrix(alpha, length, 0, sx); // this is exactly same as in 2D
+    A = createCoeffMatrix(alpha, bcLeft, bcRight, length, 0, sx); // this is exactly same as in 2D
     for (int i = 0; i < length; i++) {
-        b(i) = concentrations(0,i); // TODO check if this is really the only thing on right hand side of equation in 1D?
+        b(i) = concentrations(0,i);
+    }
+    if (bc.getBoundaryElementType(BC_SIDE_LEFT, 0) == BC_TYPE_CONSTANT) {
+        b(0) += 2 * sx * alpha(0,0) * bcLeft[0].getValue();
+    }
+    if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) {
+        b(length-1) += 2 * sx * alpha(0,length-1) * bcRight[0].getValue();
     }
-    b(0) += 2 * sx * alpha(0,0) * bcLeft(0);
-    b(length-1) += 2 * sx * alpha(0,length-1) * bcRight(0);
 
     concentrations_t1 = solve(A, b);
 
@@ -160,21 +298,21 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
 
     MatrixXd alphaX = grid.getAlphaX();
     MatrixXd alphaY = grid.getAlphaY();
-    VectorXd bcLeft = bc.getBoundarySideValues(BC_SIDE_LEFT);
+    vector<BoundaryElement> bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
-    VectorXd bcRight = bc.getBoundarySideValues(BC_SIDE_RIGHT);
+    vector<BoundaryElement> bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
-    VectorXd bcTop = bc.getBoundarySideValues(BC_SIDE_TOP);
+    vector<BoundaryElement> bcTop = bc.getBoundarySide(BC_SIDE_TOP);
-    VectorXd bcBottom = bc.getBoundarySideValues(BC_SIDE_BOTTOM);
+    vector<BoundaryElement> bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
 
     MatrixXd concentrations = grid.getConcentrations();
     for (int i = 0; i < rowMax; i++) {
       
-        A = createCoeffMatrix(alphaX, colMax, i, sx);
+        A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
         b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, 
                                     bcTop, bcBottom, colMax, i, sx, sy);
         row_t1 = solve(A, b);
         
         for (int j = 0; j < colMax; j++) {
-            concentrations_t1(i,j) = row_t1(j);
+            concentrations_t1(i,j) = row_t1(j); // TODO Eigen seq ??
         }
         
     }
@@ -184,13 +322,14 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double &timestep) {
     alphaY.transposeInPlace();
     for (int i = 0; i < colMax; i++) {
 
-        A = createCoeffMatrix(alphaY, rowMax, i, sy);
+        // swap alphas, boundary conditions and sx/sy for column-wise calculation
+        A = createCoeffMatrix(alphaY, bcLeft, bcRight, rowMax, i, sy);
         b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, 
                                     bcLeft, bcRight, rowMax, i, sy, sx);
         row_t1 = solve(A, b);
 
         for (int j = 0; j < rowMax; j++) {
-            concentrations(i,j) = row_t1(j);
+            concentrations(i,j) = row_t1(j); // TODO Eigen seq ??
         }
     }
     concentrations.transposeInPlace();