tug: cpp: Remove inner boundary

src/benchmarks/tug/cpp/include/tug/Boundary.hpp

@@ -376,174 +376,6 @@ public:
     }
   }
 
-  /**
-   *
-   * @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) {
-    if (this->dim != 1) {
-      throw std::invalid_argument(
-          "This function is only available for 1D grids.");
-    }
-    if (index >= this->cols) {
-      throw std::invalid_argument("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) {
-    if (this->dim != 2) {
-      throw std::invalid_argument(
-          "This function is only available for 2D grids.");
-    }
-    if (row >= this->rows || col >= this->cols) {
-      throw std::invalid_argument("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 {
-    if (this->dim != 1) {
-      throw std::invalid_argument(
-          "This function is only available for 1D grids.");
-    }
-    if (index >= this->cols) {
-      throw std::invalid_argument("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 {
-    if (this->dim != 2) {
-      throw std::invalid_argument(
-          "This function is only available for 2D grids.");
-    }
-    if (row >= this->rows || col >= this->cols) {
-      throw std::invalid_argument("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 {
-    if (row >= this->rows) {
-      throw std::invalid_argument("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 {
-    if (this->dim != 2) {
-      throw std::invalid_argument(
-          "This function is only available for 2D grids.");
-    }
-
-    if (col >= this->cols) {
-      throw std::invalid_argument("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;
@@ -551,9 +383,6 @@ 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_

src/benchmarks/tug/cpp/include/tug/Core/BTCS.hpp

@@ -55,7 +55,7 @@ static Eigen::SparseMatrix<T>
 createCoeffMatrix(const RowMajMat<T> &alpha,
                   const std::vector<BoundaryElement<T>> &bcLeft,
                   const std::vector<BoundaryElement<T>> &bcRight,
-                  const std::vector<std::pair<bool, T>> &inner_bc, int numCols,
+                  int numCols,
                   int rowIndex, T sx) {
 
   // square matrix of column^2 dimension for the coefficients
@@ -63,9 +63,7 @@ createCoeffMatrix(const RowMajMat<T> &alpha,
   cm.reserve(Eigen::VectorXi::Constant(numCols, 3));
 
   // left column
-  if (inner_bc[0].first) {
-    cm.insert(0, 0) = 1;
-  } else {
+  {
     switch (bcLeft[rowIndex].getType()) {
     case BC_TYPE_CONSTANT: {
       auto [centerCoeffTop, rightCoeffTop] =
@@ -91,10 +89,6 @@ createCoeffMatrix(const RowMajMat<T> &alpha,
   // inner columns
   int n = numCols - 1;
   for (int i = 1; i < n; i++) {
-    if (inner_bc[i].first) {
-      cm.insert(i, i) = 1;
-      continue;
-    }
     cm.insert(i, i - 1) =
         -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i));
     cm.insert(i, i) =
@@ -106,9 +100,7 @@ createCoeffMatrix(const RowMajMat<T> &alpha,
   }
 
   // right column
-  if (inner_bc[n].first) {
-    cm.insert(n, n) = 1;
-  } else {
+  {
     switch (bcRight[rowIndex].getType()) {
     case BC_TYPE_CONSTANT: {
       auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant(
@@ -167,7 +159,6 @@ createSolutionVector(const RowMajMat<T> &concentrations,
                      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);
@@ -176,10 +167,6 @@ createSolutionVector(const RowMajMat<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)) *
@@ -198,10 +185,6 @@ createSolutionVector(const RowMajMat<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(
@@ -225,10 +208,6 @@ createSolutionVector(const RowMajMat<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(
@@ -251,14 +230,13 @@ createSolutionVector(const RowMajMat<T> &concentrations,
 
   // first column -> additional fixed concentration change from perpendicular
   // dimension in constant bc case
-  if (bcLeft[rowIndex].getType() == BC_TYPE_CONSTANT && !inner_bc[0].first) {
+  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 in constant bc case
-  if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT &&
-      !inner_bc[length - 1].first) {
+  if (bcRight[rowIndex].getType() == BC_TYPE_CONSTANT) {
     sv(length - 1) +=
         2 * sx * alphaX(rowIndex, length - 1) * bcRight[rowIndex].getValue();
   }
@@ -367,21 +345,17 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
   const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
   const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
 
-  const auto inner_bc = bc.getInnerBoundaryRow(0);
-
   RowMajMat<T> &concentrations = grid.getConcentrations();
   int rowIndex = 0;
-  A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex,
+  A = createCoeffMatrix(alpha, bcLeft, bcRight, 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 &&
-      !inner_bc[0].first) {
+  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 &&
-      !inner_bc[length - 1].first) {
+  if (bc.getBoundaryElementType(BC_SIDE_RIGHT, 0) == BC_TYPE_CONSTANT) {
     b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue();
   }
 
@@ -421,11 +395,9 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
 
 #pragma omp parallel for private(A, b)
   for (int i = 0; i < rowMax; i++) {
-    auto inner_bc = bc.getInnerBoundaryRow(i);
-
-    A = createCoeffMatrix(alphaX, bcLeft, bcRight, inner_bc, colMax, i, sx);
+    A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx);
     b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
-                             bcTop, bcBottom, inner_bc, colMax, i, sx, sy);
+                             bcTop, bcBottom, colMax, i, sx, sy);
 
     concentrations_t1.row(i) = solverFunc(A, b);
   }
@@ -436,11 +408,10 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
 
 #pragma omp parallel for private(A, b)
   for (int i = 0; i < colMax; i++) {
-    auto inner_bc = bc.getInnerBoundaryCol(i);
     // swap alphas, boundary conditions and sx/sy for column-wise calculation
-    A = createCoeffMatrix(alphaY, bcTop, bcBottom, inner_bc, rowMax, i, sy);
+    A = createCoeffMatrix(alphaY, bcTop, bcBottom, rowMax, i, sy);
     b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
-                             bcLeft, bcRight, inner_bc, rowMax, i, sy, sx);
+                             bcLeft, bcRight, rowMax, i, sy, sx);
 
     concentrations.col(i) = solverFunc(A, b);
   }