From 2e675d8139edee4e51367cc461b2178e76cae0ec Mon Sep 17 00:00:00 2001
From: Max Luebke <mluebke@gfz-potsdam.de>
Date: Thu, 15 Dec 2022 15:56:52 +0100
Subject: [PATCH] fix: setting of inner constant cells

---
 data/SimDol2D_diffu.R      | 14 +++++++++-----
 include/poet/SimParams.hpp |  2 +-
 src/DiffusionModule.cpp    | 27 +++++++++++++--------------
 src/SimParams.cpp          | 19 ++++++++++---------
 4 files changed, 33 insertions(+), 29 deletions(-)

diff --git a/data/SimDol2D_diffu.R b/data/SimDol2D_diffu.R
index b20437fd9..c495f072a 100644
--- a/data/SimDol2D_diffu.R
+++ b/data/SimDol2D_diffu.R
@@ -95,13 +95,17 @@ vecinj_diffu <- list(
   )
 )
 
-inner_index <- c(5, 15, 25)
-inner_vecinj_index <- rep(1, 3)
+#inner_index <- c(5, 15, 25)
+#inner_vecinj_index <- rep(1, 3)
+#
+#vecinj_inner <- cbind(inner_index, inner_vecinj_index)
+vecinj_inner <- list(
+  l1 = c(1,2,2)
+)
 
-vecinj_inner <- cbind(inner_index, inner_vecinj_index)
 
 boundary <- list(
-  "N" = rep(1, n),
+  "N" = rep(0, n),
   "E" = rep(0, n),
   "S" = rep(0, n),
   "W" = rep(0, n)
@@ -112,7 +116,7 @@ diffu_list <- names(alpha_diffu)
 diffusion <- list(
   init = init_diffu,
   vecinj = do.call(rbind.data.frame, vecinj_diffu),
-  # vecinj_inner = vecinj_inner,
+  vecinj_inner = vecinj_inner,
   vecinj_index = boundary,
   alpha = alpha_diffu
 )
diff --git a/include/poet/SimParams.hpp b/include/poet/SimParams.hpp
index 11a523279..3dd126a56 100644
--- a/include/poet/SimParams.hpp
+++ b/include/poet/SimParams.hpp
@@ -87,7 +87,7 @@ using DiffusionParams = struct s_DiffusionParams {
   std::vector<std::string> prop_names;
 
   Rcpp::NumericVector alpha;
-  Rcpp::NumericMatrix vecinj_inner;
+  Rcpp::List vecinj_inner;
 
   Rcpp::DataFrame vecinj;
   Rcpp::DataFrame vecinj_index;
diff --git a/src/DiffusionModule.cpp b/src/DiffusionModule.cpp
index 2c0681fec..c993e2cbf 100644
--- a/src/DiffusionModule.cpp
+++ b/src/DiffusionModule.cpp
@@ -130,29 +130,28 @@ void DiffusionModule::initialize(poet::DiffusionParams args) {
 
   // apply inner grid constant cells
   // NOTE: opening a scope here for distinguish variable names
-  if (args.vecinj_inner.rows() != 0) {
+  if (args.vecinj_inner.size() != 0) {
     // get indices of constant grid cells
-    Rcpp::NumericVector indices_const_cells = args.vecinj_inner(Rcpp::_, 0);
-    this->index_constant_cells =
-        Rcpp::as<std::vector<uint32_t>>(indices_const_cells);
+    // Rcpp::NumericVector indices_const_cells = args.vecinj_inner(Rcpp::_, 0);
+    // this->index_constant_cells =
+    //     Rcpp::as<std::vector<uint32_t>>(indices_const_cells);
 
-    // get indices to vecinj for constant cells
-    Rcpp::NumericVector vecinj_indices = args.vecinj_inner(Rcpp::_, 1);
+    // // get indices to vecinj for constant cells
+    // Rcpp::NumericVector vecinj_indices = args.vecinj_inner(Rcpp::_, 1);
 
     // apply inner constant cells for every concentration
     for (int i = 0; i < this->prop_count; i++) {
       std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
       tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
-      for (int j = 0; j < indices_const_cells.size(); j++) {
+      for (int j = 0; j < args.vecinj_inner.size(); j++) {
+        std::vector<double> inner_tuple =
+            Rcpp::as<std::vector<double>>(args.vecinj_inner[j]);
         tug::bc::boundary_condition bc = {tug::bc::BC_TYPE_CONSTANT,
-                                          bc_vec[vecinj_indices[j] - 1]};
+                                          bc_vec[inner_tuple[0] - 1]};
 
-        uint32_t x = this->index_constant_cells[j] %
-                     this->grid.getGridCellsCount(GRID_X_DIR);
-        uint32_t y = (this->dim == this->DIM_1D
-                          ? 0
-                          : this->index_constant_cells[j] /
-                                this->grid.getGridCellsCount(GRID_Y_DIR));
+        this->index_constant_cells.push_back(inner_tuple[1]);
+        uint32_t x = inner_tuple[1];
+        uint32_t y = (this->dim == this->DIM_1D ? 0 : inner_tuple[2]);
 
         curr_bc.setInnerBC(bc, x, y);
       }
diff --git a/src/SimParams.cpp b/src/SimParams.cpp
index e90ff8d41..c41649a42 100644
--- a/src/SimParams.cpp
+++ b/src/SimParams.cpp
@@ -43,9 +43,9 @@ poet::GridParams::s_GridParams(RInside &R) {
   this->props =
       Rcpp::as<std::vector<std::string>>(R.parseEval("mysetup$grid$props"));
   this->input_script =
-  Rcpp::as<std::string>(R.parseEval("mysetup$grid$input_script"));
+      Rcpp::as<std::string>(R.parseEval("mysetup$grid$input_script"));
   this->database_path =
-  Rcpp::as<std::string>(R.parseEval("mysetup$grid$database"));
+      Rcpp::as<std::string>(R.parseEval("mysetup$grid$database"));
 }
 
 poet::DiffusionParams::s_DiffusionParams(RInside &R) {
@@ -53,9 +53,10 @@ poet::DiffusionParams::s_DiffusionParams(RInside &R) {
       R.parseEval("names(mysetup$diffusion$init)"));
   this->alpha =
       Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$diffusion$alpha"));
-  if (Rcpp::as<bool>(R.parseEval("exists('mysetup$diffusion$vecinj_inner')"))) {
-    this->vecinj_inner = Rcpp::as<Rcpp::NumericMatrix>(
-        R.parseEval("mysetup$diffusion$vecinj_inner"));
+  if (Rcpp::as<bool>(
+          R.parseEval("'vecinj_inner' %in% names(mysetup$diffusion)"))) {
+    this->vecinj_inner =
+        Rcpp::as<Rcpp::List>(R.parseEval("mysetup$diffusion$vecinj_inner"));
   }
   this->vecinj =
       Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj"));
@@ -64,10 +65,10 @@ poet::DiffusionParams::s_DiffusionParams(RInside &R) {
 }
 
 poet::ChemistryParams::s_ChemistryParams(RInside &R) {
-  this->database_path = Rcpp::as<std::string>(
-      R.parseEval("mysetup$chemistry$database"));
-  this->input_script = Rcpp::as<std::string>(
-      R.parseEval("mysetup$chemistry$input_script"));
+  this->database_path =
+      Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$database"));
+  this->input_script =
+      Rcpp::as<std::string>(R.parseEval("mysetup$chemistry$input_script"));
 }
 
 SimParams::SimParams(int world_rank_, int world_size_) {