Merge branch 'ml/row-major-mat' into 'main'

feat: Add support for setting concentrations from a pointer

See merge request naaice/tug!32

include/tug/Core/BTCS.hpp

@@ -10,6 +10,7 @@
 #ifndef BTCS_H_
 #define BTCS_H_
 
+#include "Matrix.hpp"
 #include "TugUtils.hpp"
 
 #include <cstddef>
@@ -51,7 +52,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 RowMajMat<T> &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 +161,8 @@ 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 RowMajMat<T> &concentrations,
-                     const Eigen::MatrixX<T> &alphaX,
+                     const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY,
-                     const Eigen::MatrixX<T> &alphaY,
                      const std::vector<BoundaryElement<T>> &bcLeft,
                      const std::vector<BoundaryElement<T>> &bcRight,
                      const std::vector<BoundaryElement<T>> &bcTop,
@@ -369,7 +369,7 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
 
   const auto inner_bc = bc.getInnerBoundaryRow(0);
 
-  Eigen::MatrixX<T> concentrations = grid.getConcentrations();
+  RowMajMat<T> &concentrations = grid.getConcentrations();
   int rowIndex = 0;
   A = createCoeffMatrix(alpha, bcLeft, bcRight, inner_bc, length, rowIndex,
                         sx); // this is exactly same as in 2D
@@ -385,13 +385,13 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
     b(length - 1) += 2 * sx * alpha(0, length - 1) * bcRight[0].getValue();
   }
 
-  concentrations_t1 = solverFunc(A, b);
+  concentrations = solverFunc(A, b);
 
-  for (int j = 0; j < length; j++) {
+  // for (int j = 0; j < length; j++) {
-    concentrations(0, j) = concentrations_t1(j);
+  //   concentrations(0, j) = concentrations_t1(j);
-  }
+  // }
 
-  grid.setConcentrations(concentrations);
+  // grid.setConcentrations(concentrations);
 }
 
 // BTCS solution for 2D grid
@@ -405,24 +405,22 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
   T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
   T sy = timestep / (2 * grid.getDeltaRow() * grid.getDeltaRow());
 
-  Eigen::MatrixX<T> concentrations_t1 =
+  RowMajMat<T> concentrations_t1(rowMax, colMax);
-      Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
-  Eigen::VectorX<T> row_t1(colMax);
 
   Eigen::SparseMatrix<T> A;
   Eigen::VectorX<T> b;
 
-  auto alphaX = grid.getAlphaX();
+  RowMajMat<T> alphaX = grid.getAlphaX();
-  auto alphaY = grid.getAlphaY();
+  RowMajMat<T> alphaY = grid.getAlphaY();
 
   const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT);
   const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT);
   const auto &bcTop = bc.getBoundarySide(BC_SIDE_TOP);
   const auto &bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM);
 
-  Eigen::MatrixX<T> concentrations = grid.getConcentrations();
+  RowMajMat<T> &concentrations = grid.getConcentrations();
 
-#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
+#pragma omp parallel for num_threads(numThreads) private(A, b)
   for (int i = 0; i < rowMax; i++) {
     auto inner_bc = bc.getInnerBoundaryRow(i);
 
@@ -430,17 +428,14 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
     b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
                              bcTop, bcBottom, inner_bc, colMax, i, sx, sy);
 
-    row_t1 = solverFunc(A, b);
+    concentrations_t1.row(i) = solverFunc(A, b);
-
-    concentrations_t1.row(i) = row_t1;
   }
 
   concentrations_t1.transposeInPlace();
-  concentrations.transposeInPlace();
   alphaX.transposeInPlace();
   alphaY.transposeInPlace();
 
-#pragma omp parallel for num_threads(numThreads) private(A, b, row_t1)
+#pragma omp parallel for num_threads(numThreads) private(A, b)
   for (int i = 0; i < colMax; i++) {
     auto inner_bc = bc.getInnerBoundaryCol(i);
     // swap alphas, boundary conditions and sx/sy for column-wise calculation
@@ -448,14 +443,8 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
     b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
                              bcLeft, bcRight, inner_bc, rowMax, i, sy, sx);
 
-    row_t1 = solverFunc(A, b);
+    concentrations.col(i) = solverFunc(A, b);
-
-    concentrations.row(i) = row_t1;
   }
-
-  concentrations.transposeInPlace();
-
-  grid.setConcentrations(concentrations);
 }
 
 // entry point for EigenLU solver; differentiate between 1D and 2D grid

include/tug/Core/FTCS.hpp

@@ -228,8 +228,7 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
   T deltaCol = grid.getDeltaCol();
 
   // matrix for concentrations at time t+1
-  Eigen::MatrixX<T> concentrations_t1 =
+  RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(1, colMax, 0);
-      Eigen::MatrixX<T>::Constant(1, colMax, 0);
 
   // only one row in 1D case -> row constant at index 0
   int row = 0;
@@ -270,8 +269,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
   T deltaCol = grid.getDeltaCol();
 
   // matrix for concentrations at time t+1
