Merge branch 'boundary' into 'dev'

New use of boundary conditions

See merge request mluebke/diffusion!4

src/BTCSDiffusion.cpp

@@ -3,46 +3,101 @@
 #include <Eigen/SparseLU>
 
 #include <algorithm>
+#include <cassert>
 #include <iomanip>
 #include <iostream>
+#include <ostream>
 #include <tuple>
 #include <vector>
 
-const int BTCSDiffusion::BC_NEUMANN = 0;
+const int BTCSDiffusion::BC_CONSTANT = 0;
-const int BTCSDiffusion::BC_DIRICHLET = 1;
+const int BTCSDiffusion::BC_CLOSED = 1;
+const int BTCSDiffusion::BC_FLUX = 2;
 
-BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) {
+BTCSDiffusion::BTCSDiffusion(unsigned int dim) : grid_dim(dim) {
-  this->grid_dim = 1;
+  assert(dim <= 3);
-  this->dx = 1. / (x - 1);
 
-  // per default use Neumann condition with gradient of 0 at the end of the grid
+  grid_cells.resize(dim, 1);
-  this->bc.resize(2, std::tuple<bctype, double>(BTCSDiffusion::BC_NEUMANN, 0.));
+  domain_size.resize(dim, 1);
+  deltas.resize(dim, 1);
 }
-BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) {
 
-  // this->grid_dim = 2;
+std::vector<int> BTCSDiffusion::getNumberOfGridCells() {
+  return this->grid_cells;
+}
+std::vector<int> BTCSDiffusion::getSpatialDiscretization() {
+  return this->domain_size;
+}
+void BTCSDiffusion::setNumberOfGridCells(std::vector<int> &n_grid) {
+  grid_cells = n_grid;
+  assert(grid_cells.size() == grid_dim);
+  updateInternals();
+}
+void BTCSDiffusion::setSpatialDiscretization(std::vector<int> &s_grid) {
+  domain_size = s_grid;
+  assert(domain_size.size() == grid_dim);
+  updateInternals();
+}
+
+void BTCSDiffusion::updateInternals() {
+  for (int i = 0; i < grid_dim; i++) {
+    deltas[i] = (double)domain_size[i] / grid_cells[i];
+  }
+
+  switch (grid_dim) {
+  case 1:
+    bc.resize(2, {BTCSDiffusion::BC_CLOSED, 0});
+    break;
+  case 2:
+    bc.resize(2 * grid_cells[0] + 2 * grid_cells[1],
+              {BTCSDiffusion::BC_CLOSED, 0});
+    break;
+  case 3:
+    // TODO
+    break;
+  }
+}
+// BTCSDiffusion::BTCSDiffusion(int x) : n_x(x) {
+//   this->grid_dim = 1;
+//   this->dx = 1. / (x - 1);
 
-  // this->bc.reserve(x * 2 + y * 2);
 //   // per default use Neumann condition with gradient of 0 at the end of the
-  // grid std::fill(this->bc.begin(), this->bc.end(), -1);
+//   grid this->bc.resize(2, std::tuple<bctype,
-}
+//   double>(BTCSDiffusion::BC_CONSTANT, 0.));
-BTCSDiffusion::BTCSDiffusion(int x, int y, int z) : n_x(x), n_y(y), n_z(z) {
+// }
+// BTCSDiffusion::BTCSDiffusion(int x, int y) : n_x(x), n_y(y) {
 
-  // this->grid_dim = 3;
+//   // this->grid_dim = 2;
-  // TODO: reserve memory for boundary conditions
-}
 
