feat: allow undefined boundary conditions

It is now possible to not define any boundary condition object.

In this case a grid with closed boundary conditions is assumed during
diffusion simulation.

refactor: various adjustments according to clang-tidy

include/BoundaryCondition.hpp

@@ -166,7 +166,7 @@ private:
 
   uint8_t dim;
 
-  uint32_t sizes[2];
+  std::array<uint32_t, 2> sizes;
   uint32_t maxsize;
   uint32_t maxindex;
 
@@ -260,6 +260,6 @@ public:
   }
 };
 
-} // namespace boundary_condition
+} // namespace bc
 } // namespace tug
 #endif // BOUNDARYCONDITION_H_

include/Diffusion.hpp

@@ -12,14 +12,17 @@
 namespace tug {
 namespace diffusion {
 
+constexpr uint8_t MAX_ARR_SIZE = 3;
+
 /**
  * Defines grid dimensions and boundary conditions.
  */
-typedef struct {
+typedef struct tug_grid_s {
-  uint32_t
+  std::array<uint32_t, MAX_ARR_SIZE>
-      grid_cells[3];     /**< Count of grid cells in each of the 3 directions.*/
+      grid_cells; /**< Count of grid cells in each of the 3 directions.*/
-  double domain_size[3]; /**< Domain sizes in each of the 3 directions.*/
+  std::array<double, MAX_ARR_SIZE>
-  bc::BoundaryCondition *bc; /**< Boundary conditions for the grid.*/
+      domain_size; /**< Domain sizes in each of the 3 directions.*/
+  bc::BoundaryCondition *bc = NULL; /**< Boundary conditions for the grid.*/
 } TugGrid;
 
 /**
@@ -28,7 +31,8 @@ typedef struct {
  */
 typedef struct tug_input_s {
   double time_step; /**< Time step which should be simulated by diffusion.*/
-  Eigen::VectorXd (*solver)(Eigen::SparseMatrix<double>, Eigen::VectorXd) =
+  Eigen::VectorXd (*solver)(const Eigen::SparseMatrix<double> &,
+                            const Eigen::VectorXd &) =
       tug::solver::ThomasAlgorithm; /**< Solver function to use.*/
   TugGrid grid;                     /**< Grid specification.*/
 
@@ -86,8 +90,9 @@ typedef struct tug_input_s {
    * \param f_in Pointer to function which takes a sparse matrix and a vector as
    * input and returns another vector.
    */
-  void setSolverFunction(Eigen::VectorXd (*f_in)(Eigen::SparseMatrix<double>,
+  void
-                                                 Eigen::VectorXd)) {
+  setSolverFunction(Eigen::VectorXd (*f_in)(const Eigen::SparseMatrix<double> &,
+                                            const Eigen::VectorXd &)) {
     solver = f_in;
   }
 } TugInput;

include/Solver.hpp

@@ -18,8 +18,8 @@ namespace solver {
  *
  * \return Solution represented as vector.
  */
-auto EigenLU(const Eigen::SparseMatrix<double> A_matrix,
+auto EigenLU(const Eigen::SparseMatrix<double> &A_matrix,
-             const Eigen::VectorXd b_vector) -> Eigen::VectorXd;
+             const Eigen::VectorXd &b_vector) -> Eigen::VectorXd;
 
 /**
  * Solving linear equation system with brutal implementation of the Thomas
@@ -32,8 +32,8 @@ auto EigenLU(const Eigen::SparseMatrix<double> A_matrix,
  *
  * \return Solution represented as vector.
  */
-auto ThomasAlgorithm(const Eigen::SparseMatrix<double> A_matrix,
+auto ThomasAlgorithm(const Eigen::SparseMatrix<double> &A_matrix,
-                     const Eigen::VectorXd b_vector) -> Eigen::VectorXd;
+                     const Eigen::VectorXd &b_vector) -> Eigen::VectorXd;
 
 } // namespace solver
 } // namespace tug

src/BTCS.cpp

@@ -1,5 +1,6 @@
 #include <Diffusion.hpp>
 #include <Solver.hpp>
+#include <array>
 #include <iostream>
 
 #include <Eigen/Dense>
@@ -7,6 +8,7 @@
 #include <Eigen/src/Core/Matrix.h>
 #include <bits/stdint-uintn.h>
 #include <chrono>
+#include <vector>
 
 #include "BoundaryCondition.hpp"
 #include "TugUtils.hpp"
@@ -17,12 +19,16 @@
 #define omp_get_thread_num() 0
 #endif
 
-#define init_delta(in, out, dim)                                               \
+inline auto
-  ({                                                                           \
+init_delta(const std::array<double, tug::diffusion::MAX_ARR_SIZE> &domain_size,
-    for (uint8_t i = 0; i < dim; i++) {                                        \
+           const std::array<uint32_t, tug::diffusion::MAX_ARR_SIZE> &grid_cells,
-      out[i] = (double)in.domain_size[i] / in.grid_cells[i];                   \
+           const uint8_t dim) -> std::vector<double> {
-    }                                                                          \
+  std::vector<double> out(dim);
-  })
+  for (uint8_t i = 0; i < dim; i++) {
+    out[i] = (double)(domain_size.at(i) / grid_cells.at(i));
+  }
+  return out;
+}
 
 namespace {
 enum { GRID_1D = 1, GRID_2D, GRID_3D };
@@ -30,10 +36,10 @@ enum { GRID_1D = 1, GRID_2D, GRID_3D };
 constexpr int BTCS_MAX_DEP_PER_CELL = 3;
 constexpr int BTCS_2D_DT_SIZE = 2;
 
-typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
+using DMatrixRowMajor =
-    DMatrixRowMajor;
+    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-typedef Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>
+using DVectorRowMajor =
-    DVectorRowMajor;
+    Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>;
 
 inline auto getBCFromFlux(tug::bc::boundary_condition bc, double neighbor_c,
                           double neighbor_alpha) -> double {
@@ -192,9 +198,10 @@ auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha,
 
   for (int j = 0; j < size; j++) {
     if (bc_inner[j].type != tug::bc::BC_UNSET) {
-      if (bc_inner[j].type != tug::bc::BC_TYPE_CONSTANT)
+      if (bc_inner[j].type != tug::bc::BC_TYPE_CONSTANT) {
         throw_invalid_argument("Inner boundary conditions with other type than "
                                "BC_TYPE_CONSTANT are currently not supported.");
+      }
       b_vector[j + 1] = bc_inner[j].value;
       continue;
     }
@@ -215,10 +222,10 @@ auto fillVectorFromRow(const DVectorRowMajor &c, const DVectorRowMajor &alpha,
 
 auto setupBTCSAndSolve(
     DVectorRowMajor &c, const tug::bc::bc_tuple bc_ghost,
-    const tug::bc::bc_vec bc_inner, const DVectorRowMajor &alpha, double dx,
+    const tug::bc::bc_vec &bc_inner, const DVectorRowMajor &alpha, double dx,
     double time_step, int size, const DVectorRowMajor &d_ortho,
-    Eigen::VectorXd (*solver)(Eigen::SparseMatrix<double>, Eigen::VectorXd))
+    Eigen::VectorXd (*solver)(const Eigen::SparseMatrix<double> &,
-    -> DVectorRowMajor {
+                              const Eigen::VectorXd &)) -> DVectorRowMajor {
 
   const Eigen::SparseMatrix<double> A_matrix =
       fillMatrixFromRow(alpha, bc_inner, size, dx, time_step);
@@ -241,15 +248,18 @@ auto tug::diffusion::BTCS_1D(const tug::diffusion::TugInput &input_param,
 
   auto start = time_marker();
 
-  const tug::bc::BoundaryCondition bc = *input_param.grid.bc;
-  double deltas[GRID_1D];
-
-  init_delta(input_param.grid, deltas, GRID_1D);
-
   uint32_t size = input_param.grid.grid_cells[0];
+
+  auto deltas = init_delta(input_param.grid.domain_size,
+                           input_param.grid.grid_cells, GRID_1D);
   double dx = deltas[0];
+
   double time_step = input_param.time_step;
 
+  const tug::bc::BoundaryCondition bc =
+      (input_param.grid.bc != nullptr ? *input_param.grid.bc
+                                      : tug::bc::BoundaryCondition(size));
+
   Eigen::Map<DVectorRowMajor> c_in(field, size);
   Eigen::Map<const DVectorRowMajor> alpha_in(alpha, size);
 
@@ -271,17 +281,21 @@ auto tug::diffusion::ADI_2D(const tug::diffusion::TugInput &input_param,
 
   auto start = time_marker();
 
-  tug::bc::BoundaryCondition bc = *input_param.grid.bc;
-  double deltas[GRID_2D];
-
-  init_delta(input_param.grid, deltas, GRID_2D);
-
   uint32_t n_cols = input_param.grid.grid_cells[0];
   uint32_t n_rows = input_param.grid.grid_cells[1];
+
+  auto deltas = init_delta(input_param.grid.domain_size,
+                           input_param.grid.grid_cells, GRID_2D);
   double dx = deltas[0];
   double dy = deltas[1];
+
   double local_dt = input_param.time_step / BTCS_2D_DT_SIZE;
 
+  tug::bc::BoundaryCondition bc =
+      (input_param.grid.bc != nullptr
+           ? *input_param.grid.bc
+           : tug::bc::BoundaryCondition(n_cols, n_rows));
+
   Eigen::Map<DMatrixRowMajor> c_in(field, n_rows, n_cols);
   Eigen::Map<const DMatrixRowMajor> alpha_in(alpha, n_rows, n_cols);
 

src/Solver.cpp

@@ -4,8 +4,8 @@
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
-auto tug::solver::EigenLU(const Eigen::SparseMatrix<double> A_matrix,
+auto tug::solver::EigenLU(const Eigen::SparseMatrix<double> &A_matrix,
-                          const Eigen::VectorXd b_vector) -> Eigen::VectorXd {
+                          const Eigen::VectorXd &b_vector) -> Eigen::VectorXd {
 
   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>
       solver;
@@ -18,8 +18,8 @@ auto tug::solver::EigenLU(const Eigen::SparseMatrix<double> A_matrix,
   return solver.solve(b_vector);
 }
 
-auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix<double> A_matrix,
+auto tug::solver::ThomasAlgorithm(const Eigen::SparseMatrix<double> &A_matrix,
-                                  const Eigen::VectorXd b_vector)
+                                  const Eigen::VectorXd &b_vector)
     -> Eigen::VectorXd {
   uint32_t n = b_vector.size();
 

test/testDiffusion.cpp

@@ -15,13 +15,12 @@ using namespace tug::diffusion;
 
 static std::vector<double> alpha(N *M, 1e-3);
 
-static TugInput setupDiffu(BoundaryCondition &bc) {
+static TugInput setupDiffu() {
   TugInput diffu;
 
   diffu.setTimestep(1);
   diffu.setGridCellN(N, M);
   diffu.setDomainSize(N, M);
-  diffu.setBoundaryCondition(bc);
 
   return diffu;
 }
@@ -33,7 +32,7 @@ TEST_CASE("closed boundaries - 1 concentration to 1 - rest 0") {
 
   BoundaryCondition bc(N, M);
 
-  TugInput diffu = setupDiffu(bc);
+  TugInput diffu = setupDiffu();
 
   uint32_t iterations = 1000;
   double sum = 0;
@@ -68,7 +67,8 @@ TEST_CASE("constant boundaries (0) - 1 concentration to 1 - rest 0") {
   bc.setSide(BC_SIDE_TOP, input);
   bc.setSide(BC_SIDE_BOTTOM, input);
 
-  TugInput diffu = setupDiffu(bc);
+  TugInput diffu = setupDiffu();
+  diffu.setBoundaryCondition(bc);
 
   uint32_t max_iterations = 20000;
   bool reached = false;
@@ -106,7 +106,8 @@ TEST_CASE(
   bc.setSide(BC_SIDE_TOP, top);
   bc.setSide(BC_SIDE_BOTTOM, bottom);
 
-  TugInput diffu = setupDiffu(bc);
+  TugInput diffu = setupDiffu();
+  diffu.setBoundaryCondition(bc);
 
   uint32_t max_iterations = 100;
 
@@ -136,7 +137,8 @@ TEST_CASE("2D closed boundaries, 1 constant cell in the middle") {
   field[MID] = val;
   bc(BC_INNER, MID) = {BC_TYPE_CONSTANT, val};
 
-  TugInput diffu = setupDiffu(bc);
+  TugInput diffu = setupDiffu();
+  diffu.setBoundaryCondition(bc);
 
   uint32_t max_iterations = 100;