max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 01:18:22 +01:00
Code Issues Packages Projects Releases Wiki Activity

refactor: Update testDiffusion.cpp and Diffusion.hpp

Browse Source 
Refactor testDiffusion.cpp and Diffusion.hpp to improve code readability and maintainability. Remove unnecessary exception throwing and replace with assert statements for invalid arguments.
This commit is contained in:
Max Lübke 2024-06-13 09:16:51 +02:00
parent 13f6556f54
commit d3843fb2a3
16 changed files with 368 additions and 445 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitlab-ci.yml
View File

@ -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:

3
examples/BTCS_2D_proto_example.cpp
View File

@ -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

3
examples/FTCS_2D_proto_closed_mdl.cpp
View File

@ -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

3
examples/FTCS_2D_proto_example_mdl.cpp
View File

@ -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

126
include/tug/Boundary.hpp
View File

@ -8,6 +8,7 @@
#define BOUNDARY_H_
#include "Grid.hpp"
#include "tug/Core/TugUtils.hpp"
#include <cstddef>
#include <cstdint>
@ -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<T> objects.
*/
const std::vector<BoundaryElement<T>> &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<T> 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<bool, T> 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<bool, T> 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<std::pair<bool, T>> 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<std::pair<bool, T>>(this->cols,
@ -499,14 +463,8 @@ public:
* condition
*/
std::vector<std::pair<bool, T>> 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<std::pair<bool, T>>(this->rows,

6
include/tug/Core/Numeric/BTCS.hpp
View File

@ -354,15 +354,15 @@ template <class T>
static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &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<T> concentrations_t1(length);
Eigen::SparseMatrix<T> A;
Eigen::VectorX<T> 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);

35
include/tug/Core/Numeric/FTCS.hpp
View File

@ -9,7 +9,9 @@
#define FTCS_H_
#include "../TugUtils.hpp"
#include "tug/Core/Matrix.hpp"
#include <cstring>
#include <tug/Boundary.hpp>
#ifdef _OPENMP
@ -225,6 +227,7 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
int colMax = grid.getCol();
T deltaCol = grid.getDeltaCol();
RowMajMat<T> &concentrations_grid = grid.getConcentrations();
// matrix for concentrations at time t+1
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(1, colMax, 0);
@ -234,7 +237,7 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &bc, T timestep,
T deltaRow = grid.getDeltaRow();
T deltaCol = grid.getDeltaCol();
RowMajMat<T> &concentrations_grid = grid.getConcentrations();
// matrix for concentrations at time t+1
RowMajMat<T> concentrations_t1 = RowMajMat<T>::Constant(rowMax, colMax, 0);
@ -275,7 +281,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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<T> &grid, Boundary<T> &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));
// }
}

16
include/tug/Diffusion.hpp
View File

