From d3843fb2a334a86c352f10f4c3c978792cfb2b53 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Thu, 13 Jun 2024 09:16:51 +0200 Subject: [PATCH] refactor: Update testDiffusion.cpp and Diffusion.hpp Refactor testDiffusion.cpp and Diffusion.hpp to improve code readability and maintainability. Remove unnecessary exception throwing and replace with assert statements for invalid arguments. --- .gitlab-ci.yml | 2 +- examples/BTCS_2D_proto_example.cpp | 3 +- examples/FTCS_2D_proto_closed_mdl.cpp | 3 +- examples/FTCS_2D_proto_example_mdl.cpp | 3 +- include/tug/Boundary.hpp | 126 ++++------ include/tug/Core/Numeric/BTCS.hpp | 6 +- include/tug/Core/Numeric/FTCS.hpp | 35 +-- include/tug/Diffusion.hpp | 16 +- include/tug/Grid.hpp | 324 +++++++++++-------------- naaice/BTCS_2D_NAAICE.cpp | 5 +- naaice/NAAICE_dble_vs_float.cpp | 5 +- test/CMakeLists.txt | 3 +- test/setup.cpp | 9 +- test/testBoundary.cpp | 31 ++- test/testDiffusion.cpp | 28 +-- test/testGrid.cpp | 214 ++++++++-------- 16 files changed, 368 insertions(+), 445 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b7ed776..4a34f82 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,7 +12,7 @@ test: stage: test script: - mkdir build && cd build - - cmake -DCMAKE_BUILD_TYPE=Release -DTUG_ENABLE_TESTING=ON -G Ninja .. + - cmake -DCMAKE_BUILD_TYPE=Debug -DTUG_ENABLE_TESTING=ON -G Ninja .. - ninja - ctest --output-junit test_results.xml artifacts: diff --git a/examples/BTCS_2D_proto_example.cpp b/examples/BTCS_2D_proto_example.cpp index 9cefd10..c09d93d 100644 --- a/examples/BTCS_2D_proto_example.cpp +++ b/examples/BTCS_2D_proto_example.cpp @@ -14,7 +14,6 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int row = 40; int col = 50; - Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -24,7 +23,7 @@ int main(int argc, char *argv[]) { // #row,#col,value grid.setConcentrations(concentrations); MatrixXd concentrations = MatrixXd::Constant(row, col, 0); concentrations(10, 10) = 2000; - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: // MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value diff --git a/examples/FTCS_2D_proto_closed_mdl.cpp b/examples/FTCS_2D_proto_closed_mdl.cpp index a5bd8f4..69b31c3 100644 --- a/examples/FTCS_2D_proto_closed_mdl.cpp +++ b/examples/FTCS_2D_proto_closed_mdl.cpp @@ -31,7 +31,6 @@ int main(int argc, char *argv[]) { // create a grid with a 20 x 20 field int n2 = row / 2 - 1; - Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -44,7 +43,7 @@ int main(int argc, char *argv[]) { concentrations(n2, n2 + 1) = 1; concentrations(n2 + 1, n2) = 1; concentrations(n2 + 1, n2 + 1) = 1; - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value diff --git a/examples/FTCS_2D_proto_example_mdl.cpp b/examples/FTCS_2D_proto_example_mdl.cpp index 3f750eb..41ddb91 100644 --- a/examples/FTCS_2D_proto_example_mdl.cpp +++ b/examples/FTCS_2D_proto_example_mdl.cpp @@ -22,7 +22,6 @@ int main(int argc, char *argv[]) { int row = 64; int col = 64; int n2 = row / 2 - 1; - Grid64 grid(row, col); // (optional) set the domain, e.g.: // grid.setDomain(20, 20); @@ -35,7 +34,7 @@ int main(int argc, char *argv[]) { concentrations(n2, n2 + 1) = 1; concentrations(n2 + 1, n2) = 1; concentrations(n2 + 1, n2 + 1) = 1; - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); // (optional) set alphas of the grid, e.g.: MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value diff --git a/include/tug/Boundary.hpp b/include/tug/Boundary.hpp index ebcde98..cdfe3be 100644 --- a/include/tug/Boundary.hpp +++ b/include/tug/Boundary.hpp @@ -8,6 +8,7 @@ #define BOUNDARY_H_ #include "Grid.hpp" +#include "tug/Core/TugUtils.hpp" #include #include @@ -161,13 +162,11 @@ public: * @param side Side to be set to closed, e.g. BC_SIDE_LEFT. */ void setBoundarySideClosed(BC_SIDE side) { - if (this->dim == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw std::invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } + tug_assert((this->dim > 1) || + ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)), + "For the " + "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT " + "borders exist."); const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; const int n = is_vertical ? this->rows : this->cols; @@ -186,13 +185,11 @@ public: * page. */ void setBoundarySideConstant(BC_SIDE side, double value) { - if (this->dim == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw std::invalid_argument( - "For the one-dimensional case, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } + tug_assert((this->dim > 1) || + ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)), + "For the " + "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT " + "borders exist."); const bool is_vertical = side == BC_SIDE_LEFT || side == BC_SIDE_RIGHT; const int n = is_vertical ? this->rows : this->cols; @@ -212,10 +209,9 @@ public: */ void setBoundaryElemenClosed(BC_SIDE side, int index) { // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument( - "Index is selected either too large or too small."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); + this->boundaries[side][index].setType(BC_TYPE_CLOSED); } @@ -233,10 +229,8 @@ public: */ void setBoundaryElementConstant(BC_SIDE side, int index, double value) { // tests whether the index really points to an element of the boundary side. - if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument( - "Index is selected either too large or too small."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); this->boundaries[side][index].setType(BC_TYPE_CONSTANT); this->boundaries[side][index].setValue(value); } @@ -251,13 +245,11 @@ public: * BoundaryElement objects. */ const std::vector> &getBoundarySide(BC_SIDE side) const { - if (this->dim == 1) { - if ((side == BC_SIDE_BOTTOM) || (side == BC_SIDE_TOP)) { - throw std::invalid_argument( - "For the one-dimensional trap, only the BC_SIDE_LEFT and " - "BC_SIDE_RIGHT borders exist."); - } - } + tug_assert((this->dim > 1) || + ((side == BC_SIDE_LEFT) || (side == BC_SIDE_RIGHT)), + "For the " + "one-dimensional case, only the BC_SIDE_LEFT and BC_SIDE_RIGHT " + "borders exist."); return this->boundaries[side]; } @@ -295,10 +287,8 @@ public: * object. */ BoundaryElement getBoundaryElement(BC_SIDE side, int index) const { - if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument( - "Index is selected either too large or too small."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); return this->boundaries[side][index]; } @@ -313,10 +303,8 @@ public: * @return Boundary Type of the corresponding boundary condition. */ BC_TYPE getBoundaryElementType(BC_SIDE side, int index) const { - if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument( - "Index is selected either too large or too small."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); return this->boundaries[side][index].getType(); } @@ -333,14 +321,12 @@ public: * object. */ T getBoundaryElementValue(BC_SIDE side, int index) const { - if ((boundaries[side].size() < index) || index < 0) { - throw std::invalid_argument( - "Index is selected either too large or too small."); - } - if (boundaries[side][index].getType() != BC_TYPE_CONSTANT) { - throw std::invalid_argument( - "A value can only be output if it is a constant boundary condition."); - } + tug_assert(boundaries[side].size() > index && index >= 0, + "Index is selected either too large or too small."); + tug_assert( + boundaries[side][index].getType() == BC_TYPE_CONSTANT, + "A value can only be output if it is a constant boundary condition."); + return this->boundaries[side][index].getValue(); } @@ -382,13 +368,8 @@ public: * @param value Value of the inner constant boundary condition */ void setInnerBoundary(std::uint32_t index, T value) { - if (this->dim != 1) { - throw std::invalid_argument( - "This function is only available for 1D grids."); - } - if (index >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 1, "This function is only available for 1D grids."); + tug_assert(index < this->cols, "Index is out of bounds."); this->inner_boundary[std::make_pair(0, index)] = value; } @@ -401,13 +382,8 @@ public: * @param value Value of the inner constant boundary condition */ void setInnerBoundary(std::uint32_t row, std::uint32_t col, T value) { - if (this->dim != 2) { - throw std::invalid_argument( - "This function is only available for 2D grids."); - } - if (row >= this->rows || col >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 2, "This function is only available for 2D grids."); + tug_assert(row < this->rows && col < this->cols, "Index is out of bounds."); this->inner_boundary[std::make_pair(row, col)] = value; } @@ -420,13 +396,8 @@ public: * set or not) and value of the inner constant boundary condition */ std::pair getInnerBoundary(std::uint32_t index) const { - if (this->dim != 1) { - throw std::invalid_argument( - "This function is only available for 1D grids."); - } - if (index >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 1, "This function is only available for 1D grids."); + tug_assert(index < this->cols, "Index is out of bounds."); auto it = this->inner_boundary.find(std::make_pair(0, index)); if (it == this->inner_boundary.end()) { @@ -445,13 +416,8 @@ public: */ std::pair getInnerBoundary(std::uint32_t row, std::uint32_t col) const { - if (this->dim != 2) { - throw std::invalid_argument( - "This function is only available for 2D grids."); - } - if (row >= this->rows || col >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 2, "This function is only available for 2D grids."); + tug_assert(row < this->rows && col < this->cols, "Index is out of bounds."); auto it = this->inner_boundary.find(std::make_pair(row, col)); if (it == this->inner_boundary.end()) { @@ -471,9 +437,7 @@ public: * condition */ std::vector> getInnerBoundaryRow(std::uint32_t row) const { - if (row >= this->rows) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(row < this->rows, "Index is out of bounds."); if (inner_boundary.empty()) { return std::vector>(this->cols, @@ -499,14 +463,8 @@ public: * condition */ std::vector> getInnerBoundaryCol(std::uint32_t col) const { - if (this->dim != 2) { - throw std::invalid_argument( - "This function is only available for 2D grids."); - } - - if (col >= this->cols) { - throw std::invalid_argument("Index is out of bounds."); - } + tug_assert(this->dim == 2, "This function is only available for 2D grids."); + tug_assert(col < this->cols, "Index is out of bounds."); if (inner_boundary.empty()) { return std::vector>(this->rows, diff --git a/include/tug/Core/Numeric/BTCS.hpp b/include/tug/Core/Numeric/BTCS.hpp index 3aca4a8..460119d 100644 --- a/include/tug/Core/Numeric/BTCS.hpp +++ b/include/tug/Core/Numeric/BTCS.hpp @@ -354,15 +354,15 @@ template static void BTCS_1D(Grid &grid, Boundary &bc, T timestep, Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, Eigen::VectorX &b)) { - int length = grid.getLength(); - T sx = timestep / (grid.getDelta() * grid.getDelta()); + int length = grid.getCol(); + T sx = timestep / (grid.getDeltaCol() * grid.getDeltaCol()); Eigen::VectorX concentrations_t1(length); Eigen::SparseMatrix A; Eigen::VectorX b(length); - const auto &alpha = grid.getAlpha(); + const auto &alpha = grid.getAlphaX(); const auto &bcLeft = bc.getBoundarySide(BC_SIDE_LEFT); const auto &bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); diff --git a/include/tug/Core/Numeric/FTCS.hpp b/include/tug/Core/Numeric/FTCS.hpp index 41eeb80..9e899aa 100644 --- a/include/tug/Core/Numeric/FTCS.hpp +++ b/include/tug/Core/Numeric/FTCS.hpp @@ -9,7 +9,9 @@ #define FTCS_H_ #include "../TugUtils.hpp" +#include "tug/Core/Matrix.hpp" +#include #include #ifdef _OPENMP @@ -225,6 +227,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { int colMax = grid.getCol(); T deltaCol = grid.getDeltaCol(); + RowMajMat &concentrations_grid = grid.getConcentrations(); // matrix for concentrations at time t+1 RowMajMat concentrations_t1 = RowMajMat::Constant(1, colMax, 0); @@ -234,7 +237,7 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { // inner cells // independent of boundary condition type for (int col = 1; col < colMax - 1; col++) { - concentrations_t1(row, col) = grid.getConcentrations()(row, col) + + concentrations_t1(row, col) = concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChange(grid, row, col)); } @@ -242,19 +245,20 @@ static void FTCS_1D(Grid &grid, Boundary &bc, T timestep) { // left boundary; hold column constant at index 0 int col = 0; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeLeftBoundary(grid, bc, row, col)); // right boundary; hold column constant at max index col = colMax - 1; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeRightBoundary(grid, bc, row, col)); // overwrite obsolete concentrations - grid.setConcentrations(concentrations_t1); + std::memcpy(concentrations_grid.data(), concentrations_t1.data(), + colMax * sizeof(T)); } // FTCS solution for 2D grid @@ -266,6 +270,8 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, T deltaRow = grid.getDeltaRow(); T deltaCol = grid.getDeltaCol(); + RowMajMat &concentrations_grid = grid.getConcentrations(); + // matrix for concentrations at time t+1 RowMajMat concentrations_t1 = RowMajMat::Constant(rowMax, colMax, 0); @@ -275,7 +281,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int row = 1; row < rowMax - 1; row++) { for (int col = 1; col < colMax - 1; col++) { - concentrations_t1(row, col) = grid.getConcentrations()(row, col) + + concentrations_t1(row, col) = concentrations_grid(row, col) + timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)) + timestep / (deltaCol * deltaCol) * @@ -290,7 +296,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int row = 1; row < rowMax - 1; row++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); @@ -302,7 +308,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int row = 1; row < rowMax - 1; row++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * (calcVerticalChange(grid, row, col)); @@ -314,7 +320,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int col = 1; col < colMax - 1; col++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaRow * deltaRow) * (calcVerticalChangeTopBoundary(grid, bc, row, col)) + timestep / (deltaCol * deltaCol) * @@ -327,7 +333,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, #pragma omp parallel for num_threads(numThreads) for (int col = 1; col < colMax - 1; col++) { concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaRow * deltaRow) * (calcVerticalChangeBottomBoundary(grid, bc, row, col)) + timestep / (deltaCol * deltaCol) * @@ -339,7 +345,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = 0; col = 0; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * @@ -350,7 +356,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = 0; col = colMax - 1; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * @@ -361,7 +367,7 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = rowMax - 1; col = 0; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeLeftBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * @@ -372,14 +378,15 @@ static void FTCS_2D(Grid &grid, Boundary &bc, T timestep, row = rowMax - 1; col = colMax - 1; concentrations_t1(row, col) = - grid.getConcentrations()(row, col) + + concentrations_grid(row, col) + timestep / (deltaCol * deltaCol) * (calcHorizontalChangeRightBoundary(grid, bc, row, col)) + timestep / (deltaRow * deltaRow) * (calcVerticalChangeBottomBoundary(grid, bc, row, col)); // overwrite obsolete concentrations - grid.setConcentrations(concentrations_t1); + std::memcpy(concentrations_grid.data(), concentrations_t1.data(), + rowMax * colMax * sizeof(T)); // } } diff --git a/include/tug/Diffusion.hpp b/include/tug/Diffusion.hpp index 4d2f0f0..675605c 100644 --- a/include/tug/Diffusion.hpp +++ b/include/tug/Diffusion.hpp @@ -96,17 +96,15 @@ public: * @param timestep Valid timestep greater than zero. */ void setTimestep(T timestep) { - if (timestep <= 0) { - throw_invalid_argument("Timestep has to be greater than zero."); - } + tug_assert(timestep > 0, "Timestep has to be greater than zero."); if constexpr (approach == FTCS_APPROACH || approach == CRANK_NICOLSON_APPROACH) { T cfl; if (grid.getDim() == 1) { - const T deltaSquare = grid.getDelta(); - const T maxAlpha = grid.getAlpha().maxCoeff(); + const T deltaSquare = grid.getDeltaCol(); + const T maxAlpha = grid.getAlphaX().maxCoeff(); // Courant-Friedrichs-Lewy condition cfl = deltaSquare / (4 * maxAlpha); @@ -272,12 +270,8 @@ public: * parameters. */ void run() { - if (this->timestep == -1) { - throw_invalid_argument("Timestep is not set!"); - } - if (this->iterations == -1) { - throw_invalid_argument("Number of iterations are not set!"); - } + tug_assert(this->timestep > 0, "Timestep is not set!"); + tug_assert(this->iterations > 0, "Number of iterations are not set!"); std::string filename; if (this->console_output > CONSOLE_OUTPUT_OFF) { diff --git a/include/tug/Grid.hpp b/include/tug/Grid.hpp index b8e471a..6e73742 100644 --- a/include/tug/Grid.hpp +++ b/include/tug/Grid.hpp @@ -1,5 +1,4 @@ -#ifndef GRID_H_ -#define GRID_H_ +#pragma once /** * @file Grid.hpp @@ -9,11 +8,12 @@ */ #include "Core/Matrix.hpp" +#include "tug/Core/TugUtils.hpp" #include #include #include #include -#include +#include namespace tug { @@ -26,83 +26,97 @@ namespace tug { template class Grid { public: /** - * @brief Constructs a new 1D-Grid object of a given length, which holds a - * matrix with concentrations and a respective matrix of alpha coefficients. - * The domain length is per default the same as the length. The - * concentrations are all 20 by default and the alpha coefficients are 1. + * @brief Construct a new Grid object. * - * @param length Length of the 1D-Grid. Must be greater than 3. + * Constructs a new Grid object with given concentrations, defined by an + * Eigen::Matrix. The domain length is per default the same as the length. The + * alpha coefficients are set to 1. The dimensions of the grid are determined + * by the given matrix, which can also be an Eigen::Vector for a 1D-Grid. + * + * @param concentrations An Eigen3 MatrixX holding the concentrations. */ - Grid(int length) : col(length), domainCol(length) { - if (length <= 3) { - throw std::invalid_argument( - "Given grid length too small. Must be greater than 3."); + Grid(const RowMajMat &concentrations) { + if (concentrations.rows() == 1) { + this->dim = 1; + this->domainCol = static_cast(concentrations.cols()); + this->deltaCol = static_cast(this->domainCol) / + static_cast(concentrations.cols()); // -> 1 + + this->concentrations = concentrations; + return; } - this->dim = 1; - this->deltaCol = - static_cast(this->domainCol) / static_cast(this->col); // -> 1 + if (concentrations.cols() == 1) { + this->dim = 1; + this->domainCol = static_cast(concentrations.rows()); + this->deltaCol = static_cast(this->domainCol) / + static_cast(concentrations.rows()); // -> 1 - this->concentrations = RowMajMat::Constant(1, col, MAT_INIT_VAL); - this->alphaX = RowMajMat::Constant(1, col, MAT_INIT_VAL); - } - - /** - * @brief Constructs a new 2D-Grid object of given dimensions, which holds a - * matrix with concentrations and the respective matrices of alpha coefficient - * for each direction. The domain in x- and y-direction is per default equal - * to the col length and row length, respectively. The concentrations are all - * 20 by default across the entire grid and the alpha coefficients 1 in both - * directions. - * - * @param row Length of the 2D-Grid in y-direction. Must be greater than 3. - * @param col Length of the 2D-Grid in x-direction. Must be greater than 3. - */ - Grid(int _row, int _col) - : row(_row), col(_col), domainRow(_row), domainCol(_col) { - if (row <= 1 || col <= 1) { - throw std::invalid_argument( - "At least one dimension is 1. Use 1D grid for better results."); + this->concentrations = concentrations.transpose(); + return; } this->dim = 2; - this->deltaRow = - static_cast(this->domainRow) / static_cast(this->row); // -> 1 - this->deltaCol = - static_cast(this->domainCol) / static_cast(this->col); // -> 1 - - this->concentrations = RowMajMat::Constant(row, col, MAT_INIT_VAL); - this->alphaX = RowMajMat::Constant(row, col, MAT_INIT_VAL); - this->alphaY = RowMajMat::Constant(row, col, MAT_INIT_VAL); - } - - /** - * @brief Sets the concentrations matrix for a 1D or 2D-Grid. - * - * @param concentrations An Eigen3 MatrixX holding the concentrations. - * Matrix must have correct dimensions as defined in row and col. (Or length, - * in 1D case). - */ - void setConcentrations(const RowMajMat &concentrations) { - if (concentrations.rows() != this->row || - concentrations.cols() != this->col) { - throw std::invalid_argument( - "Given matrix of concentrations mismatch with Grid dimensions!"); - } + this->domainRow = static_cast(concentrations.rows()); + this->domainCol = static_cast(concentrations.cols()); + this->deltaRow = static_cast(this->domainRow) / + static_cast(concentrations.rows()); // -> 1 + this->deltaCol = static_cast(this->domainCol) / + static_cast(concentrations.cols()); // -> 1 this->concentrations = concentrations; + // this->alphaX = RowMajMat::Constant(concentrations.rows(), + // concentrations.cols(), + // MAT_INIT_VAL); + // this->alphaY = RowMajMat::Constant(concentrations.rows(), + // concentrations.cols(), + // MAT_INIT_VAL); } /** - * @brief Sets the concentrations matrix for a 1D or 2D-Grid. + * @brief Construct a new Grid object. * - * @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! + * Constructs a new 1D Grid object with given concentrations, defined by a + * pointer to consecutive memory and the length of the array. The domain + * length is per default the same as the count of grid cells (length of + * array). The memory region is mapped internally, changes will affect the + * original array and the memory shall not be freed. There is no check for + * correct memory size! + * + * @param concentrations Pointer to consecutive memory holding concentrations. + * @param length Length of the array/the 1D grid. */ - void setConcentrations(T *concentrations) { - tug::RowMajMatMap map(concentrations, this->row, this->col); - this->concentrations = map; + Grid(T *concentrations, std::size_t length) : dim(1) { + this->domainCol = static_cast(length); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(length); // -> 1 + + this->concentrations = RowMajMatMap(concentrations, 1, length); + } + + /** + * @brief Construct a new Grid object. + * + * Constructs a new 2D Grid object with given concentrations, defined by a + * pointer to consecutive memory and the number of rows and columns. The + * domain size is per default the same as the number of rows and columns. The + * memory region is mapped internally, changes will affect the original array + * and the memory shall not be freed. There is no check for correct memory + * size! + * + * @param concentrations Pointer to consecutive memory holding concentrations. + * @param row Number of rows. + * @param col Number of columns. + */ + Grid(T *concentrations, std::size_t row, std::size_t col) : dim(2) { + this->domainRow = static_cast(row); // -> 1 + this->domainCol = static_cast(col); // -> 1 + this->deltaCol = + static_cast(this->domainCol) / static_cast(col); // -> 1 + this->deltaRow = + static_cast(this->domainRow) / static_cast(row); // -> 1 + + this->concentrations = RowMajMatMap(concentrations, row, col); } /** @@ -113,6 +127,17 @@ public: */ auto &getConcentrations() { return this->concentrations; } + void initAlpha() { + this->alphaX = RowMajMat::Constant( + this->concentrations.rows(), this->concentrations.cols(), MAT_INIT_VAL); + if (dim > 1) { + + this->alphaY = + RowMajMat::Constant(this->concentrations.rows(), + this->concentrations.cols(), MAT_INIT_VAL); + } + } + /** * @brief Set the alpha coefficients of a 1D-Grid. Grid must be one * dimensional. @@ -121,15 +146,12 @@ public: * coefficients. Matrix columns must have same size as length of grid. */ void setAlpha(const RowMajMat &alpha) { - if (dim != 1) { - throw std::invalid_argument( - "Grid is not one dimensional, you should probably " - "use 2D setter function!"); - } - if (alpha.rows() != 1 || alpha.cols() != this->col) { - throw std::invalid_argument( - "Given matrix of alpha coefficients mismatch with Grid dimensions!"); - } + tug_assert(dim == 1, + "Grid is not one dimensional, use 2D setter function!"); + + tug_assert( + alpha.rows() == 1 && alpha.cols() == this->concentrations.cols(), + "Given matrix of alpha coefficients mismatch with Grid dimensions!"); this->alphaX = alpha; } @@ -143,11 +165,9 @@ public: * 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!"); - } + tug_assert(dim == 1, + "Grid is not one dimensional, use 2D setter function!"); + RowMajMatMap map(alpha, 1, this->col); this->alphaX = map; } @@ -162,21 +182,21 @@ public: * y-direction. Matrix must be of same size as the grid. */ void setAlpha(const RowMajMat &alphaX, const RowMajMat &alphaY) { - if (dim != 2) { - throw std::invalid_argument( - "Grid is not two dimensional, you should probably " - "use 1D setter function!"); - } - if (alphaX.rows() != this->row || alphaX.cols() != this->col) { - throw std::invalid_argument( - "Given matrix of alpha coefficients in x-direction " - "mismatch with GRid dimensions!"); - } - if (alphaY.rows() != this->row || alphaY.cols() != this->col) { - throw std::invalid_argument( - "Given matrix of alpha coefficients in y-direction " - "mismatch with GRid dimensions!"); - } + tug_assert(dim == 2, + "Grid is not two dimensional, use 1D setter function!"); + + tug_assert(alphaX.rows() == this->concentrations.rows(), + "Alpha in x-direction " + "has wrong number of rows!"); + tug_assert(alphaX.cols() == this->concentrations.cols(), + "Alpha in x-direction has wrong number of columns!"); + + tug_assert(alphaY.rows() == this->concentrations.rows(), + "Alpha in y-direction " + "has wrong number of rows!"); + + tug_assert(alphaY.cols() == this->concentrations.cols(), + "Alpha in y-direction has wrong number of columns!"); this->alphaX = alphaX; this->alphaY = alphaY; @@ -194,33 +214,14 @@ public: * 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!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!"); + RowMajMatMap mapX(alphaX, this->row, this->col); RowMajMatMap 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 auto &getAlpha() const { - if (dim != 1) { - throw std::invalid_argument( - "Grid is not one dimensional, you should probably " - "use either getAlphaX() or getAlphaY()!"); - } - - return this->alphaX; - } - /** * @brief Gets the matrix of alpha coefficients in x-direction of a 2D-Grid. * Grid must be two dimensional. @@ -228,12 +229,7 @@ public: * @return A matrix holding the alpha coefficients in x-direction. */ const auto &getAlphaX() const { - - if (dim != 2) { - throw std::invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } - + tug_assert(this->alphaX.size() > 0, "AlphaX is empty!"); return this->alphaX; } @@ -244,11 +240,8 @@ public: * @return A matrix holding the alpha coefficients in y-direction. */ const auto &getAlphaY() const { - - if (dim != 2) { - throw std::invalid_argument( - "Grid is not two dimensional, you should probably use getAlpha()!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, there is no alphaY!"); + tug_assert(this->alphaY.size() > 0, "AlphaY is empty!"); return this->alphaY; } @@ -260,34 +253,19 @@ public: */ int getDim() const { return this->dim; } - /** - * @brief Gets length of 1D grid. Must be one dimensional grid. - * - * @return Length of 1D grid. - */ - int getLength() const { - if (dim != 1) { - throw std::invalid_argument( - "Grid is not one dimensional, you should probably " - "use getRow() or getCol()!"); - } - - return col; - } - /** * @brief Gets the number of rows of the grid. * * @return Number of rows. */ - int getRow() const { return this->row; } + int getRow() const { return this->concentrations.rows(); } /** * @brief Gets the number of columns of the grid. * * @return Number of columns. */ - int getCol() const { return this->col; } + int getCol() const { return this->concentrations.cols(); } /** * @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional. @@ -295,17 +273,12 @@ public: * @param domainLength A double value of the domain length. Must be positive. */ void setDomain(double domainLength) { - if (dim != 1) { - throw std::invalid_argument( - "Grid is not one dimensional, you should probaly " - "use the 2D domain setter!"); - } - if (domainLength <= 0) { - throw std::invalid_argument("Given domain length is not positive!"); - } + tug_assert(dim == 1, "Grid is not one dimensional, use 2D domain setter!"); + tug_assert(domainLength > 0, "Given domain length is not positive!"); this->domainCol = domainLength; - this->deltaCol = double(this->domainCol) / double(this->col); + this->deltaCol = + double(this->domainCol) / double(this->concentrations.cols()); } /** @@ -317,35 +290,18 @@ public: * be positive. */ void setDomain(double domainRow, double domainCol) { - if (dim != 2) { - throw std::invalid_argument( - "Grid is not two dimensional, you should probably " - "use the 1D domain setter!"); - } - if (domainRow <= 0 || domainCol <= 0) { - throw std::invalid_argument("Given domain size is not positive!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, use 1D domain setter!"); + tug_assert(domainCol > 0, + "Given domain size in x-direction is not positive!"); + tug_assert(domainRow > 0, + "Given domain size in y-direction is not positive!"); this->domainRow = domainRow; this->domainCol = domainCol; - this->deltaRow = double(this->domainRow) / double(this->row); - this->deltaCol = double(this->domainCol) / double(this->col); - } - - /** - * @brief Gets the delta value for 1D-Grid. Grid must be one dimensional. - * - * @return Delta value. - */ - T getDelta() const { - - if (dim != 1) { - throw std::invalid_argument( - "Grid is not one dimensional, you should probably " - "use the 2D delta getters"); - } - - return this->deltaCol; + this->deltaRow = + double(this->domainRow) / double(this->concentrations.rows()); + this->deltaCol = + double(this->domainCol) / double(this->concentrations.cols()); } /** @@ -361,23 +317,18 @@ public: * @return Delta value in y-direction. */ T getDeltaRow() const { - if (dim != 2) { - throw std::invalid_argument( - "Grid is not two dimensional, meaning there is no " - "delta in y-direction!"); - } + tug_assert(dim == 2, "Grid is not two dimensional, there is no delta in " + "y-direction!"); return this->deltaRow; } private: - int col; // number of grid columns - int row{1}; // number of grid rows int dim; // 1D or 2D T domainCol; // number of domain columns T domainRow{0}; // number of domain rows T deltaCol; // delta in x-direction (between columns) - T deltaRow{0}; // delta in y-direction (between rows) + T deltaRow; // delta in y-direction (between rows) RowMajMat concentrations; // Matrix holding grid concentrations RowMajMat alphaX; // Matrix holding alpha coefficients in x-direction @@ -389,4 +340,3 @@ private: using Grid64 = Grid; using Grid32 = Grid; } // namespace tug -#endif // GRID_H_ diff --git a/naaice/BTCS_2D_NAAICE.cpp b/naaice/BTCS_2D_NAAICE.cpp index f344893..15f199f 100644 --- a/naaice/BTCS_2D_NAAICE.cpp +++ b/naaice/BTCS_2D_NAAICE.cpp @@ -104,15 +104,14 @@ int main(int argc, char *argv[]) { // create a grid with a 5 x 10 field constexpr int row = 5; constexpr int col = 10; - Grid64 grid(row, col); // (optional) set the domain, e.g.: - grid.setDomain(0.005, 0.01); const auto init_values_vec = CSVToVector(INPUT_CONC_FILE); Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col); - grid.setConcentrations(concentrations); + Grid64 grid(concentrations); + grid.setDomain(0.005, 0.01); const double sum_init = concentrations.sum(); // // (optional) set alphas of the grid, e.g.: diff --git a/naaice/NAAICE_dble_vs_float.cpp b/naaice/NAAICE_dble_vs_float.cpp index c432f0f..2238375 100644 --- a/naaice/NAAICE_dble_vs_float.cpp +++ b/naaice/NAAICE_dble_vs_float.cpp @@ -114,15 +114,14 @@ template int doWork(int ngrid) { std::cout << name << " grid: " << ngrid << std::endl; - Grid grid(ngrid, ngrid); // Grid64 grid(ngrid, ngrid); // (optional) set the domain, e.g.: - grid.setDomain(0.1, 0.1); Eigen::MatrixX initConc64 = Eigen::MatrixX::Constant(ngrid, ngrid, 0); initConc64(50, 50) = 1E-6; - grid.setConcentrations(initConc64); + Grid grid(initConc64); + grid.setDomain(0.1, 0.1); const T sum_init64 = initConc64.sum(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 78d67e2..961b035 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -10,12 +10,13 @@ FetchContent_MakeAvailable(googletest) add_executable(testTug +setup.cpp testDiffusion.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp ) -target_link_libraries(testTug tug GTest::gtest_main) +target_link_libraries(testTug tug GTest::gtest) include(GoogleTest) gtest_discover_tests(testTug) diff --git a/test/setup.cpp b/test/setup.cpp index 0a3f254..b885e52 100644 --- a/test/setup.cpp +++ b/test/setup.cpp @@ -1,2 +1,7 @@ -#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN -#include +#include + +int main(int argc, char **argv) { + ::testing::InitGoogleTest(&argc, argv); + GTEST_FLAG_SET(death_test_style, "threadsafe"); + return RUN_ALL_TESTS(); +} \ No newline at end of file diff --git a/test/testBoundary.cpp b/test/testBoundary.cpp index 0b2601b..8b256f6 100644 --- a/test/testBoundary.cpp +++ b/test/testBoundary.cpp @@ -1,3 +1,5 @@ +#include "gtest/gtest.h" +#include #include #include #include @@ -15,7 +17,7 @@ BOUNDARY_TEST(Element) { EXPECT_NO_THROW(BoundaryElement()); EXPECT_EQ(boundaryElementClosed.getType(), BC_TYPE_CLOSED); EXPECT_DOUBLE_EQ(boundaryElementClosed.getValue(), -1); - EXPECT_ANY_THROW(boundaryElementClosed.setValue(0.2)); + EXPECT_THROW(boundaryElementClosed.setValue(0.2), std::invalid_argument); BoundaryElement boundaryElementConstant = BoundaryElement(0.1); EXPECT_NO_THROW(BoundaryElement(0.1)); @@ -26,8 +28,10 @@ BOUNDARY_TEST(Element) { } BOUNDARY_TEST(Class) { - Grid grid1D = Grid64(10); - Grid grid2D = Grid64(10, 12); + Eigen::VectorXd conc(10); + Grid grid1D = Grid64(conc); + Eigen::MatrixXd conc2D(10, 12); + Grid grid2D = Grid64(conc2D); Boundary boundary1D = Boundary(grid1D); Boundary boundary2D = Boundary(grid2D); vector> boundary1DVector(1, BoundaryElement(1.0)); @@ -48,20 +52,25 @@ BOUNDARY_TEST(Class) { EXPECT_EQ(boundary1D.getBoundarySide(BC_SIDE_RIGHT).size(), 1); EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CLOSED); - EXPECT_ANY_THROW(boundary1D.getBoundarySide(BC_SIDE_TOP)); - EXPECT_ANY_THROW(boundary1D.getBoundarySide(BC_SIDE_BOTTOM)); + EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_TOP), + ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*"); + EXPECT_DEATH(boundary1D.getBoundarySide(BC_SIDE_BOTTOM), + ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*"); EXPECT_NO_THROW(boundary1D.setBoundarySideClosed(BC_SIDE_LEFT)); - EXPECT_ANY_THROW(boundary1D.setBoundarySideClosed(BC_SIDE_TOP)); + EXPECT_DEATH(boundary1D.setBoundarySideClosed(BC_SIDE_TOP), + ".*BC_SIDE_LEFT .* BC_SIDE_RIGHT.*"); EXPECT_NO_THROW(boundary1D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); EXPECT_DOUBLE_EQ(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - EXPECT_ANY_THROW(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2)); + EXPECT_DEATH(boundary1D.getBoundaryElementValue(BC_SIDE_LEFT, 2), + ".*Index is selected either too large or too small.*"); EXPECT_EQ(boundary1D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); EXPECT_EQ(boundary1D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); EXPECT_NO_THROW(boundary1D.setInnerBoundary(0, inner_condition_value)); - EXPECT_ANY_THROW(boundary1D.setInnerBoundary(0, 0, inner_condition_value)); + EXPECT_DEATH(boundary1D.setInnerBoundary(0, 0, inner_condition_value), + ".*only available for 2D grids.*"); EXPECT_EQ(boundary1D.getInnerBoundary(0), innerBoundary); EXPECT_FALSE(boundary1D.getInnerBoundary(1).first); } @@ -80,13 +89,15 @@ BOUNDARY_TEST(Class) { EXPECT_NO_THROW(boundary2D.setBoundarySideClosed(BC_SIDE_TOP)); EXPECT_NO_THROW(boundary2D.setBoundarySideConstant(BC_SIDE_LEFT, 1.0)); EXPECT_DOUBLE_EQ(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 0), 1.0); - EXPECT_ANY_THROW(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12)); + EXPECT_DEATH(boundary2D.getBoundaryElementValue(BC_SIDE_LEFT, 12), + ".*too large or too small.*"); EXPECT_EQ(boundary2D.getBoundaryElementType(BC_SIDE_LEFT, 0), BC_TYPE_CONSTANT); EXPECT_EQ(boundary2D.getBoundaryElement(BC_SIDE_LEFT, 0).getType(), boundary1DVector[0].getType()); - EXPECT_ANY_THROW(boundary2D.setInnerBoundary(0, inner_condition_value)); + EXPECT_DEATH(boundary2D.setInnerBoundary(0, inner_condition_value), + ".* 1D .*"); EXPECT_NO_THROW(boundary2D.setInnerBoundary(0, 1, inner_condition_value)); EXPECT_EQ(boundary2D.getInnerBoundary(0, 1), innerBoundary); EXPECT_FALSE(boundary2D.getInnerBoundary(0, 2).first); diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 79d0506..b650d3e 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -1,4 +1,5 @@ #include "TestUtils.hpp" +#include "gtest/gtest.h" #include #include #include @@ -22,12 +23,11 @@ Grid64 setupSimulation(double timestep, int iterations) { int domain_col = 10; // Grid - Grid grid = Grid64(row, col); - grid.setDomain(domain_row, domain_col); - MatrixXd concentrations = MatrixXd::Constant(row, col, 0); concentrations(5, 5) = 1; - grid.setConcentrations(concentrations); + + Grid grid = Grid64(concentrations); + grid.setDomain(domain_row, domain_col); MatrixXd alpha = MatrixXd::Constant(row, col, 1); for (int i = 0; i < 5; i++) { @@ -96,7 +96,8 @@ DIFFUSION_TEST(EqualityBTCS) { DIFFUSION_TEST(InitializeEnvironment) { int rc = 12; - Grid64 grid(rc, rc); + Eigen::MatrixXd concentrations(rc, rc); + Grid64 grid(concentrations); Boundary boundary(grid); EXPECT_NO_THROW(Diffusion sim(grid, boundary)); @@ -104,7 +105,9 @@ DIFFUSION_TEST(InitializeEnvironment) { DIFFUSION_TEST(SimulationEnvironment) { int rc = 12; - Grid64 grid(rc, rc); + Eigen::MatrixXd concentrations(rc, rc); + Grid64 grid(concentrations); + grid.initAlpha(); Boundary boundary(grid); Diffusion sim(grid, boundary); @@ -116,7 +119,7 @@ DIFFUSION_TEST(SimulationEnvironment) { EXPECT_NO_THROW(sim.setTimestep(0.1)); EXPECT_DOUBLE_EQ(sim.getTimestep(), 0.1); - EXPECT_THROW(sim.setTimestep(-0.3), std::invalid_argument); + EXPECT_DEATH(sim.setTimestep(-0.3), ".* greater than zero.*"); } DIFFUSION_TEST(ClosedBoundaries) { @@ -124,13 +127,12 @@ DIFFUSION_TEST(ClosedBoundaries) { constexpr std::uint32_t nrows = 5; constexpr std::uint32_t ncols = 5; - tug::Grid64 grid(nrows, ncols); - auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0); auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - grid.setConcentrations(concentrations); + tug::Grid64 grid(concentrations); + grid.setAlpha(alphax, alphay); tug::Boundary bc(grid); @@ -153,13 +155,11 @@ DIFFUSION_TEST(ConstantInnerCell) { constexpr std::uint32_t nrows = 5; constexpr std::uint32_t ncols = 5; - tug::Grid64 grid(nrows, ncols); - auto concentrations = Eigen::MatrixXd::Constant(nrows, ncols, 1.0); auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5); - grid.setConcentrations(concentrations); + tug::Grid64 grid(concentrations); grid.setAlpha(alphax, alphay); tug::Boundary bc(grid); @@ -179,4 +179,4 @@ DIFFUSION_TEST(ConstantInnerCell) { EXPECT_FALSE((grid.getConcentrations().array() > 1.0).any()); EXPECT_FALSE((grid.getConcentrations().array() < 0.0).any()); -} \ No newline at end of file +} diff --git a/test/testGrid.cpp b/test/testGrid.cpp index 4539389..e8b7741 100644 --- a/test/testGrid.cpp +++ b/test/testGrid.cpp @@ -1,6 +1,6 @@ +#include "gtest/gtest.h" #include #include -#include #include @@ -10,41 +10,26 @@ using namespace tug; #define GRID_TEST(x) TEST(Grid, x) -GRID_TEST(InvalidConstructorParams) { - EXPECT_ANY_THROW(Grid64(2)); - EXPECT_ANY_THROW(Grid64(1, 4)); - EXPECT_ANY_THROW(Grid64(4, 1)); -} - GRID_TEST(Grid64OneDimensional) { int l = 12; - Grid64 grid(l); + Eigen::VectorXd conc(l); + Grid64 grid(conc); + grid.initAlpha(); { EXPECT_EQ(grid.getDim(), 1); - EXPECT_EQ(grid.getLength(), l); + EXPECT_EQ(grid.getCol(), l); EXPECT_EQ(grid.getCol(), l); EXPECT_EQ(grid.getRow(), 1); EXPECT_EQ(grid.getConcentrations().rows(), 1); EXPECT_EQ(grid.getConcentrations().cols(), l); - EXPECT_EQ(grid.getAlpha().rows(), 1); - EXPECT_EQ(grid.getAlpha().cols(), l); + EXPECT_EQ(grid.getAlphaX().rows(), 1); + EXPECT_EQ(grid.getAlphaX().cols(), l); EXPECT_EQ(grid.getDeltaCol(), 1); - EXPECT_ANY_THROW(grid.getAlphaX()); - EXPECT_ANY_THROW(grid.getAlphaY()); - EXPECT_ANY_THROW(grid.getDeltaRow()); - } - - { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(1, l, 3); - EXPECT_NO_THROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(2, l, 4); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); + EXPECT_DEATH(grid.getAlphaY(), ".* no alphaY!.*"); + EXPECT_DEATH(grid.getDeltaRow(), ".* not two dimensional, .*"); } { @@ -52,16 +37,16 @@ GRID_TEST(Grid64OneDimensional) { MatrixXd alpha = MatrixXd::Constant(1, l, 3); EXPECT_NO_THROW(grid.setAlpha(alpha)); - EXPECT_ANY_THROW(grid.setAlpha(alpha, alpha)); + EXPECT_DEATH(grid.setAlpha(alpha, alpha), ".* is not two dimensional, .*"); grid.setAlpha(alpha); - EXPECT_EQ(grid.getAlpha(), alpha); - EXPECT_ANY_THROW(grid.getAlphaX()); - EXPECT_ANY_THROW(grid.getAlphaY()); + EXPECT_EQ(grid.getAlphaX(), alpha); + EXPECT_NO_THROW(grid.getAlphaX()); + EXPECT_DEATH(grid.getAlphaY(), ".* no alphaY!.*"); // false alpha matrix MatrixXd wAlpha = MatrixXd::Constant(3, l, 2); - EXPECT_ANY_THROW(grid.setAlpha(wAlpha)); + EXPECT_DEATH(grid.setAlpha(wAlpha), ".* mismatch with Grid dimensions!.*"); } { @@ -70,33 +55,30 @@ GRID_TEST(Grid64OneDimensional) { EXPECT_NO_THROW(grid.setDomain(d)); // set 2D domain - EXPECT_ANY_THROW(grid.setDomain(d, d)); + EXPECT_DEATH(grid.setDomain(d, d), ".* not two dimensional, .*"); grid.setDomain(d); EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(d) / double(l)); - EXPECT_ANY_THROW(grid.getDeltaRow()); + EXPECT_DEATH(grid.getDeltaRow(), ".* not two dimensional, .*"); // set too small domain - d = 0; - EXPECT_ANY_THROW(grid.setDomain(d)); - d = -2; - EXPECT_ANY_THROW(grid.setDomain(d)); + EXPECT_DEATH(grid.setDomain(-2), "Given domain length .*"); } } GRID_TEST(Grid64Quadratic) { int rc = 12; - Grid64 grid(rc, rc); + Eigen::MatrixXd conc(rc, rc); + Grid64 grid(conc); + grid.initAlpha(); { EXPECT_EQ(grid.getDim(), 2); - EXPECT_ANY_THROW(grid.getLength()); EXPECT_EQ(grid.getCol(), rc); EXPECT_EQ(grid.getRow(), rc); EXPECT_EQ(grid.getConcentrations().rows(), rc); EXPECT_EQ(grid.getConcentrations().cols(), rc); - EXPECT_ANY_THROW(grid.getAlpha()); EXPECT_EQ(grid.getAlphaX().rows(), rc); EXPECT_EQ(grid.getAlphaX().cols(), rc); @@ -105,39 +87,25 @@ GRID_TEST(Grid64Quadratic) { EXPECT_EQ(grid.getDeltaRow(), 1); EXPECT_EQ(grid.getDeltaCol(), 1); } - - { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(rc, rc, 2); - EXPECT_NO_THROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(rc, rc + 3, 1); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc + 3, rc, 4); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(rc + 2, rc + 4, 6); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - } - { // correct alpha matrices MatrixXd alphax = MatrixXd::Constant(rc, rc, 2); MatrixXd alphay = MatrixXd::Constant(rc, rc, 4); EXPECT_NO_THROW(grid.setAlpha(alphax, alphay)); - EXPECT_ANY_THROW(grid.setAlpha(alphax)); + EXPECT_DEATH(grid.setAlpha(alphax), ".* 2D setter function!.*"); grid.setAlpha(alphax, alphay); EXPECT_EQ(grid.getAlphaX(), alphax); EXPECT_EQ(grid.getAlphaY(), alphay); - EXPECT_ANY_THROW(grid.getAlpha()); // false alpha matrices alphax = MatrixXd::Constant(rc + 3, rc + 1, 3); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); + EXPECT_DEATH(grid.setAlpha(alphax, alphay), + ".*has wrong number of rows!.*"); alphay = MatrixXd::Constant(rc + 2, rc + 1, 3); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); + EXPECT_DEATH(grid.setAlpha(alphax, alphay), + ".*has wrong number of rows!.*"); } { @@ -145,7 +113,7 @@ GRID_TEST(Grid64Quadratic) { int dc = 9; // set 1D domain - EXPECT_ANY_THROW(grid.setDomain(dr)); + EXPECT_DEATH(grid.setDomain(dr), ".* 2D domain setter!.*"); // set 2D domain EXPECT_NO_THROW(grid.setDomain(dr, dc)); @@ -156,29 +124,29 @@ GRID_TEST(Grid64Quadratic) { // set too small domain dr = 0; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); + EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*"); dr = 8; dc = 0; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); + EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*"); dr = -2; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); + EXPECT_DEATH(grid.setDomain(dr, dc), ".* not positive!.*"); } } GRID_TEST(Grid64NonQuadratic) { int r = 12; int c = 15; - Grid64 grid(r, c); + Eigen::MatrixXd conc(r, c); + Grid64 grid(conc); + grid.initAlpha(); { EXPECT_EQ(grid.getDim(), 2); - EXPECT_ANY_THROW(grid.getLength()); EXPECT_EQ(grid.getCol(), c); EXPECT_EQ(grid.getRow(), r); EXPECT_EQ(grid.getConcentrations().rows(), r); EXPECT_EQ(grid.getConcentrations().cols(), c); - EXPECT_ANY_THROW(grid.getAlpha()); EXPECT_EQ(grid.getAlphaX().rows(), r); EXPECT_EQ(grid.getAlphaX().cols(), c); @@ -188,75 +156,109 @@ GRID_TEST(Grid64NonQuadratic) { EXPECT_EQ(grid.getDeltaCol(), 1); } - { - // correct concentrations matrix - MatrixXd concentrations = MatrixXd::Constant(r, c, 2); - EXPECT_NO_THROW(grid.setConcentrations(concentrations)); - - // false concentrations matrix - MatrixXd wConcentrations = MatrixXd::Constant(r, c + 3, 6); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r + 3, c, 3); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - wConcentrations = MatrixXd::Constant(r + 2, c + 4, 2); - EXPECT_ANY_THROW(grid.setConcentrations(wConcentrations)); - } - { // correct alpha matrices MatrixXd alphax = MatrixXd::Constant(r, c, 2); MatrixXd alphay = MatrixXd::Constant(r, c, 4); EXPECT_NO_THROW(grid.setAlpha(alphax, alphay)); - EXPECT_ANY_THROW(grid.setAlpha(alphax)); - grid.setAlpha(alphax, alphay); EXPECT_EQ(grid.getAlphaX(), alphax); EXPECT_EQ(grid.getAlphaY(), alphay); - EXPECT_ANY_THROW(grid.getAlpha()); - - // false alpha matrices - alphax = MatrixXd::Constant(r + 3, c + 1, 3); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); - alphay = MatrixXd::Constant(r + 2, c + 1, 5); - EXPECT_ANY_THROW(grid.setAlpha(alphax, alphay)); } { int dr = 8; int dc = 9; - // set 1D domain - EXPECT_ANY_THROW(grid.setDomain(dr)); - // set 2D domain EXPECT_NO_THROW(grid.setDomain(dr, dc)); grid.setDomain(dr, dc); - EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(c)); - EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(r)); - - // set too small domain - dr = 0; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - dr = 8; - dc = -1; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); - dr = -2; - EXPECT_ANY_THROW(grid.setDomain(dr, dc)); + EXPECT_EQ(grid.getDeltaCol(), double(dc) / double(c)); + EXPECT_EQ(grid.getDeltaRow(), double(dr) / double(r)); } { + int r = 4; + int c = 5; std::vector concentrations(r * c); for (int i = 0; i < r * c; i++) { concentrations[i] = i; } + Grid64 grid(concentrations.data(), r, c); + grid.initAlpha(); - grid.setConcentrations(concentrations.data()); + { + EXPECT_EQ(grid.getDim(), 2); + EXPECT_EQ(grid.getCol(), c); + EXPECT_EQ(grid.getRow(), r); - EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0); - EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1); - EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c); + EXPECT_EQ(grid.getConcentrations().rows(), r); + EXPECT_EQ(grid.getConcentrations().cols(), c); + + EXPECT_EQ(grid.getAlphaX().rows(), r); + EXPECT_EQ(grid.getAlphaX().cols(), c); + EXPECT_EQ(grid.getAlphaY().rows(), r); + EXPECT_EQ(grid.getAlphaY().cols(), c); + EXPECT_EQ(grid.getDeltaRow(), 1); + EXPECT_EQ(grid.getDeltaCol(), 1); + } + + { + // correct alpha matrices + MatrixXd alphax = MatrixXd::Constant(r, c, 2); + MatrixXd alphay = MatrixXd::Constant(r, c, 4); + EXPECT_NO_THROW(grid.setAlpha(alphax, alphay)); + + EXPECT_DEATH(grid.setAlpha(alphax), ".* 2D setter function!.*"); + + grid.setAlpha(alphax, alphay); + EXPECT_EQ(grid.getAlphaX(), alphax); + EXPECT_EQ(grid.getAlphaY(), alphay); + + // false alpha matrices + alphax = MatrixXd::Constant(r + 3, c + 1, 3); + EXPECT_DEATH(grid.setAlpha(alphax, alphay), + ".*has wrong number of rows!.*"); + alphay = MatrixXd::Constant(r + 2, c + 1, 5); + EXPECT_DEATH(grid.setAlpha(alphax, alphay), + ".*has wrong number of rows!.*"); + + { + int dr = 8; + int dc = 9; + + // set 1D domain + EXPECT_DEATH(grid.setDomain(dr), ".* 2D domain setter!.*"); + + // set 2D domain + EXPECT_NO_THROW(grid.setDomain(dr, dc)); + + grid.setDomain(dr, dc); + EXPECT_DOUBLE_EQ(grid.getDeltaCol(), double(dc) / double(c)); + EXPECT_DOUBLE_EQ(grid.getDeltaRow(), double(dr) / double(r)); + } + + { + auto &concentrations = grid.getConcentrations(); + + for (int i = 0; i < r; i++) { + for (int j = 0; j < c; j++) { + concentrations(i, j) = i * c + j; + } + } + + EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0); + EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1); + EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c); + + EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 0), 0); + EXPECT_DOUBLE_EQ(grid.getConcentrations()(0, 1), 1); + EXPECT_DOUBLE_EQ(grid.getConcentrations()(1, 0), c); + EXPECT_DOUBLE_EQ(grid.getConcentrations()(2, 1), 2 * c + 1); + } + } } -} +} \ No newline at end of file