-  Eigen::MatrixX<T> concentrations_t1 =
+  RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
-      Eigen::MatrixX<T>::Constant(rowMax, colMax, 0);
 
   // inner cells
   // these are independent of the boundary condition type

include/tug/Core/Matrix.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <Eigen/Core>
+
+namespace tug {
+/**
+ * @brief Alias template for a row-major matrix using Eigen library.
+ *
+ * This alias template defines a type `RowMajMat` which represents a row-major
+ * matrix using the Eigen library. It is a template that takes a type `T` as its
+ * template parameter. The matrix is dynamically sized with `Eigen::Dynamic` for
+ * both rows and columns. The matrix is stored in row-major order.
+ *
+ * @tparam T The type of the matrix elements.
+ */
+template <typename T>
+using RowMajMat =
+    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+template <typename T> using RowMajMatMap = Eigen::Map<RowMajMat<T>>;
+} // namespace tug

include/tug/Grid.hpp

@@ -8,8 +8,11 @@
  *
  */
 
+#include "Core/Matrix.hpp"
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <Eigen/src/Core/Matrix.h>
+#include <Eigen/src/Core/util/Constants.h>
 #include <stdexcept>
 
 namespace tug {
@@ -40,8 +43,8 @@ public:
     this->deltaCol =
         static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
 
-    this->concentrations = Eigen::MatrixX<T>::Constant(1, col, MAT_INIT_VAL);
+    this->concentrations = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
-    this->alphaX = Eigen::MatrixX<T>::Constant(1, col, MAT_INIT_VAL);
+    this->alphaX = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
   }
 
   /**
@@ -68,9 +71,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<T>::Constant(row, col, MAT_INIT_VAL);
-    this->alphaX = Eigen::MatrixX<T>::Constant(row, col, MAT_INIT_VAL);
+    this->alphaX = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
-    this->alphaY = Eigen::MatrixX<T>::Constant(row, col, MAT_INIT_VAL);
+    this->alphaY = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
   }
 
   /**
@@ -80,7 +83,7 @@ public:
    * Matrix must have correct dimensions as defined in row and col. (Or length,
    * in 1D case).
    */
-  void setConcentrations(const Eigen::MatrixX<T> &concentrations) {
+  void setConcentrations(const RowMajMat<T> &concentrations) {
     if (concentrations.rows() != this->row ||
         concentrations.cols() != this->col) {
       throw std::invalid_argument(
@@ -90,13 +93,25 @@ 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) {
+    tug::RowMajMatMap<T> 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; }
+  auto &getConcentrations() { return this->concentrations; }
 
   /**
    * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
@@ -105,7 +120,7 @@ public:
    * @param alpha An Eigen3 MatrixX<T> with 1 row holding the alpha
    * coefficients. Matrix columns must have same size as length of grid.
    */
-  void setAlpha(const Eigen::MatrixX<T> &alpha) {
+  void setAlpha(const RowMajMat<T> &alpha) {
     if (dim != 1) {
       throw std::invalid_argument(
           "Grid is not one dimensional, you should probably "
@@ -119,6 +134,24 @@ 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!");
+    }
+    RowMajMatMap<T> map(alpha, 1, this->col);
+    this->alphaX = map;
+  }
+
   /**
    * @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
    * dimensional.
@@ -128,8 +161,7 @@ public:
    * @param alphaY An Eigen3 MatrixX<T> holding the alpha coefficients in
    * y-direction. Matrix must be of same size as the grid.
    */
-  void setAlpha(const Eigen::MatrixX<T> &alphaX,
+  void setAlpha(const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY) {
-                const Eigen::MatrixX<T> &alphaY) {
     if (dim != 2) {
       throw std::invalid_argument(
           "Grid is not two dimensional, you should probably "
@@ -150,13 +182,36 @@ 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!");
+    }
+    RowMajMatMap<T> mapX(alphaX, this->row, this->col);
+    RowMajMatMap<T> 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 +227,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 +243,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 +371,17 @@ 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<T> concentrations; // Matrix holding grid concentrations
-  Eigen::MatrixX<T> alphaY; // Matrix holding alpha coefficients in y-direction
+  RowMajMat<T> alphaX; // Matrix holding alpha coefficients in x-direction
+  RowMajMat<T> alphaY; // Matrix holding alpha coefficients in y-direction
 
   static constexpr T MAT_INIT_VAL = 0;
 };

include/tug/Simulation.hpp

@@ -446,9 +446,9 @@ public:
       // TODO this implementation is very inefficient!
       // a separate implementation that sets up a specific tridiagonal matrix
       // for Crank-Nicolson would be better
-      Eigen::MatrixX<T> concentrations;
+      RowMajMat<T> concentrations;
-      Eigen::MatrixX<T> concentrationsFTCS;
+      RowMajMat<T> concentrationsFTCS;
-      Eigen::MatrixX<T> concentrationsResult;
+      RowMajMat<T> concentrationsResult;
       for (int i = 0; i < iterations * innerIterations; i++) {
         if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
           printConcentrationsConsole();

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);
+  }
 }