feat: Add support for setting concentrations from a pointer

refactor: Use Row-major matrix internally

include/tug/Core/BTCS.hpp

@@ -51,7 +51,7 @@ constexpr std::pair<T, T> calcBoundaryCoeffClosed(T alpha_center, T alpha_side,
 // creates coefficient matrix for next time step from alphas in x-direction
 template <class T>
 static Eigen::SparseMatrix<T>
-createCoeffMatrix(const Eigen::MatrixX<T> &alpha,
+createCoeffMatrix(const typename Grid<T>::RowMajMat &alpha,
                   const std::vector<BoundaryElement<T>> &bcLeft,
                   const std::vector<BoundaryElement<T>> &bcRight,
                   const std::vector<std::pair<bool, T>> &inner_bc, int numCols,
@@ -160,9 +160,9 @@ constexpr T calcExplicitConcentrationsBoundaryConstant(T conc_center, T conc_bc,
 // concentrations
 template <class T>
 static Eigen::VectorX<T>
-createSolutionVector(const Eigen::MatrixX<T> &concentrations,
+createSolutionVector(const typename Grid<T>::RowMajMat &concentrations,
-                     const Eigen::MatrixX<T> &alphaX,
+                     const typename Grid<T>::RowMajMat &alphaX,
-                     const Eigen::MatrixX<T> &alphaY,
+                     const typename Grid<T>::RowMajMat &alphaY,
                      const std::vector<BoundaryElement<T>> &bcLeft,
                      const std::vector<BoundaryElement<T>> &bcRight,
                      const std::vector<BoundaryElement<T>> &bcTop,

include/tug/Grid.hpp

@@ -10,6 +10,8 @@
 
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <Eigen/src/Core/Matrix.h>
+#include <Eigen/src/Core/util/Constants.h>
 #include <stdexcept>
 
 namespace tug {
@@ -22,6 +24,9 @@ namespace tug {
  */
 template <class T> class Grid {
 public:
+  using RowMajMat =
+      Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
   /**
    * @brief Constructs a new 1D-Grid object of a given length, which holds a
    * matrix with concentrations and a respective matrix of alpha coefficients.
@@ -68,9 +73,9 @@ public:
     this->deltaCol =
         static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
 
-    this->concentrations = Eigen::MatrixX<T>::Constant(row, col, MAT_INIT_VAL);
+    this->concentrations = RowMajMat::Constant(row, col, MAT_INIT_VAL);
-    this->alphaX = Eigen::MatrixX<T>::Constant(row, col, MAT_INIT_VAL);
+    this->alphaX = RowMajMat::Constant(row, col, MAT_INIT_VAL);
-    this->alphaY = Eigen::MatrixX<T>::Constant(row, col, MAT_INIT_VAL);
+    this->alphaY = RowMajMat::Constant(row, col, MAT_INIT_VAL);
   }
 
   /**
@@ -90,13 +95,27 @@ public:
     this->concentrations = concentrations;
   }
 
+  /**
+   * @brief Sets the concentrations matrix for a 1D or 2D-Grid.
+   *
+   * @param concentrations A pointer to an array holding the concentrations.
+   * Array must have correct dimensions as defined in row and col. (Or length,
+   * in 1D case). There is no check for correct dimensions, so be careful!
+   */
+  void setConcentrations(T *concentrations) {
+    Eigen::Map<
+        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+        map(concentrations, this->row, this->col);
+    this->concentrations = map;
+  }
+
   /**
    * @brief Gets the concentrations matrix for a Grid.
    *
    * @return An Eigen3 matrix holding the concentrations and having
    * the same dimensions as the grid.
    */
-  const Eigen::MatrixX<T> &getConcentrations() { return this->concentrations; }
+  const auto &getConcentrations() { return this->concentrations; }
 
   /**
    * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
@@ -119,6 +138,26 @@ public:
     this->alphaX = alpha;
   }
 
+  /**
+   * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
+   * dimensional.
+   *
+   * @param alpha A pointer to an array holding the alpha coefficients. Array
+   * must have correct dimensions as defined in length. There is no check for
+   * correct dimensions, so be careful!
+   */
+  void setAlpha(T *alpha) {
+    if (dim != 1) {
+      throw std::invalid_argument(
+          "Grid is not one dimensional, you should probably "
+          "use 2D setter function!");
+    }
+    Eigen::Map<
+        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+        map(alpha, 1, this->col);
+    this->alphaX = map;
+  }
+
   /**
    * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
    * dimensional.
@@ -150,13 +189,40 @@ public:
     this->alphaY = alphaY;
   }
 
+  /**
+   * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
+   * dimensional.
+   *
+   * @param alphaX A pointer to an array holding the alpha coefficients in
+   * x-direction. Array must have correct dimensions as defined in row and col.
+   * There is no check for correct dimensions, so be careful!
+   * @param alphaY A pointer to an array holding the alpha coefficients in
+   * y-direction. Array must have correct dimensions as defined in row and col.
+   * There is no check for correct dimensions, so be careful!
+   */
+  void setAlpha(T *alphaX, T *alphaY) {
+    if (dim != 2) {
+      throw std::invalid_argument(
+          "Grid is not two dimensional, you should probably "
+          "use 1D setter function!");
+    }
+    Eigen::Map<
+        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+        mapX(alphaX, this->row, this->col);
+    Eigen::Map<
+        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+        mapY(alphaY, this->row, this->col);
+    this->alphaX = mapX;
+    this->alphaY = mapY;
+  }
+
   /**
    * @brief Gets the matrix of alpha coefficients of a 1D-Grid. Grid must be one
    * dimensional.
    *
    * @return A matrix with 1 row holding the alpha coefficients.
    */
-  const Eigen::MatrixX<T> &getAlpha() const {
+  const auto &getAlpha() const {
     if (dim != 1) {
       throw std::invalid_argument(
           "Grid is not one dimensional, you should probably "
@@ -172,7 +238,7 @@ public:
    *
    * @return A matrix holding the alpha coefficients in x-direction.
    */
-  const Eigen::MatrixX<T> &getAlphaX() const {
+  const auto &getAlphaX() const {
 
     if (dim != 2) {
       throw std::invalid_argument(
@@ -188,7 +254,7 @@ public:
    *
    * @return A matrix holding the alpha coefficients in y-direction.
    */
-  const Eigen::MatrixX<T> &getAlphaY() const {
+  const auto &getAlphaY() const {
 
     if (dim != 2) {
       throw std::invalid_argument(
@@ -316,16 +382,19 @@ public:
   }
 
 private:
-  int col;                          // number of grid columns
+  int col;        // number of grid columns
-  int row{1};                       // number of grid rows
+  int row{1};     // number of grid rows
-  int dim;                          // 1D or 2D
+  int dim;        // 1D or 2D
-  T domainCol;                      // number of domain columns
+  T domainCol;    // number of domain columns
-  T domainRow{0};                   // number of domain rows
+  T domainRow{0}; // number of domain rows
-  T deltaCol;                       // delta in x-direction (between columns)
+  T deltaCol;     // delta in x-direction (between columns)
-  T deltaRow{0};                    // delta in y-direction (between rows)
+  T deltaRow{0};  // delta in y-direction (between rows)
-  Eigen::MatrixX<T> concentrations; // Matrix holding grid concentrations
+
-  Eigen::MatrixX<T> alphaX; // Matrix holding alpha coefficients in x-direction
+  RowMajMat concentrations; // Matrix holding grid concentrations
-  Eigen::MatrixX<T> alphaY; // Matrix holding alpha coefficients in y-direction
+  RowMajMat alphaX;         // Matrix holding alpha coefficients in x-direction
+  RowMajMat alphaY;         // Matrix holding alpha coefficients in y-direction
+
+  using RowMajMatMap = Eigen::Map<RowMajMat>;
 
   static constexpr T MAT_INIT_VAL = 0;
 };

test/testGrid.cpp

@@ -1,6 +1,7 @@
 #include <Eigen/Core>
 #include <doctest/doctest.h>
 #include <tug/Grid.hpp>
+#include <vector>
 
 using namespace Eigen;
 using namespace std;
@@ -250,4 +251,18 @@ TEST_CASE("2D Grid64 non-quadratic") {
     dr = -2;
     CHECK_THROWS(grid.setDomain(dr, dc));
   }
+
+  SUBCASE("set concentration from pointer") {
+    std::vector<double> concentrations(r * c);
+
+    for (int i = 0; i < r * c; i++) {
+      concentrations[i] = i;
+    }
+
+    grid.setConcentrations(concentrations.data());
+
+    CHECK_EQ(grid.getConcentrations()(0, 0), 0);
+    CHECK_EQ(grid.getConcentrations()(0, 1), 1);
+    CHECK_EQ(grid.getConcentrations()(1, 0), c);
+  }
 }