@ -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) {

324
include/tug/Grid.hpp
View File

@ -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 <Eigen/Core>
#include <Eigen/Sparse>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
#include <stdexcept>
#include <cstddef>
namespace tug {
@ -26,83 +26,97 @@ namespace tug {
template <class T> 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<T> 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<T> &concentrations) {
if (concentrations.rows() == 1) {
this->dim = 1;
this->domainCol = static_cast<T>(concentrations.cols());
this->deltaCol = static_cast<T>(this->domainCol) /
static_cast<T>(concentrations.cols()); // -> 1
this->concentrations = concentrations;
return;
}
this->dim = 1;
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
if (concentrations.cols() == 1) {
this->dim = 1;
this->domainCol = static_cast<T>(concentrations.rows());
this->deltaCol = static_cast<T>(this->domainCol) /
static_cast<T>(concentrations.rows()); // -> 1
this->concentrations = RowMajMat<T>::Constant(1, col, MAT_INIT_VAL);
this->alphaX = RowMajMat<T>::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<T>(this->domainRow) / static_cast<T>(this->row); // -> 1
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(this->col); // -> 1
this->concentrations = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
this->alphaX = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
this->alphaY = RowMajMat<T>::Constant(row, col, MAT_INIT_VAL);
}
/**
* @brief Sets the concentrations matrix for a 1D or 2D-Grid.
*
* @param concentrations An Eigen3 MatrixX<T> holding the concentrations.
* Matrix must have correct dimensions as defined in row and col. (Or length,
* in 1D case).
*/
void setConcentrations(const RowMajMat<T> &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<T>(concentrations.rows());
this->domainCol = static_cast<T>(concentrations.cols());
this->deltaRow = static_cast<T>(this->domainRow) /
static_cast<T>(concentrations.rows()); // -> 1
this->deltaCol = static_cast<T>(this->domainCol) /
static_cast<T>(concentrations.cols()); // -> 1
this->concentrations = concentrations;
// this->alphaX = RowMajMat<T>::Constant(concentrations.rows(),
// concentrations.cols(),
// MAT_INIT_VAL);
// this->alphaY = RowMajMat<T>::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<T> map(concentrations, this->row, this->col);
this->concentrations = map;
Grid(T *concentrations, std::size_t length) : dim(1) {
this->domainCol = static_cast<T>(length); // -> 1
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(length); // -> 1
this->concentrations = RowMajMatMap<T>(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<T>(row); // -> 1
this->domainCol = static_cast<T>(col); // -> 1
this->deltaCol =
static_cast<T>(this->domainCol) / static_cast<T>(col); // -> 1
this->deltaRow =
static_cast<T>(this->domainRow) / static_cast<T>(row); // -> 1
this->concentrations = RowMajMatMap<T>(concentrations, row, col);
}
/**
@ -113,6 +127,17 @@ public:
*/
auto &getConcentrations() { return this->concentrations; }
void initAlpha() {
this->alphaX = RowMajMat<T>::Constant(
this->concentrations.rows(), this->concentrations.cols(), MAT_INIT_VAL);
if (dim > 1) {
this->alphaY =
RowMajMat<T>::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<T> &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<T> 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<T> &alphaX, const RowMajMat<T> &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<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 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<T> concentrations; // Matrix holding grid concentrations
RowMajMat<T> alphaX; // Matrix holding alpha coefficients in x-direction
@ -389,4 +340,3 @@ private:
using Grid64 = Grid<double>;
using Grid32 = Grid<float>;
} // namespace tug
#endif // GRID_H_

5
naaice/BTCS_2D_NAAICE.cpp
View File

@ -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<double>(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.:

5
naaice/NAAICE_dble_vs_float.cpp
View File

@ -114,15 +114,14 @@ template <class T, tug::APPROACH app> int doWork(int ngrid) {
std::cout << name << " grid: " << ngrid << std::endl;
Grid<T> grid(ngrid, ngrid);
// Grid64 grid(ngrid, ngrid);
// (optional) set the domain, e.g.:
grid.setDomain(0.1, 0.1);
Eigen::MatrixX<T> initConc64 = Eigen::MatrixX<T>::Constant(ngrid, ngrid, 0);
initConc64(50, 50) = 1E-6;
grid.setConcentrations(initConc64);
Grid<T> grid(initConc64);
grid.setDomain(0.1, 0.1);
const T sum_init64 = initConc64.sum();

3
test/CMakeLists.txt
View File

@ -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)

9
test/setup.cpp
View File

@ -1,2 +1,7 @@
#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
#include <doctest/doctest.h>
#include <gtest/gtest.h>
int main(int argc, char **argv) {
::testing::InitGoogleTest(&argc, argv);
GTEST_FLAG_SET(death_test_style, "threadsafe");
return RUN_ALL_TESTS();
}

31
test/testBoundary.cpp
View File

@ -1,3 +1,5 @@
#include "gtest/gtest.h"
#include <stdexcept>
#include <tug/Boundary.hpp>
#include <utility>
#include <vector>
@ -15,7 +17,7 @@ BOUNDARY_TEST(Element) {
EXPECT_NO_THROW(BoundaryElement<double>());
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<BoundaryElement<double>> 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);

28
test/testDiffusion.cpp
View File

@ -1,4 +1,5 @@
#include "TestUtils.hpp"
#include "gtest/gtest.h"
#include <gtest/gtest.h>
#include <stdexcept>
#include <tug/Diffusion.hpp>
@ -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<double, tug::FTCS_APPROACH> 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());
}
}

214
test/testGrid.cpp
View File

@ -1,6 +1,6 @@
#include "gtest/gtest.h"
#include <Eigen/Core>
#include <tug/Grid.hpp>
#include <vector>
#include <gtest/gtest.h>
@ -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<double> 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);
}
}
}
}
}