Use simulate_base for actual solving of les.

src/BTCSDiffusion.cpp

@@ -76,74 +76,22 @@ void Diffusion::BTCSDiffusion::updateInternals() {
     deltas[i] = (double)domain_size[i] / grid_cells[i];
   }
 }
+void Diffusion::BTCSDiffusion::simulate_base(DVectorRowMajor &c,
+                                             const BCVectorRowMajor &bc,
+                                             const DVectorRowMajor &alpha,
+                                             double dx, double time_step,
+                                             int size,
+                                             const DVectorRowMajor &t0_c) {
 
-void Diffusion::BTCSDiffusion::simulate_base(
+  reserveMemory(size, BTCS_MAX_DEP_PER_CELL);
-    DVectorRowMajor &c, Eigen::Map<const BCVectorRowMajor> &bc,
-    Eigen::Map<const DVectorRowMajor> &alpha, double dx, double time_step,
-    int size, DVectorRowMajor &t0_c) {
 
-  // The sizes for matrix and vectors of the equation system is defined by the
+  fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step);
-  // actual size of the input vector and if the system is (partially) closed.
+  fillVectorFromRowADI(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size,
-  // Then we will need ghost nodes. So this variable will give the count of
+                       dx, time_step);
-  // ghost nodes.
-  // int bc_offset = !left_is_constant + !right_is_constant;
-  // ;
-
-  // set sizes of private and yet allocated vectors
-  // b_vector.resize(size + bc_offset);
-  // x_vector.resize(size + bc_offset);
-
-  // /*
-  //  * Begin to solve the equation system using LU solver of Eigen.
-  //  *
-  //  * But first fill the A matrix and b vector.
-  //  */
-
-  // // Set boundary condition for ghost nodes (for closed or flux system) or
-  // outer
-  // // inlet nodes (constant boundary condition)
-  // A_matrix.resize(size + bc_offset, size + bc_offset);
-  // A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3));
-
-  // A_matrix.insert(0, 0) = 1;
-  // b_vector[0] =
-  //     (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
-
-  // A_matrix.insert(size + 1, size + 1) = 1;
-  // b_vector[size + 1] =
-  //     (right_is_constant ? right.value
-  //                        : getBCFromFlux(right, c[size - 1], alpha[size -
-  //                        1]));
-
-  // Start filling the A matrix
-  // =i= is used for equation system matrix and vector indexing
-  // and =j= for indexing of c,alpha and bc
-  // for (int i = 1, j = i + !(left_is_constant); i < size - right_is_constant;
-  //      i++, j++) {
-
-  //   // if current grid cell is considered as constant boundary conditon
-  //   if (bc[j].type == Diffusion::BC_CONSTANT) {
-  //     A_matrix.insert(i, i) = 1;
-  //     b_vector[i] = bc[j].value;
-  //     continue;
-  //   }
-
-  //   double sx = (alpha[j] * time_step) / (dx * dx);
-
-  //   A_matrix.insert(i, i) = -1. - 2. * sx;
-  //   A_matrix.insert(i, i - 1) = sx;
-  //   A_matrix.insert(i, i + 1) = sx;
-
-  //   b_vector[i] = -c[j];
-  // }
-
-  fillMatrixFromRow(alpha, bc, size, dx, time_step);
-  fillVectorFromRowADI(c, alpha, bc, t0_c, size, dx, time_step);
 
   solveLES();
 
-  // write back result to input/output vector
+  c = x_vector.segment(1, size);
-  // c = x_vector.segment(!left_is_constant, c.size());
 }
 
 inline void Diffusion::BTCSDiffusion::reserveMemory(int size,
@@ -165,15 +113,12 @@ void Diffusion::BTCSDiffusion::simulate1D(
   double dx = this->deltas[0];
   double time_step = this->time_step;
 
-  reserveMemory(size, BTCS_MAX_DEP_PER_CELL);
+  DVectorRowMajor input_field = c.row(0);
 
-  fillMatrixFromRow(alpha.row(0), bc.row(0), size, dx, time_step);
+  simulate_base(input_field, bc, alpha, dx, time_step, size,
-  fillVectorFromRowADI(c, alpha, bc, Eigen::VectorXd::Constant(size, 0), size,
+                Eigen::VectorXd::Constant(size, 0));
-                       dx, time_step);
 
-  solveLES();
+  c.row(0) << input_field;
-
-  c = x_vector.segment(1, size);
 }
 
 void Diffusion::BTCSDiffusion::simulate2D(

src/BTCSDiffusion.hpp

@@ -115,9 +115,9 @@ private:
                         Eigen::RowMajor>
       BCVectorRowMajor;
 
-  void simulate_base(DVectorRowMajor &c, Eigen::Map<const BCVectorRowMajor> &bc,
+  void simulate_base(DVectorRowMajor &c, const BCVectorRowMajor &bc,
-                     Eigen::Map<const DVectorRowMajor> &alpha, double dx,
+                     const DVectorRowMajor &alpha, double dx,
-                     double time_step, int size, DVectorRowMajor &t0_c);
+                     double time_step, int size, const DVectorRowMajor &t0_c);
 
   void simulate1D(Eigen::Map<DVectorRowMajor> &c,
                   Eigen::Map<const DVectorRowMajor> &alpha,