-void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
+//   // this->bc.reserve(x * 2 + y * 2);
-                               double bc_right,
+//   // // per default use Neumann condition with gradient of 0 at the end of
+//   the
+//   // grid std::fill(this->bc.begin(), this->bc.end(), -1);
+// }
+// BTCSDiffusion::BTCSDiffusion(int x, int y, int z) : n_x(x), n_y(y), n_z(z) {
+
+//   // this->grid_dim = 3;
+//   // TODO: reserve memory for boundary conditions
+// }
+
+void BTCSDiffusion::simulate1D(std::vector<double> &c, boundary_condition left,
+                               boundary_condition right,
                                const std::vector<double> &alpha, double dx,
                                int size) {
 
+  bool left_is_constant = (left.type == BTCSDiffusion::BC_CONSTANT);
+  bool right_is_constant = (right.type == BTCSDiffusion::BC_CONSTANT);
+  int loop_end = size + !right_is_constant;
+
   // we need 2 more grid cells for ghost cells
-  size = size + 2;
+  // size = size + 2;
+
+  int bc_offset = !left_is_constant + !right_is_constant;
+  ;
 
   // set sizes of private and yet allocated vectors
-  b_vector.resize(size);
+  b_vector.resize(size + bc_offset);
-  x_vector.resize(size);
+  x_vector.resize(size + bc_offset);
 
   /*
    * Begin to solve the equation system using LU solver of Eigen.
@@ -53,25 +108,33 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
    * TODO: remove output
    */
 
-  A_matrix.resize(size, size);
+  A_matrix.resize(size + bc_offset, size + bc_offset);
-  A_matrix.reserve(Eigen::VectorXi::Constant(size, 3));
+  A_matrix.reserve(Eigen::VectorXi::Constant(size + bc_offset, 3));
 
   A_matrix.insert(0, 0) = 1;
-  A_matrix.insert(size - 1, size - 1) = 1;
+  b_vector[0] =
+      (left_is_constant ? left.value : getBCFromFlux(left, c[0], alpha[0]));
 
-  b_vector[0] = bc_left;
+  A_matrix.insert((size + bc_offset) - 1, (size + bc_offset) - 1) = 1;
-  b_vector[size - 1] = bc_right;
+  b_vector[size + bc_offset - 1] =
+      (right_is_constant ? right.value
+                         : getBCFromFlux(right, c[size - 1], alpha[size - 1]));
 
-  for (int i = 1; i < this->n_x + 1; i++) {
+  // A_matrix.insert(0, 0) = 1;
-    double sx = (alpha[i - 1] * time_step) / (dx * dx);
+  // A_matrix.insert(size + 1, size + 1) = 1;
+
+  for (int i = 1; i < size - right_is_constant; i++) {
+    double sx = (alpha[i + !(left_is_constant)] * 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[i - 1];
+    b_vector[i] = -c[i + !(left_is_constant)];
   }
 
+  std::cout << b_vector << "\n" << A_matrix << std::endl;
+
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       solver;
   solver.analyzePattern(A_matrix);
@@ -85,7 +148,7 @@ void BTCSDiffusion::simulate1D(std::vector<double> &c, double bc_left,
   std::cout << std::setprecision(10) << x_vector << std::endl << std::endl;
 
   for (int i = 0; i < c.size(); i++) {
-    c[i] = x_vector[i + 1];
+    c[i] = x_vector[i + !left_is_constant];
   }
 }
 
@@ -96,24 +159,25 @@ void BTCSDiffusion::setTimestep(double time_step) {
 void BTCSDiffusion::simulate(std::vector<double> &c,
                              const std::vector<double> &alpha) {
   if (this->grid_dim == 1) {
-    double bc_left = getBCFromTuple(0, c[0], alpha[0]);
+    // double bc_left = getBCFromTuple(0, c[0], alpha[0]);
-    double bc_right =
+    // double bc_right =
-        getBCFromTuple(1, c[c.size() - 1], alpha[alpha.size() - 1]);
+    //     getBCFromTuple(1, c[c.size() - 1], alpha[alpha.size() - 1]);
 
-    simulate1D(c, bc_left, bc_right, alpha, this->dx, this->n_x);
+    simulate1D(c, bc[0], bc[1], alpha, this->deltas[0], this->grid_cells[0]);
   }
 }
 
-double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
+inline double BTCSDiffusion::getBCFromFlux(boundary_condition bc,
+                                           double neighbor_c,
                                            double neighbor_alpha) {
-  double val = -1;
-  int type = std::get<0>(bc[index]);
 
-  if (type == BTCSDiffusion::BC_NEUMANN) {
+  double val;
-    val = neighbor_c + (this->time_step / (dx * dx)) * neighbor_alpha *
+
-                           std::get<1>(bc[index]);
+  if (bc.type == BTCSDiffusion::BC_CLOSED) {
-  } else if (type == BTCSDiffusion::BC_DIRICHLET) {
+    val = neighbor_c;
-    val = std::get<1>(bc[index]);
+  } else if (bc.type == BTCSDiffusion::BC_FLUX) {
+    // TODO
+    //  val = bc[index].value;
   } else {
     // TODO: implement error handling here. Type was set to wrong value.
   }
@@ -122,6 +186,10 @@ double BTCSDiffusion::getBCFromTuple(int index, double neighbor_c,
 }
 
 void BTCSDiffusion::setBoundaryCondition(int index, double val, bctype type) {
-  std::get<0>(bc[index]) = type;
+
-  std::get<1>(bc[index]) = val;
+  bc[index].type = type;
+  bc[index].value = val;
+
+  // std::get<0>(bc[index]) = type;
+  // std::get<1>(bc[index]) = val;
 }

src/BTCSDiffusion.hpp

@@ -5,17 +5,33 @@
 #include <tuple>
 #include <vector>
 
-/*!
- * Datatype to fill the sparse matrix which is used to solve the equation
- * system.
- */
-typedef Eigen::Triplet<double> T;
-
 /*!
  * Defines both types of boundary condition as a datatype.
  */
 typedef int bctype;
 
+/*!
+ * Class implementing a solution for a 1/2/3D diffusion equation using backward
+ * euler.
+ */
+class BTCSDiffusion {
+
+public:
+  /*!
+   * Defines a constant/Dirichlet boundary condition.
+   */
+  static const int BC_CONSTANT;
+
+  /*!
+   * Defines a closed/Neumann boundary condition.
+   */
+  static const int BC_CLOSED;
+
+  /*!
+   * Defines a flux/Cauchy boundary condition.
+   */
+  static const int BC_FLUX;
+
   /*!
    * A boundary condition consists of two features. A type and the according
    * value. Here we can differentiate between:
@@ -25,47 +41,56 @@ typedef int bctype;
    * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
    * the boundary condition
    */
-typedef std::vector<std::tuple<bctype, double>> boundary_condition;
+  typedef struct boundary_condition {
+    bctype type;
+    double value;
+  } boundary_condition;
 
   /*!
- * Class implementing a solution for a 1/2/3D diffusion equation using backward
+   * A boundary condition consists of two features. A type and the according
- * euler.
+   * value. Here we can differentiate between:
+   *
+   * - Neumann boundary conditon: type BC_NEUMANN with the value defining the
+   * gradient
+   * - Dirichlet boundary condition: type BC_DIRICHLET with the actual value of
+   * the boundary condition
    */
-class BTCSDiffusion {
+  // typedef std::vector<std::tuple<bctype, double>> boundary_condition;
 
-public:
   /*!
-   * Defines a Neumann boundary condition.
+   * Datatype to fill the sparse matrix which is used to solve the equation
+   * system.
    */
-  static const int BC_NEUMANN;
+  typedef Eigen::Triplet<double> T;
-  /*!
-   * Defines a Dirichlet boundary condition.
-   */
-  static const int BC_DIRICHLET;
 
   /*!
    * Create 1D-diffusion module.
    *
    * @param x Count of cells in x direction.
    */
-  explicit BTCSDiffusion(int x);
+  BTCSDiffusion(unsigned int dim);
 
-  /*!
+  std::vector<int> getNumberOfGridCells();
-   * Currently not implemented: Create 2D-diffusion module.
+  std::vector<int> getSpatialDiscretization();
-   *
+  void setNumberOfGridCells(std::vector<int> &n_grid);
-   * @param x Count of cells in x direction.
+  void setSpatialDiscretization(std::vector<int> &s_grid);
-   * @param y Count of cells in y direction.
-   */
-  explicit BTCSDiffusion(int x, int y);
 
-  /*!
+  // /*!
-   * Currently not implemented: Create 3D-diffusion module.
+  //  * Currently not implemented: Create 2D-diffusion module.
-   *
+  //  *
-   * @param x Count of cells in x direction.
+  //  * @param x Count of cells in x direction.
-   * @param y Count of cells in y direction.
+  //  * @param y Count of cells in y direction.
-   * @param z Count of cells in z direction.
+  //  */
-   */
+  // explicit BTCSDiffusion(int x, int y);
-  explicit BTCSDiffusion(int x, int y, int z);
+
+  // /*!
+  //  * Currently not implemented: Create 3D-diffusion module.
+  //  *
+  //  * @param x Count of cells in x direction.
+  //  * @param y Count of cells in y direction.
+  //  * @param z Count of cells in z direction.
+  //  */
+  // explicit BTCSDiffusion(int x, int y, int z);
 
   /*!
    * With given ghost zones simulate diffusion. Only 1D allowed at this moment.
@@ -96,14 +121,17 @@ public:
   void setBoundaryCondition(int index, double val, bctype type);
 
 private:
-  void simulate1D(std::vector<double> &c, double bc_left, double bc_right,
+  void simulate1D(std::vector<double> &c, boundary_condition left,
-                  const std::vector<double> &alpha, double dx, int size);
+                  boundary_condition right, const std::vector<double> &alpha,
+                  double dx, int size);
   void simulate2D(std::vector<double> &c);
   void simulate3D(std::vector<double> &c);
 
-  double getBCFromTuple(int index, double nearest_value, double neighbor_alpha);
+  inline double getBCFromFlux(boundary_condition bc, double nearest_value, double neighbor_alpha);
 
-  boundary_condition bc;
+  void updateInternals();
+
+  std::vector<boundary_condition> bc;
 
   Eigen::SparseMatrix<double> A_matrix;
   Eigen::VectorXd b_vector;
@@ -112,12 +140,9 @@ private:
   double time_step;
 
   int grid_dim;
-  int n_x;
+  std::vector<int> grid_cells;
-  double dx;
+  std::vector<int> domain_size;
-  int n_y;
+  std::vector<double> deltas;
-  double dy;
-  int n_z;
-  double dz;
 };
 
 #endif // BTCSDIFFUSION_H_

src/main.cpp

@@ -6,19 +6,27 @@ using namespace std;
 
 int main(int argc, char *argv[]) {
 
-  // count of grid cells
+  // dimension of grid
-  int x = 20;
+  int dim = 1;
+
+  int n = 20;
 
   // create input + diffusion coefficients for each grid cell
-  std::vector<double> alpha(x, 1 * pow(10, -1));
+  std::vector<double> alpha(n, 1 * pow(10, -1));
-  std::vector<double> input(x, 1 * std::pow(10, -6));
+  std::vector<double> field(n, 1 * std::pow(10, -6));
 
   // create instance of diffusion module
-  BTCSDiffusion diffu(x);
+  BTCSDiffusion diffu(dim);
+
+  std::vector<int> vec_n = diffu.getNumberOfGridCells();
+
+  vec_n[0] = n;
+
+  diffu.setNumberOfGridCells(vec_n);
 
   // set the boundary condition for the left ghost cell to dirichlet
   diffu.setBoundaryCondition(0, 5. * std::pow(10, -6),
-                             BTCSDiffusion::BC_DIRICHLET);
+                             BTCSDiffusion::BC_CONSTANT);
 
   // set timestep for simulation to 1 second
   diffu.setTimestep(1.);
@@ -26,7 +34,7 @@ int main(int argc, char *argv[]) {
   // loop 100 times
   // output is currently generated by the method itself
   for (int i = 0; i < 100; i++) {
-    diffu.simulate(input, alpha);
+    diffu.simulate(field, alpha);
   }
 
   return 0;