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

[wip] save state

Browse Source
This commit is contained in:
Max Lübke 2024-06-17 08:47:11 +02:00
parent 8c18e82b15
commit 0e840c51b3
17 changed files with 387 additions and 534 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/BTCS_2D_proto_example.cpp
View File

@ -23,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;
Grid64 grid(concentrations);
UniformGrid64 grid(concentrations);
// (optional) set alphas of the grid, e.g.:
// MatrixXd alphax = MatrixXd::Constant(20,20,1); // row,col,value

2
examples/FTCS_2D_proto_closed_mdl.cpp
View File

@ -43,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;
Grid64 grid(concentrations);
UniformGrid64 grid(concentrations);
// (optional) set alphas of the grid, e.g.:
MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value

2
examples/FTCS_2D_proto_example_mdl.cpp
View File

@ -34,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;
Grid64 grid(concentrations);
UniformGrid64 grid(concentrations);
// (optional) set alphas of the grid, e.g.:
MatrixXd alphax = MatrixXd::Constant(row, col, 1E-4); // row,col,value

8
include/tug/Boundary.hpp
View File

@ -7,7 +7,7 @@
#ifndef BOUNDARY_H_
#define BOUNDARY_H_
#include "Grid.hpp"
#include "UniformGrid.hpp"
#include "tug/Core/TugUtils.hpp"
#include <cstddef>
@ -115,7 +115,7 @@ public:
*
* @param length Length of the grid
*/
Boundary(std::uint32_t length) : Boundary(Grid<T>(length)){};
Boundary(std::uint32_t length) : Boundary(UnfiormGrid<T>(length)){};
/**
* @brief Creates a boundary object for a 2D grid
@ -124,7 +124,7 @@ public:
* @param cols Number of columns of the grid
*/
Boundary(std::uint32_t rows, std::uint32_t cols)
: Boundary(Grid<T>(rows, cols)){};
: Boundary(UnfiormGrid<T>(rows, cols)){};
/**
* @brief Creates a boundary object based on the passed grid object and
@ -133,7 +133,7 @@ public:
* @param grid Grid object on the basis of which the simulation takes place
* and from which the dimensions (in 2D case) are taken.
*/
Boundary(const Grid<T> &grid)
Boundary(const UnfiormGrid<T> &grid)
: dim(grid.getDim()), cols(grid.getCol()), rows(grid.getRow()) {
if (this->dim == 1) {
this->boundaries = std::vector<std::vector<BoundaryElement<T>>>(

45
include/tug/Core/Matrix.hpp
View File

@ -1,6 +1,8 @@
#pragma once
#include <Eigen/Core>
#include <Eigen/src/Core/Matrix.h>
#include <Eigen/src/Core/util/Constants.h>
namespace tug {
/**
@ -13,9 +15,48 @@ namespace tug {
*
* @tparam T The type of the matrix elements.
*/
// template <typename T>
// using RowMajMat =
// Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
template <typename T>
using RowMajMat =
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
class RowMajMat
: public Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> {
protected:
std::uint8_t dim;
public:
RowMajMat(Eigen::Index rows, Eigen::Index cols, T initial_value) : dim(2) {
*this = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
Eigen::RowMajor>::Constant(rows, cols, initial_value);
};
RowMajMat(Eigen::Index n_cells, T initial_value) : dim(1) {
*this = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
Eigen::RowMajor>::Constant(1, n_cells, initial_value);
};
/**
* @brief Gets the number of rows of the grid.
*
* @return Number of rows.
*/
Eigen::Index getRow() const { return this->rows(); }
/**
* @brief Gets the number of columns of the grid.
*
* @return Number of columns.
*/
Eigen::Index getCol() const { return this->cols(); }
/**
* @brief Gets the dimensions of the grid.
*
* @return Dimensions, either 1 or 2.
*/
int getDim() const { return this->dim; }
};
template <typename T> using RowMajMatMap = Eigen::Map<RowMajMat<T>>;
} // namespace tug

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

@ -351,7 +351,7 @@ static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
// BTCS solution for 1D grid
template <class T>
static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
static void BTCS_1D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b)) {
int length = grid.getCol();
@ -396,7 +396,7 @@ static void BTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep,
// BTCS solution for 2D grid
template <class T>
static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
static void BTCS_2D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b),
int numThreads) {
@ -449,7 +449,8 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
// entry point for EigenLU solver; differentiate between 1D and 2D grid
template <class T>
void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
void BTCS_LU(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, EigenLUAlgorithm);
} else if (grid.getDim() == 2) {
@ -462,7 +463,8 @@ void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
// entry point for Thomas algorithm solver; differentiate 1D and 2D grid
template <class T>
void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep, int numThreads) {
void BTCS_Thomas(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
int numThreads) {
if (grid.getDim() == 1) {
BTCS_1D(grid, bc, timestep, ThomasAlgorithm);
} else if (grid.getDim() == 2) {

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

@ -24,7 +24,7 @@ namespace tug {
// calculates horizontal change on one cell independent of boundary type
template <class T>
static inline T calcHorizontalChange(Grid<T> &grid, int &row, int &col) {
static inline T calcHorizontalChange(UnfiormGrid<T> &grid, int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
@ -41,7 +41,7 @@ static inline T calcHorizontalChange(Grid<T> &grid, int &row, int &col) {
// calculates vertical change on one cell independent of boundary type
template <class T>
static inline T calcVerticalChange(Grid<T> &grid, int &row, int &col) {
static inline T calcVerticalChange(UnfiormGrid<T> &grid, int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
@ -58,7 +58,7 @@ static inline T calcVerticalChange(Grid<T> &grid, int &row, int &col) {
// calculates horizontal change on one cell with a constant left boundary
template <class T>
static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid,
static inline T calcHorizontalChangeLeftBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc,
int &row, int &col) {
@ -75,8 +75,8 @@ static inline T calcHorizontalChangeLeftBoundaryConstant(Grid<T> &grid,
// calculates horizontal change on one cell with a closed left boundary
template <class T>
static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
static inline T calcHorizontalChangeLeftBoundaryClosed(UnfiormGrid<T> &grid,
int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaX()(row, col + 1),
grid.getAlphaX()(row, col)) *
@ -86,8 +86,9 @@ static inline T calcHorizontalChangeLeftBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the left edge of grid
template <class T>
static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &bc,
int &row, int &col) {
static inline T calcHorizontalChangeLeftBoundary(UnfiormGrid<T> &grid,
Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CONSTANT) {
return calcHorizontalChangeLeftBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_LEFT, row) == BC_TYPE_CLOSED) {
@ -99,7 +100,7 @@ static inline T calcHorizontalChangeLeftBoundary(Grid<T> &grid, Boundary<T> &bc,
// calculates horizontal change on one cell with a constant right boundary
template <class T>
static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid,
static inline T calcHorizontalChangeRightBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc,
int &row, int &col) {
@ -116,8 +117,8 @@ static inline T calcHorizontalChangeRightBoundaryConstant(Grid<T> &grid,
// calculates horizontal change on one cell with a closed right boundary
template <class T>
static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
static inline T calcHorizontalChangeRightBoundaryClosed(UnfiormGrid<T> &grid,
int &row, int &col) {
return -(calcAlphaIntercell(grid.getAlphaX()(row, col - 1),
grid.getAlphaX()(row, col)) *
@ -127,7 +128,7 @@ static inline T calcHorizontalChangeRightBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the right edge of grid
template <class T>
static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid,
static inline T calcHorizontalChangeRightBoundary(UnfiormGrid<T> &grid,
Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_RIGHT, row) == BC_TYPE_CONSTANT) {
@ -141,7 +142,7 @@ static inline T calcHorizontalChangeRightBoundary(Grid<T> &grid,
// calculates vertical change on one cell with a constant top boundary
template <class T>
static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid,
static inline T calcVerticalChangeTopBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc, int &row,
int &col) {
@ -158,8 +159,8 @@ static inline T calcVerticalChangeTopBoundaryConstant(Grid<T> &grid,
// calculates vertical change on one cell with a closed top boundary
template <class T>
static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
static inline T calcVerticalChangeTopBoundaryClosed(UnfiormGrid<T> &grid,
int &row, int &col) {
return calcAlphaIntercell(grid.getAlphaY()(row + 1, col),
grid.getAlphaY()(row, col)) *
@ -169,8 +170,9 @@ static inline T calcVerticalChangeTopBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the top edge of grid
template <class T>
static inline T calcVerticalChangeTopBoundary(Grid<T> &grid, Boundary<T> &bc,
int &row, int &col) {
static inline T calcVerticalChangeTopBoundary(UnfiormGrid<T> &grid,
Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeTopBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_TOP, col) == BC_TYPE_CLOSED) {
@ -182,7 +184,7 @@ static inline T calcVerticalChangeTopBoundary(Grid<T> &grid, Boundary<T> &bc,
// calculates vertical change on one cell with a constant bottom boundary
template <class T>
static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid,
static inline T calcVerticalChangeBottomBoundaryConstant(UnfiormGrid<T> &grid,
Boundary<T> &bc,
int &row, int &col) {
@ -199,8 +201,8 @@ static inline T calcVerticalChangeBottomBoundaryConstant(Grid<T> &grid,
// calculates vertical change on one cell with a closed bottom boundary
template <class T>
static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &grid, int &row,
int &col) {
static inline T calcVerticalChangeBottomBoundaryClosed(UnfiormGrid<T> &grid,
int &row, int &col) {
return -(calcAlphaIntercell(grid.getAlphaY()(row, col),
grid.getAlphaY()(row - 1, col)) *
@ -210,8 +212,9 @@ static inline T calcVerticalChangeBottomBoundaryClosed(Grid<T> &grid, int &row,
// checks boundary condition type for a cell on the bottom edge of grid
template <class T>
static inline T calcVerticalChangeBottomBoundary(Grid<T> &grid, Boundary<T> &bc,
int &row, int &col) {
static inline T calcVerticalChangeBottomBoundary(UnfiormGrid<T> &grid,
Boundary<T> &bc, int &row,
int &col) {
if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CONSTANT) {
return calcVerticalChangeBottomBoundaryConstant(grid, bc, row, col);
} else if (bc.getBoundaryElementType(BC_SIDE_BOTTOM, col) == BC_TYPE_CLOSED) {
@ -223,7 +226,7 @@ static inline T calcVerticalChangeBottomBoundary(Grid<T> &grid, Boundary<T> &bc,
// FTCS solution for 1D grid
template <class T>
static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
static void FTCS_1D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep) {
int colMax = grid.getCol();
T deltaCol = grid.getDeltaCol();
@ -263,7 +266,7 @@ static void FTCS_1D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
// FTCS solution for 2D grid
template <class T>
static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
static void FTCS_2D(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep,
int numThreads) {
int rowMax = grid.getRow();
int colMax = grid.getCol();
@ -392,7 +395,7 @@ static void FTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
// entry point; differentiate between 1D and 2D grid
template <class T>
void FTCS(Grid<T> &grid, Boundary<T> &bc, T timestep, int &numThreads) {
void FTCS(UnfiormGrid<T> &grid, Boundary<T> &bc, T timestep, int &numThreads) {
if (grid.getDim() == 1) {
FTCS_1D(grid, bc, timestep);
} else if (grid.getDim() == 2) {

8
include/tug/Diffusion.hpp
View File

@ -9,7 +9,7 @@
#pragma once
#include "Boundary.hpp"
#include "Grid.hpp"
#include "UniformGrid.hpp"
#include <algorithm>
#include <filesystem>
#include <fstream>
@ -69,7 +69,7 @@ private:
int innerIterations{1};
int numThreads{omp_get_num_procs()};
Grid<T> &grid;
UnfiormGrid<T> &grid;
Boundary<T> &bc;
const std::vector<std::string> approach_names = {"FTCS", "BTCS", "CRNI"};
@ -87,7 +87,7 @@ public:
* @param bc Valid boundary condition object
* @param approach Approach to solving the problem. Either FTCS or BTCS.
*/
Diffusion(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
Diffusion(UnfiormGrid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
/**
* @brief Setting the time step for each iteration step. Time step must be
@ -96,7 +96,7 @@ public:
* @param timestep Valid timestep greater than zero.
*/
void setTimestep(T timestep) {
tug_assert(timestep > 0, "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) {

342
include/tug/Grid.hpp
View File

@ -1,342 +0,0 @@
#pragma once
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients.
*
*/
#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 <cstddef>
namespace tug {
/**
* @brief Holds a matrix with concenctration and respective matrix/matrices of
* alpha coefficients.
*
* @tparam T Type to be used for matrices, e.g. double or float
*/
template <class T> class Grid {
public:
/**
* @brief Construct a new Grid object.
*
* 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(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;
}
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 = concentrations.transpose();
return;
}
this->dim = 2;
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 Construct a new Grid object.
*
* 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.
*/
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);
}
/**
* @brief Gets the concentrations matrix for a Grid.
*
* @return An Eigen3 matrix holding the concentrations and having
* the same dimensions as the grid.
*/
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.
*
* @param alpha An Eigen3 MatrixX<T> with 1 row holding the alpha
* coefficients. Matrix columns must have same size as length of grid.
*/
void setAlpha(const RowMajMat<T> &alpha) {
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;
}
/**
* @brief Set the alpha coefficients of a 1D-Grid. Grid must be one
* dimensional.
*
* @param alpha A pointer to an array holding the alpha coefficients. Array
* must have correct dimensions as defined in length. There is no check for
* correct dimensions, so be careful!
*/
void setAlpha(T *alpha) {
tug_assert(dim == 1,
"Grid is not one dimensional, use 2D setter function!");
RowMajMatMap<T> map(alpha, 1, this->col);
this->alphaX = map;
}
/**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
* dimensional.
*
* @param alphaX An Eigen3 MatrixX<T> holding the alpha coefficients in
* x-direction. Matrix must be of same size as the grid.
* @param alphaY An Eigen3 MatrixX<T> holding the alpha coefficients in
* y-direction. Matrix must be of same size as the grid.
*/
void setAlpha(const RowMajMat<T> &alphaX, const RowMajMat<T> &alphaY) {
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;
}
/**
* @brief Set the alpha coefficients of a 2D-Grid. Grid must be two
* dimensional.
*
* @param alphaX A pointer to an array holding the alpha coefficients in
* x-direction. Array must have correct dimensions as defined in row and col.
* There is no check for correct dimensions, so be careful!
* @param alphaY A pointer to an array holding the alpha coefficients in
* y-direction. Array must have correct dimensions as defined in row and col.
* There is no check for correct dimensions, so be careful!
*/
void setAlpha(T *alphaX, T *alphaY) {
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 in x-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return A matrix holding the alpha coefficients in x-direction.
*/
const auto &getAlphaX() const {
tug_assert(this->alphaX.size() > 0, "AlphaX is empty!");
return this->alphaX;
}
/**
* @brief Gets the matrix of alpha coefficients in y-direction of a 2D-Grid.
* Grid must be two dimensional.
*
* @return A matrix holding the alpha coefficients in y-direction.
*/
const auto &getAlphaY() const {
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;
}
/**
* @brief Gets the dimensions of the grid.
*
* @return Dimensions, either 1 or 2.
*/
int getDim() const { return this->dim; }
/**
* @brief Gets the number of rows of the grid.
*
* @return Number of rows.
*/
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->concentrations.cols(); }
/**
* @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional.
*
* @param domainLength A double value of the domain length. Must be positive.
*/
void setDomain(double domainLength) {
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->concentrations.cols());
}
/**
* @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional.
*
* @param domainRow A double value of the domain size in y-direction. Must
* be positive.
* @param domainCol A double value of the domain size in x-direction. Must
* be positive.
*/
void setDomain(double domainRow, double domainCol) {
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->concentrations.rows());
this->deltaCol =
double(this->domainCol) / double(this->concentrations.cols());
}
/**
* @brief Gets the delta value in x-direction.
*
* @return Delta value in x-direction.
*/
T getDeltaCol() const { return this->deltaCol; }
/**
* @brief Gets the delta value in y-direction. Must be two dimensional grid.
*
* @return Delta value in y-direction.
*/
T getDeltaRow() const {
tug_assert(dim == 2, "Grid is not two dimensional, there is no delta in "
"y-direction!");
return this->deltaRow;
}
private:
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; // delta in y-direction (between rows)
RowMajMat<T> concentrations; // Matrix holding grid concentrations
RowMajMat<T> alphaX; // Matrix holding alpha coefficients in x-direction
RowMajMat<T> alphaY; // Matrix holding alpha coefficients in y-direction
static constexpr T MAT_INIT_VAL = 0;
};
using Grid64 = Grid<double>;
using Grid32 = Grid<float>;
} // namespace tug

30
include/tug/UniformAlpha.hpp Normal file
View File

@ -0,0 +1,30 @@
#pragma once
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients.
*
*/
#include "tug/UniformGrid.hpp"
#include <Eigen/Core>
namespace tug {
/**
* @brief Holds a matrix with concenctration and respective matrix/matrices of
* alpha coefficients.
*
* @tparam T Type to be used for matrices, e.g. double or float
*/
template <class T> class UniformAlpha : public RowMajMat<T> {
public:
UniformAlpha(const UniformGrid<T> &grid, T init_value = 0)
: RowMajMat<T>::Constant(grid.rows(), grid.cols(), init_value) {}
UniformAlpha(const UniformGrid<T> &grid, T *alpha)
: tug::RowMajMatMap<T>(alpha, grid.rows(), grid.cols()) {}
};
} // namespace tug

117
include/tug/UniformGrid.hpp Normal file
View File

@ -0,0 +1,117 @@
#pragma once
/**
* @file Grid.hpp
* @brief API of Grid class, that holds a matrix with concenctrations and a
* respective matrix/matrices of alpha coefficients.
*
*/
#include "Core/Matrix.hpp"
#include "tug/Core/TugUtils.hpp"
#include <Eigen/Core>
namespace tug {
/**
* @brief Holds a matrix with concenctration and respective matrix/matrices of
* alpha coefficients.
*
* @tparam T Type to be used for matrices, e.g. double or float
*/
template <class T> class UniformGrid : public RowMajMat<T> {
public:
UniformGrid(Eigen::Index n_cells, T spatial_size, T initial_value)
: RowMajMat<T>(n_cells, initial_value), domainCol(spatial_size) {
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(n_cells);
}
UniformGrid(Eigen::Index row, Eigen::Index col, T spatial_row_size,
T spatial_col_size, T initial_value)
: RowMajMat<T>(row, col, initial_value), domainRow(spatial_row_size),
domainCol(spatial_col_size) {
this->deltaRow = static_cast<T>(this->domainRow) / static_cast<T>(row);
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(col);
}
UniformGrid(T *concentrations, Eigen::Index n_cells, T spatial_size)
: RowMajMatMap<T>(concentrations, 1, n_cells), domainCol(spatial_size) {
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(n_cells);
}
UniformGrid(T *concentrations, Eigen::Index row, Eigen::Index col,
T spatial_row_size, T spatial_col_size)
: RowMajMatMap<T>(concentrations, row, col), domainRow(spatial_row_size),
domainCol(spatial_col_size) {
this->deltaRow = static_cast<T>(this->domainRow) / static_cast<T>(row);
this->deltaCol = static_cast<T>(this->domainCol) / static_cast<T>(col);
}
/**
* @brief Sets the domain length of a 1D-Grid. Grid must be one dimensional.
*
* @param domainLength A double value of the domain length. Must be positive.
*/
void setDomain(double domainLength) {
tug_assert(this->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->cols());
}
/**
* @brief Sets the domain size of a 2D-Grid. Grid must be two dimensional.
*
* @param domainRow A double value of the domain size in y-direction. Must
* be positive.
* @param domainCol A double value of the domain size in x-direction. Must
* be positive.
*/
void setDomain(double domainRow, double domainCol) {
tug_assert(this->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->concentrations.rows());
this->deltaCol =
double(this->domainCol) / double(this->concentrations.cols());
}
/**
* @brief Gets the delta value in x-direction.
*
* @return Delta value in x-direction.
*/
T getDeltaCol() const { return this->deltaCol; }
/**
* @brief Gets the delta value in y-direction. Must be two dimensional grid.
*
* @return Delta value in y-direction.
*/
T getDeltaRow() const {
tug_assert(this->dim == 2,
"Grid is not two dimensional, there is no delta in "
"y-direction!");
return this->deltaRow;
}
private:
T domainCol; // number of domain columns
T domainRow{0}; // number of domain rows
T deltaCol; // delta in x-direction (between columns)
T deltaRow; // delta in y-direction (between rows)
};
using UniformGrid64 = UniformGrid<double>;
using UniformGrid32 = UniformGrid<float>;
} // namespace tug

2
naaice/BTCS_2D_NAAICE.cpp
View File

@ -109,7 +109,7 @@ int main(int argc, char *argv[]) {
const auto init_values_vec = CSVToVector<double>(INPUT_CONC_FILE);
Eigen::MatrixXd concentrations = rmVecTocmMatrix(init_values_vec, row, col);
Grid64 grid(concentrations);
UniformGrid64 grid(concentrations);
grid.setDomain(0.005, 0.01);
const double sum_init = concentrations.sum();

5
naaice/NAAICE_dble_vs_float.cpp
View File

@ -1,3 +1,4 @@
#include "tug/UniformGrid.hpp"
#include <Eigen/Eigen>
#include <chrono>
#include <cstdint>
@ -120,8 +121,8 @@ template <class T, tug::APPROACH app> int doWork(int ngrid) {
Eigen::MatrixX<T> initConc64 = Eigen::MatrixX<T>::Constant(ngrid, ngrid, 0);
initConc64(50, 50) = 1E-6;
Grid<T> grid(initConc64);
grid.setDomain(0.1, 0.1);
UniformGrid<T> grid(initConc64.data(), ngrid, ngrid);
const T sum_init64 = initConc64.sum();

3
test/CMakeLists.txt
View File

@ -1,6 +1,7 @@
find_package(doctest REQUIRED)
add_executable(testTug setup.cpp testDiffusion.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp)
# add_executable(testTug setup.cpp testDiffusion.cpp testGrid.cpp testFTCS.cpp testBoundary.cpp)
add_executable(testTug setup.cpp testGrid.cpp)
target_link_libraries(testTug doctest::doctest tug)
# get relative path of the CSV file

4
test/testBoundary.cpp
View File

@ -34,9 +34,9 @@ TEST_CASE("BoundaryElement") {
TEST_CASE("Boundary Class") {
Eigen::VectorXd conc(10);
Grid grid1D = Grid64(conc);
UnfiormGrid grid1D = UniformGrid64(conc);
Eigen::MatrixXd conc2D(10, 12);
Grid grid2D = Grid64(conc2D);
UnfiormGrid grid2D = UniformGrid64(conc2D);
Boundary boundary1D = Boundary(grid1D);
Boundary boundary2D = Boundary(grid2D);
vector<BoundaryElement<double>> boundary1DVector(1, BoundaryElement(1.0));

16
test/testDiffusion.cpp
View File

@ -12,7 +12,7 @@ using namespace Eigen;
using namespace std;
using namespace tug;
Grid64 setupSimulation(double timestep, int iterations) {
UniformGrid64 setupSimulation(double timestep, int iterations) {
int row = 11;
int col = 11;
int domain_row = 10;
@ -22,7 +22,7 @@ Grid64 setupSimulation(double timestep, int iterations) {
MatrixXd concentrations = MatrixXd::Constant(row, col, 0);
concentrations(5, 5) = 1;
Grid grid = Grid64(concentrations);
UnfiormGrid grid = UniformGrid64(concentrations);
grid.setDomain(domain_row, domain_col);
MatrixXd alpha = MatrixXd::Constant(row, col, 1);
@ -55,7 +55,7 @@ TEST_CASE("equality to reference matrix with FTCS") {
MatrixXd reference = CSV2Eigen(test_path);
cout << "FTCS Test: " << endl;
Grid grid = setupSimulation(timestep, iterations); // Boundary
UnfiormGrid grid = setupSimulation(timestep, iterations); // Boundary
Boundary bc = Boundary(grid);
// Simulation
@ -76,7 +76,7 @@ TEST_CASE("equality to reference matrix with BTCS") {
MatrixXd reference = CSV2Eigen(test_path);
cout << "BTCS Test: " << endl;
Grid grid = setupSimulation(timestep, iterations); // Boundary
UnfiormGrid grid = setupSimulation(timestep, iterations); // Boundary
Boundary bc = Boundary(grid);
// Simulation
@ -93,7 +93,7 @@ TEST_CASE("equality to reference matrix with BTCS") {
TEST_CASE("Initialize environment") {
int rc = 12;
Eigen::MatrixXd concentrations(rc, rc);
Grid64 grid(concentrations);
UniformGrid64 grid(concentrations);
Boundary boundary(grid);
CHECK_NOTHROW(Diffusion sim(grid, boundary));
@ -102,7 +102,7 @@ TEST_CASE("Initialize environment") {
TEST_CASE("Simulation environment") {
int rc = 12;
Eigen::MatrixXd concentrations(rc, rc);
Grid64 grid(concentrations);
UniformGrid64 grid(concentrations);
grid.initAlpha();
Boundary boundary(grid);
Diffusion<double, tug::FTCS_APPROACH> sim(grid, boundary);
@ -129,7 +129,7 @@ TEST_CASE("Closed Boundaries - no change expected") {
auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
tug::Grid64 grid(concentrations);
tug::UniformGrid64 grid(concentrations);
grid.setAlpha(alphax, alphay);
@ -158,7 +158,7 @@ TEST_CASE("Constant inner cell - 'absorbing' concentration") {
auto alphax = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
auto alphay = Eigen::MatrixXd::Constant(nrows, ncols, 1E-5);
tug::Grid64 grid(concentrations);
tug::UniformGrid64 grid(concentrations);
grid.setAlpha(alphax, alphay);
tug::Boundary bc(grid);

274
test/testGrid.cpp
View File

@ -1,7 +1,8 @@
#include "tug/UniformGrid.hpp"
#include <Eigen/Core>
#include <Eigen/src/Core/Matrix.h>
#include <doctest/doctest.h>
#include <tug/Grid.hpp>
#include <tug/UniformGrid.hpp>
#include <vector>
using namespace Eigen;
@ -11,30 +12,29 @@ using namespace tug;
TEST_CASE("1D Grid64") {
int l = 12;
Eigen::VectorXd conc(l);
Grid64 grid(conc);
grid.initAlpha();
UniformGrid64 grid(l, l, 1);
SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 1);
CHECK_EQ(grid.getCol(), l);
CHECK_EQ(grid.getCol(), l);
CHECK_EQ(grid.getRow(), 1);
CHECK_EQ(grid.getConcentrations().rows(), 1);
CHECK_EQ(grid.getConcentrations().cols(), l);
CHECK_EQ(grid.getAlphaX().rows(), 1);
CHECK_EQ(grid.getAlphaX().cols(), l);
CHECK_EQ(grid.getDeltaCol(), 1);
// CHECK_EQ(grid.getConcentrations().rows(), 1);
// CHECK_EQ(grid.getConcentrations().cols(), l);
// CHECK_EQ(grid.getAlphaX().rows(), 1);
// CHECK_EQ(grid.getAlphaX().cols(), l);
}
SUBCASE("setting alpha") {
// correct alpha matrix
MatrixXd alpha = MatrixXd::Constant(1, l, 3);
CHECK_NOTHROW(grid.setAlpha(alpha));
// SUBCASE("setting alpha") {
// // correct alpha matrix
// MatrixXd alpha = MatrixXd::Constant(1, l, 3);
// CHECK_NOTHROW(grid.setAlpha(alpha));
grid.setAlpha(alpha);
CHECK_EQ(grid.getAlphaX(), alpha);
}
// grid.setAlpha(alpha);
// CHECK_EQ(grid.getAlphaX(), alpha);
// }
SUBCASE("setting domain") {
int d = 8;
@ -46,153 +46,153 @@ TEST_CASE("1D Grid64") {
}
}
TEST_CASE("2D Grid64 quadratic") {
int rc = 12;
Eigen::MatrixXd conc(rc, rc);
Grid64 grid(conc);
grid.initAlpha();
// TEST_CASE("2D Grid64 quadratic") {
// int rc = 12;
// Eigen::MatrixXd conc(rc, rc);
// Grid64 grid(conc);
// grid.initAlpha();
SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 2);
CHECK_EQ(grid.getCol(), rc);
CHECK_EQ(grid.getRow(), rc);
// SUBCASE("correct construction") {
// CHECK_EQ(grid.getDim(), 2);
// CHECK_EQ(grid.getCol(), rc);
// CHECK_EQ(grid.getRow(), rc);
CHECK_EQ(grid.getConcentrations().rows(), rc);
CHECK_EQ(grid.getConcentrations().cols(), rc);
// CHECK_EQ(grid.getConcentrations().rows(), rc);
// CHECK_EQ(grid.getConcentrations().cols(), rc);
CHECK_EQ(grid.getAlphaX().rows(), rc);
CHECK_EQ(grid.getAlphaX().cols(), rc);
CHECK_EQ(grid.getAlphaY().rows(), rc);
CHECK_EQ(grid.getAlphaY().cols(), rc);
CHECK_EQ(grid.getDeltaRow(), 1);
CHECK_EQ(grid.getDeltaCol(), 1);
}
// CHECK_EQ(grid.getAlphaX().rows(), rc);
// CHECK_EQ(grid.getAlphaX().cols(), rc);
// CHECK_EQ(grid.getAlphaY().rows(), rc);
// CHECK_EQ(grid.getAlphaY().cols(), rc);
// CHECK_EQ(grid.getDeltaRow(), 1);
// CHECK_EQ(grid.getDeltaCol(), 1);
// }
SUBCASE("setting alphas") {
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
// SUBCASE("setting alphas") {
// // correct alpha matrices
// MatrixXd alphax = MatrixXd::Constant(rc, rc, 2);
// MatrixXd alphay = MatrixXd::Constant(rc, rc, 4);
// CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
grid.setAlpha(alphax, alphay);
CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_EQ(grid.getAlphaY(), alphay);
}
// grid.setAlpha(alphax, alphay);
// CHECK_EQ(grid.getAlphaX(), alphax);
// CHECK_EQ(grid.getAlphaY(), alphay);
// }
SUBCASE("setting domain") {
int dr = 8;
int dc = 9;
// SUBCASE("setting domain") {
// int dr = 8;
// int dc = 9;
// set 2D domain
CHECK_NOTHROW(grid.setDomain(dr, dc));
// // set 2D domain
// CHECK_NOTHROW(grid.setDomain(dr, dc));
grid.setDomain(dr, dc);
CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc));
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc));
}
}
// grid.setDomain(dr, dc);
// CHECK_EQ(grid.getDeltaCol(), double(dc) / double(rc));
// CHECK_EQ(grid.getDeltaRow(), double(dr) / double(rc));
// }
// }
TEST_CASE("2D Grid64 non-quadratic") {
int r = 12;
int c = 15;
Eigen::MatrixXd conc(r, c);
Grid64 grid(conc);
grid.initAlpha();
// TEST_CASE("2D Grid64 non-quadratic") {
// int r = 12;
// int c = 15;
// Eigen::MatrixXd conc(r, c);
// Grid64 grid(conc);
// grid.initAlpha();
SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 2);
CHECK_EQ(grid.getCol(), c);
CHECK_EQ(grid.getRow(), r);
// SUBCASE("correct construction") {
// CHECK_EQ(grid.getDim(), 2);
// CHECK_EQ(grid.getCol(), c);
// CHECK_EQ(grid.getRow(), r);
CHECK_EQ(grid.getConcentrations().rows(), r);
CHECK_EQ(grid.getConcentrations().cols(), c);
// CHECK_EQ(grid.getConcentrations().rows(), r);
// CHECK_EQ(grid.getConcentrations().cols(), c);
CHECK_EQ(grid.getAlphaX().rows(), r);
CHECK_EQ(grid.getAlphaX().cols(), c);
CHECK_EQ(grid.getAlphaY().rows(), r);
CHECK_EQ(grid.getAlphaY().cols(), c);
CHECK_EQ(grid.getDeltaRow(), 1);
CHECK_EQ(grid.getDeltaCol(), 1);
}
// CHECK_EQ(grid.getAlphaX().rows(), r);
// CHECK_EQ(grid.getAlphaX().cols(), c);
// CHECK_EQ(grid.getAlphaY().rows(), r);
// CHECK_EQ(grid.getAlphaY().cols(), c);
// CHECK_EQ(grid.getDeltaRow(), 1);
// CHECK_EQ(grid.getDeltaCol(), 1);
// }
SUBCASE("setting alphas") {
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(r, c, 2);
MatrixXd alphay = MatrixXd::Constant(r, c, 4);
CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
// SUBCASE("setting alphas") {
// // correct alpha matrices
// MatrixXd alphax = MatrixXd::Constant(r, c, 2);
// MatrixXd alphay = MatrixXd::Constant(r, c, 4);
// CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
grid.setAlpha(alphax, alphay);
CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_EQ(grid.getAlphaY(), alphay);
}
// grid.setAlpha(alphax, alphay);
// CHECK_EQ(grid.getAlphaX(), alphax);
// CHECK_EQ(grid.getAlphaY(), alphay);
// }
SUBCASE("setting domain") {
int dr = 8;
int dc = 9;
// SUBCASE("setting domain") {
// int dr = 8;
// int dc = 9;
// set 2D domain
CHECK_NOTHROW(grid.setDomain(dr, dc));
// // set 2D domain
// CHECK_NOTHROW(grid.setDomain(dr, dc));
grid.setDomain(dr, dc);
CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
}
}
// grid.setDomain(dr, dc);
// CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
// CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
// }
// }
TEST_CASE("2D Grid64 non-quadratic from pointer") {
int r = 4;
int c = 5;
std::vector<double> concentrations(r * c);
// TEST_CASE("2D Grid64 non-quadratic from pointer") {
// 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();
// for (int i = 0; i < r * c; i++) {
// concentrations[i] = i;
// }
// Grid64 grid(concentrations.data(), r, c);
// grid.initAlpha();
SUBCASE("correct construction") {
CHECK_EQ(grid.getDim(), 2);
CHECK_EQ(grid.getCol(), c);
CHECK_EQ(grid.getRow(), r);
// SUBCASE("correct construction") {
// CHECK_EQ(grid.getDim(), 2);
// CHECK_EQ(grid.getCol(), c);
// CHECK_EQ(grid.getRow(), r);
CHECK_EQ(grid.getConcentrations().rows(), r);
CHECK_EQ(grid.getConcentrations().cols(), c);
// CHECK_EQ(grid.getConcentrations().rows(), r);
// CHECK_EQ(grid.getConcentrations().cols(), c);
CHECK_EQ(grid.getAlphaX().rows(), r);
CHECK_EQ(grid.getAlphaX().cols(), c);
CHECK_EQ(grid.getAlphaY().rows(), r);
CHECK_EQ(grid.getAlphaY().cols(), c);
CHECK_EQ(grid.getDeltaRow(), 1);
CHECK_EQ(grid.getDeltaCol(), 1);
}
// CHECK_EQ(grid.getAlphaX().rows(), r);
// CHECK_EQ(grid.getAlphaX().cols(), c);
// CHECK_EQ(grid.getAlphaY().rows(), r);
// CHECK_EQ(grid.getAlphaY().cols(), c);
// CHECK_EQ(grid.getDeltaRow(), 1);
// CHECK_EQ(grid.getDeltaCol(), 1);
// }
SUBCASE("setting alphas") {
// correct alpha matrices
MatrixXd alphax = MatrixXd::Constant(r, c, 2);
MatrixXd alphay = MatrixXd::Constant(r, c, 4);
CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
// SUBCASE("setting alphas") {
// // correct alpha matrices
// MatrixXd alphax = MatrixXd::Constant(r, c, 2);
// MatrixXd alphay = MatrixXd::Constant(r, c, 4);
// CHECK_NOTHROW(grid.setAlpha(alphax, alphay));
grid.setAlpha(alphax, alphay);
CHECK_EQ(grid.getAlphaX(), alphax);
CHECK_EQ(grid.getAlphaY(), alphay);
}
// grid.setAlpha(alphax, alphay);
// CHECK_EQ(grid.getAlphaX(), alphax);
// CHECK_EQ(grid.getAlphaY(), alphay);
// }
SUBCASE("setting domain") {
int dr = 8;
int dc = 9;
// SUBCASE("setting domain") {
// int dr = 8;
// int dc = 9;
// set 2D domain
CHECK_NOTHROW(grid.setDomain(dr, dc));
// // set 2D domain
// CHECK_NOTHROW(grid.setDomain(dr, dc));
grid.setDomain(dr, dc);
CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
}
// grid.setDomain(dr, dc);
// CHECK_EQ(grid.getDeltaCol(), double(dc) / double(c));
// CHECK_EQ(grid.getDeltaRow(), double(dr) / double(r));
// }
SUBCASE("correct values") {
CHECK_EQ(grid.getConcentrations()(0, 0), 0);
CHECK_EQ(grid.getConcentrations()(0, 1), 1);
CHECK_EQ(grid.getConcentrations()(1, 0), c);
CHECK_EQ(grid.getConcentrations()(2, 1), 2 * c + 1);
}
}
// SUBCASE("correct values") {
// CHECK_EQ(grid.getConcentrations()(0, 0), 0);
// CHECK_EQ(grid.getConcentrations()(0, 1), 1);
// CHECK_EQ(grid.getConcentrations()(1, 0), c);
// CHECK_EQ(grid.getConcentrations()(2, 1), 2 * c + 1);
// }
// }