From 031c1b2eef27668b130b48d1ad6e24498ff4e62d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Max=20L=C3=BCbke?= <mluebke@uni-potsdam.de>
Date: Wed, 11 Dec 2024 13:22:23 +0100
Subject: [PATCH] feat: Implement inner boundaries for FTCS

---
 include/tug/Core/Numeric/FTCS.hpp | 24 ++++++++++++++++++++++++
 test/testDiffusion.cpp            | 30 ++++++++++++++++++++++++++++++
 2 files changed, 54 insertions(+)

diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp
index ce1d5e8..3ab5cc2 100644
--- a/include/tug/Core/Numeric/FTCS.hpp
+++ b/include/tug/Core/Numeric/FTCS.hpp
@@ -11,6 +11,7 @@
 #include "tug/Core/TugUtils.hpp"
 #include <cstddef>
 #include <cstring>
+#include <iterator>
 #include <tug/Boundary.hpp>
 #include <tug/Core/Matrix.hpp>
 #include <tug/Core/Numeric/SimulationInput.hpp>
@@ -55,6 +56,27 @@ constexpr T calcChangeBoundary(T conc_c, T conc_neighbor, T alpha_center,
   tug_assert(false, "Undefined Boundary Condition Type!");
 }
 
+template <typename T>
+static inline void checkAndSetConstantInnerCells(const Boundary<T> &bc,
+                                                 Grid<T> &grid) {
+  const auto &inner_boundaries = bc.getInnerBoundaries();
+  if (inner_boundaries.empty()) {
+    return;
+  }
+
+  auto &concentrations = grid.getConcentrations();
+  const auto rows = grid.getRow();
+  const auto cols = grid.getCol();
+
+  for (const auto &[rowcol, value] : inner_boundaries) {
+    const auto &row = rowcol.first;
+    const auto &col = rowcol.second;
+    concentrations(row, col) = value;
+  }
+
+  return;
+}
+
 // FTCS solution for 1D grid
 template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   const std::size_t &colMax = input.colMax;
@@ -71,6 +93,8 @@ template <class T> static void FTCS_1D(SimulationInput<T> &input) {
   // only one row in 1D case -> row constant at index 0
   int row = 0;
 
+  const auto inner_bc = bc.getInnerBoundaryRow(0);
+
   // inner cells
   // independent of boundary condition type
   for (int col = 1; col < colMax - 1; col++) {
diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp
index c49a166..9a74670 100644
--- a/test/testDiffusion.cpp
+++ b/test/testDiffusion.cpp
@@ -215,3 +215,33 @@ DIFFUSION_TEST(ConstantInnerCell) {
 
   EXPECT_FALSE((concentrations_result.array() < 0.0).any());
 }
+
+DIFFUSION_TEST(ConstantInnerCellFTCS) {
+  constexpr std::uint32_t nrows = 5;
+  constexpr std::uint32_t ncols = 5;
+
+  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);
+
+  tug::Grid64 grid(concentrations);
+  grid.setAlpha(alphax, alphay);
+
+  tug::Boundary bc(grid);
+  // inner
+  bc.setInnerBoundary(2, 2, 0);
+
+  tug::Diffusion<double, tug::FTCS_APPROACH> sim(grid, bc);
+  sim.setTimestep(1);
+  sim.setIterations(1);
+
+  MatrixXd input_values(concentrations);
+  sim.run();
+
+  EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 2), 0);
+  EXPECT_LT(grid.getConcentrations().sum(), input_values.sum());
+
+  EXPECT_FALSE((grid.getConcentrations().array() > 1.0).any());
+
+  EXPECT_FALSE((grid.getConcentrations().array() < 0.0).any());
+}
\ No newline at end of file