diff --git a/src/algebra/DenseMatrix.hpp b/src/algebra/DenseMatrix.hpp
index 952342fedbd04246879356d67425275d471af2c9..15967441b22069df2b809252f6ccba629facf23e 100644
--- a/src/algebra/DenseMatrix.hpp
+++ b/src/algebra/DenseMatrix.hpp
@@ -55,16 +55,14 @@ class DenseMatrix // LCOV_EXCL_LINE
return A_transpose;
}
- friend PUGS_INLINE DenseMatrix
- operator*(const DataType& a, const DenseMatrix& A)
+ friend PUGS_INLINE DenseMatrix operator*(const DataType& a, const DenseMatrix& A)
{
DenseMatrix<std::remove_const_t<DataType>> aA = copy(A);
return aA *= a;
}
template <typename DataType2>
- PUGS_INLINE Vector<std::remove_const_t<DataType>>
- operator*(const Vector<DataType2>& x) const
+ PUGS_INLINE Vector<std::remove_const_t<DataType>> operator*(const Vector<DataType2>& x) const
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
@@ -82,8 +80,7 @@ class DenseMatrix // LCOV_EXCL_LINE
}
template <typename DataType2>
- PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
- operator*(const DenseMatrix<DataType2>& B) const
+ PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>> operator*(const DenseMatrix<DataType2>& B) const
{
static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
"incompatible data types");
@@ -285,15 +282,14 @@ class DenseMatrix // LCOV_EXCL_LINE
explicit DenseMatrix(size_t nb_rows) noexcept : m_nb_rows{nb_rows}, m_nb_columns{nb_rows}, m_values{nb_rows * nb_rows}
{}
- template <size_t N>
- explicit DenseMatrix(const TinyMatrix<N, DataType>& M) noexcept : m_nb_rows{N}, m_nb_columns{N}, m_values{N * N}
+ template <size_t M, size_t N>
+ explicit DenseMatrix(const TinyMatrix<M, N, DataType>& A) noexcept : m_nb_rows{M}, m_nb_columns{N}, m_values{M * N}
{
- parallel_for(
- N, PUGS_LAMBDA(const index_type i) {
- for (size_t j = 0; j < N; ++j) {
- m_values[i * N + j] = M(i, j);
- }
- });
+ for (size_t i = 0; i < M; ++i) {
+ for (size_t j = 0; j < N; ++j) {
+ m_values[i * M + j] = A(i, j);
+ }
+ }
}
DenseMatrix() noexcept : m_nb_rows{0}, m_nb_columns{0} {}
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 98030dd447ee6cac14b6de0f3c45d55b68cdf7ba..73c187772126a478610aa2c4ee7313cc3333b41f 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -10,15 +10,14 @@
#include <iostream>
-template <size_t N, typename T = double>
+template <size_t M, size_t N = M, typename T = double>
class [[nodiscard]] TinyMatrix
{
public:
using data_type = T;
private:
- T m_values[N * N];
- static_assert((N > 0), "TinyMatrix size must be strictly positive");
+ T m_values[M * N];
PUGS_FORCEINLINE
constexpr size_t _index(size_t i, size_t j) const noexcept // LCOV_EXCL_LINE (due to forced inline)
@@ -29,7 +28,7 @@ class [[nodiscard]] TinyMatrix
template <typename... Args>
PUGS_FORCEINLINE constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
{
- m_values[N * N - 1 - sizeof...(args)] = t;
+ m_values[M * N - 1 - sizeof...(args)] = t;
if constexpr (sizeof...(args) > 0) {
this->_unpackVariadicInput(std::forward<Args>(args)...);
}
@@ -39,13 +38,13 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr bool isSquare() const noexcept
{
- return true;
+ return M == N;
}
PUGS_INLINE
constexpr size_t dimension() const
{
- return N * N;
+ return M * N;
}
PUGS_INLINE
@@ -57,7 +56,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr size_t numberOfRows() const
{
- return N;
+ return M;
}
PUGS_INLINE
@@ -70,7 +69,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix operator-() const
{
TinyMatrix opposite;
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
opposite.m_values[i] = -m_values[i];
}
return opposite;
@@ -92,19 +91,19 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix& operator*=(const T& t)
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] *= t;
}
return *this;
}
- PUGS_INLINE
- constexpr TinyMatrix operator*(const TinyMatrix& B) const
+ template <size_t P>
+ PUGS_INLINE constexpr TinyMatrix<M, P, T> operator*(const TinyMatrix<N, P, T>& B) const
{
const TinyMatrix& A = *this;
- TinyMatrix AB;
- for (size_t i = 0; i < N; ++i) {
- for (size_t j = 0; j < N; ++j) {
+ TinyMatrix<M, P, T> AB;
+ for (size_t i = 0; i < M; ++i) {
+ for (size_t j = 0; j < P; ++j) {
T sum = A(i, 0) * B(0, j);
for (size_t k = 1; k < N; ++k) {
sum += A(i, k) * B(k, j);
@@ -116,11 +115,11 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyVector<N, T> operator*(const TinyVector<N, T>& x) const
+ constexpr TinyVector<M, T> operator*(const TinyVector<N, T>& x) const
{
const TinyMatrix& A = *this;
- TinyVector<N, T> Ax;
- for (size_t i = 0; i < N; ++i) {
+ TinyVector<M, T> Ax;
+ for (size_t i = 0; i < M; ++i) {
T sum = A(i, 0) * x[0];
for (size_t j = 1; j < N; ++j) {
sum += A(i, j) * x[j];
@@ -134,7 +133,7 @@ class [[nodiscard]] TinyMatrix
constexpr friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
{
os << '[';
- for (size_t i = 0; i < N; ++i) {
+ for (size_t i = 0; i < M; ++i) {
os << '(' << NaNHelper(A(i, 0));
for (size_t j = 1; j < N; ++j) {
os << ',' << NaNHelper(A(i, j));
@@ -149,7 +148,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr bool operator==(const TinyMatrix& A) const
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
if (m_values[i] != A.m_values[i])
return false;
}
@@ -166,7 +165,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix operator+(const TinyMatrix& A) const
{
TinyMatrix sum;
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
sum.m_values[i] = m_values[i] + A.m_values[i];
}
return sum;
@@ -175,9 +174,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix operator+(TinyMatrix&& A) const
{
- for (size_t i = 0; i < N * N; ++i) {
- A.m_values[i] += m_values[i];
- }
+ A += *this;
return std::move(A);
}
@@ -185,7 +182,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix operator-(const TinyMatrix& A) const
{
TinyMatrix difference;
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
difference.m_values[i] = m_values[i] - A.m_values[i];
}
return difference;
@@ -194,7 +191,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix operator-(TinyMatrix&& A) const
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
A.m_values[i] = m_values[i] - A.m_values[i];
}
return std::move(A);
@@ -203,7 +200,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix& operator+=(const TinyMatrix& A)
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
}
return *this;
@@ -212,7 +209,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr void operator+=(const volatile TinyMatrix& A) volatile
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
}
}
@@ -220,7 +217,7 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr TinyMatrix& operator-=(const TinyMatrix& A)
{
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] -= A.m_values[i];
}
return *this;
@@ -229,14 +226,14 @@ class [[nodiscard]] TinyMatrix
PUGS_INLINE
constexpr T& operator()(size_t i, size_t j) noexcept(NO_ASSERT)
{
- Assert((i < N) and (j < N));
+ Assert((i < M) and (j < N));
return m_values[_index(i, j)];
}
PUGS_INLINE
constexpr const T& operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
{
- Assert((i < N) and (j < N));
+ Assert((i < M) and (j < N));
return m_values[_index(i, j)];
}
@@ -244,7 +241,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix& operator=(ZeroType) noexcept
{
static_assert(std::is_arithmetic<T>(), "Cannot assign 'zero' value for non-arithmetic types");
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] = 0;
}
return *this;
@@ -254,7 +251,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix& operator=(IdentityType) noexcept
{
static_assert(std::is_arithmetic<T>(), "Cannot assign 'identity' value for non-arithmetic types");
- for (size_t i = 0; i < N; ++i) {
+ for (size_t i = 0; i < M; ++i) {
for (size_t j = 0; j < N; ++j) {
m_values[_index(i, j)] = (i == j) ? 1 : 0;
}
@@ -271,7 +268,7 @@ class [[nodiscard]] TinyMatrix
template <typename... Args>
PUGS_INLINE constexpr TinyMatrix(const T& t, Args&&... args) noexcept
{
- static_assert(sizeof...(args) == N * N - 1, "wrong number of parameters");
+ static_assert(sizeof...(args) == M * N - 1, "wrong number of parameters");
this->_unpackVariadicInput(t, std::forward<Args>(args)...);
}
@@ -281,7 +278,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix() noexcept
{
#ifndef NDEBUG
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] = invalid_data_v<T>;
}
#endif // NDEBUG
@@ -292,7 +289,7 @@ class [[nodiscard]] TinyMatrix
{
static_assert(std::is_arithmetic<T>(), "Cannot construct from 'zero' value "
"for non-arithmetic types");
- for (size_t i = 0; i < N * N; ++i) {
+ for (size_t i = 0; i < M * N; ++i) {
m_values[i] = 0;
}
}
@@ -302,7 +299,7 @@ class [[nodiscard]] TinyMatrix
{
static_assert(std::is_arithmetic<T>(), "Cannot construct from 'identity' "
"value for non-arithmetic types");
- for (size_t i = 0; i < N; ++i) {
+ for (size_t i = 0; i < M; ++i) {
for (size_t j = 0; j < N; ++j) {
m_values[_index(i, j)] = (i == j) ? 1 : 0;
}
@@ -319,12 +316,12 @@ class [[nodiscard]] TinyMatrix
~TinyMatrix() = default;
};
-template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N, T>
-tensorProduct(const TinyVector<N, T>& x, const TinyVector<N, T>& y)
+template <size_t M, size_t N, typename T>
+PUGS_INLINE constexpr TinyMatrix<M, N, T>
+tensorProduct(const TinyVector<M, T>& x, const TinyVector<N, T>& y)
{
- TinyMatrix<N, T> A;
- for (size_t i = 0; i < N; ++i) {
+ TinyMatrix<M, N, T> A;
+ for (size_t i = 0; i < M; ++i) {
for (size_t j = 0; j < N; ++j) {
A(i, j) = x[i] * y[j];
}
@@ -334,16 +331,17 @@ tensorProduct(const TinyVector<N, T>& x, const TinyVector<N, T>& y)
template <size_t N, typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<N, T>& A)
+det(const TinyMatrix<N, N, T>& A)
{
- static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
- static_assert(std::is_floating_point<T>::value, "determinent for arbitrary dimension N is defined for floating "
+ static_assert(std::is_arithmetic<T>::value, "determinant is not defined for non-arithmetic types");
+ static_assert(std::is_floating_point<T>::value, "determinant for arbitrary dimension N is defined for floating "
"point types only");
- TinyMatrix<N, T> M = A;
+ TinyMatrix B = A;
TinyVector<N, size_t> index;
- for (size_t i = 0; i < N; ++i)
+ for (size_t i = 0; i < N; ++i) {
index[i] = i;
+ }
T determinent = 1;
for (size_t i = 0; i < N; ++i) {
@@ -351,7 +349,7 @@ det(const TinyMatrix<N, T>& A)
size_t l = j;
const size_t J = index[j];
for (size_t k = j + 1; k < N; ++k) {
- if (std::abs(M(index[k], i)) > std::abs(M(J, i))) {
+ if (std::abs(B(index[k], i)) > std::abs(B(J, i))) {
l = k;
}
}
@@ -361,14 +359,14 @@ det(const TinyMatrix<N, T>& A)
}
}
const size_t I = index[i];
- const T Mii = M(I, i);
- if (Mii != 0) {
- const T inv_Mii = 1. / M(I, i);
+ const T Bii = B(I, i);
+ if (Bii != 0) {
+ const T inv_Bii = 1. / B(I, i);
for (size_t k = i + 1; k < N; ++k) {
const size_t K = index[k];
- const T factor = M(K, i) * inv_Mii;
+ const T factor = B(K, i) * inv_Bii;
for (size_t l = i + 1; l < N; ++l) {
- M(K, l) -= factor * M(I, l);
+ B(K, l) -= factor * B(I, l);
}
}
} else {
@@ -377,14 +375,14 @@ det(const TinyMatrix<N, T>& A)
}
for (size_t i = 0; i < N; ++i) {
- determinent *= M(index[i], i);
+ determinent *= B(index[i], i);
}
return determinent;
}
template <typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<1, T>& A)
+det(const TinyMatrix<1, 1, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0, 0);
@@ -392,7 +390,7 @@ det(const TinyMatrix<1, T>& A)
template <typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<2, T>& A)
+det(const TinyMatrix<2, 2, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0, 0) * A(1, 1) - A(1, 0) * A(0, 1);
@@ -400,56 +398,56 @@ det(const TinyMatrix<2, T>& A)
template <typename T>
PUGS_INLINE constexpr T
-det(const TinyMatrix<3, T>& A)
+det(const TinyMatrix<3, 3, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0, 0) * (A(1, 1) * A(2, 2) - A(2, 1) * A(1, 2)) - A(1, 0) * (A(0, 1) * A(2, 2) - A(2, 1) * A(0, 2)) +
A(2, 0) * (A(0, 1) * A(1, 2) - A(1, 1) * A(0, 2));
}
-template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N - 1, T>
-getMinor(const TinyMatrix<N, T>& A, size_t I, size_t J)
+template <size_t M, size_t N, typename T>
+PUGS_INLINE constexpr TinyMatrix<M - 1, N - 1, T>
+getMinor(const TinyMatrix<M, N, T>& A, size_t I, size_t J)
{
- static_assert(N >= 2, "minor calculation requires at least 2x2 matrices");
- Assert((I < N) and (J < N));
- TinyMatrix<N - 1, T> M;
+ static_assert(M >= 2 and N >= 2, "minor calculation requires at least 2x2 matrices");
+ Assert((I < M) and (J < N));
+ TinyMatrix<M - 1, N - 1, T> minor;
for (size_t i = 0; i < I; ++i) {
for (size_t j = 0; j < J; ++j) {
- M(i, j) = A(i, j);
+ minor(i, j) = A(i, j);
}
for (size_t j = J + 1; j < N; ++j) {
- M(i, j - 1) = A(i, j);
+ minor(i, j - 1) = A(i, j);
}
}
- for (size_t i = I + 1; i < N; ++i) {
+ for (size_t i = I + 1; i < M; ++i) {
for (size_t j = 0; j < J; ++j) {
- M(i - 1, j) = A(i, j);
+ minor(i - 1, j) = A(i, j);
}
for (size_t j = J + 1; j < N; ++j) {
- M(i - 1, j - 1) = A(i, j);
+ minor(i - 1, j - 1) = A(i, j);
}
}
- return M;
+ return minor;
}
template <size_t N, typename T>
-PUGS_INLINE constexpr TinyMatrix<N, T> inverse(const TinyMatrix<N, T>& A);
+PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
template <typename T>
-PUGS_INLINE constexpr TinyMatrix<1, T>
-inverse(const TinyMatrix<1, T>& A)
+PUGS_INLINE constexpr TinyMatrix<1, 1, T>
+inverse(const TinyMatrix<1, 1, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
- TinyMatrix<1, T> A_1(1. / A(0, 0));
+ TinyMatrix<1, 1, T> A_1(1. / A(0, 0));
return A_1;
}
template <size_t N, typename T>
PUGS_INLINE constexpr T
-cofactor(const TinyMatrix<N, T>& A, size_t i, size_t j)
+cofactor(const TinyMatrix<N, N, T>& A, size_t i, size_t j)
{
static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non-arithmetic types");
const T sign = ((i + j) % 2) ? -1 : 1;
@@ -458,8 +456,8 @@ cofactor(const TinyMatrix<N, T>& A, size_t i, size_t j)
}
template <typename T>
-PUGS_INLINE constexpr TinyMatrix<2, T>
-inverse(const TinyMatrix<2, T>& A)
+PUGS_INLINE constexpr TinyMatrix<2, 2, T>
+inverse(const TinyMatrix<2, 2, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
@@ -467,22 +465,22 @@ inverse(const TinyMatrix<2, T>& A)
const T determinent = det(A);
const T inv_determinent = 1. / determinent;
- TinyMatrix<2, T> A_cofactors_T(A(1, 1), -A(0, 1), -A(1, 0), A(0, 0));
+ TinyMatrix<2, 2, T> A_cofactors_T(A(1, 1), -A(0, 1), -A(1, 0), A(0, 0));
return A_cofactors_T *= inv_determinent;
}
template <typename T>
-PUGS_INLINE constexpr TinyMatrix<3, T>
-inverse(const TinyMatrix<3, T>& A)
+PUGS_INLINE constexpr TinyMatrix<3, 3, T>
+inverse(const TinyMatrix<3, 3, T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
const T determinent = det(A);
- TinyMatrix<3, T> A_cofactors_T(cofactor(A, 0, 0), cofactor(A, 1, 0), cofactor(A, 2, 0), cofactor(A, 0, 1),
- cofactor(A, 1, 1), cofactor(A, 2, 1), cofactor(A, 0, 2), cofactor(A, 1, 2),
- cofactor(A, 2, 2));
+ TinyMatrix<3, 3, T> A_cofactors_T(cofactor(A, 0, 0), cofactor(A, 1, 0), cofactor(A, 2, 0), cofactor(A, 0, 1),
+ cofactor(A, 1, 1), cofactor(A, 2, 1), cofactor(A, 0, 2), cofactor(A, 1, 2),
+ cofactor(A, 2, 2));
return A_cofactors_T *= 1. / determinent;
}
diff --git a/src/output/OutputNamedItemValueSet.hpp b/src/output/OutputNamedItemValueSet.hpp
index bf90ddc3bad9636f2339695e87a446b11d81ae1a..69d6424d5bf3d4e9be9bbfb4ed3cee3281b3a786 100644
--- a/src/output/OutputNamedItemValueSet.hpp
+++ b/src/output/OutputNamedItemValueSet.hpp
@@ -67,9 +67,9 @@ class OutputNamedItemDataSet
NodeValue<const TinyVector<1, double>>,
NodeValue<const TinyVector<2, double>>,
NodeValue<const TinyVector<3, double>>,
- NodeValue<const TinyMatrix<1, double>>,
- NodeValue<const TinyMatrix<2, double>>,
- NodeValue<const TinyMatrix<3, double>>,
+ NodeValue<const TinyMatrix<1, 1, double>>,
+ NodeValue<const TinyMatrix<2, 2, double>>,
+ NodeValue<const TinyMatrix<3, 3, double>>,
CellValue<const bool>,
CellValue<const int>,
@@ -79,9 +79,9 @@ class OutputNamedItemDataSet
CellValue<const TinyVector<1, double>>,
CellValue<const TinyVector<2, double>>,
CellValue<const TinyVector<3, double>>,
- CellValue<const TinyMatrix<1, double>>,
- CellValue<const TinyMatrix<2, double>>,
- CellValue<const TinyMatrix<3, double>>,
+ CellValue<const TinyMatrix<1, 1, double>>,
+ CellValue<const TinyMatrix<2, 2, double>>,
+ CellValue<const TinyMatrix<3, 3, double>>,
NodeArray<const bool>,
NodeArray<const int>,
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 0bbdde090daab2cb49a086d359dbbe3b4e53739a..d01476cac3da67c1400efef7e3603e9bb91e096d 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -273,7 +273,7 @@ template <size_t N, typename DataType>
void
VTKWriter::_write_node_value(std::ofstream& os,
const std::string& name,
- const NodeValue<const TinyMatrix<N, DataType>>& item_value,
+ const NodeValue<const TinyMatrix<N, N, DataType>>& item_value,
SerializedDataList& serialized_data_list) const
{
os << "<DataArray type=\"" << VTKType<DataType>::name << "\" Name=\"" << name << "\" NumberOfComponents=\"" << N * N
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 7baede4556fb37e2b3b29f449dfbdaf5dfa8b143..913a05ef1de358cfe32f9b6faf1a9e6720585238 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -61,7 +61,7 @@ class VTKWriter : public WriterBase
template <size_t N, typename DataType>
void _write_node_value(std::ofstream& os,
const std::string& name,
- const NodeValue<const TinyMatrix<N, DataType>>& item_value,
+ const NodeValue<const TinyMatrix<N, N, DataType>>& item_value,
SerializedDataList& serialized_data_list) const;
template <typename DataType>
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index 21589e732ecde0a85491438ae7dba259b8ed3060..adf3f4a7f38216558602b5e5634c89364ef3cf6e 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -89,19 +89,19 @@ WriterBase::_register(const std::string& name,
switch (data_type.numberOfRows()) {
case 1: {
_register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, double>>&>(i_discrete_function),
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, 1, double>>&>(i_discrete_function),
named_item_data_set);
break;
}
case 2: {
_register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<2, double>>&>(i_discrete_function),
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<2, 2, double>>&>(i_discrete_function),
named_item_data_set);
break;
}
case 3: {
_register(name,
- dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<3, double>>&>(i_discrete_function),
+ dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<3, 3, double>>&>(i_discrete_function),
named_item_data_set);
break;
}
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index 504d694ad86a08811a4a87e2c8098dab77cd607d..46a2b1d6dc5923d8ee12668d88484875f404944c 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -12,7 +12,7 @@
template <size_t N, typename T>
class TinyVector;
-template <size_t N, typename T>
+template <size_t M, size_t N, typename T>
class TinyMatrix;
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
@@ -30,10 +30,10 @@ inline constexpr bool is_trivially_castable<TinyVector<N, T>> = is_trivially_cas
template <size_t N, typename T>
inline constexpr bool is_trivially_castable<const TinyVector<N, T>> = is_trivially_castable<T>;
-template <size_t N, typename T>
-inline constexpr bool is_trivially_castable<TinyMatrix<N, T>> = is_trivially_castable<T>;
-template <size_t N, typename T>
-inline constexpr bool is_trivially_castable<const TinyMatrix<N, T>> = is_trivially_castable<T>;
+template <size_t M, size_t N, typename T>
+inline constexpr bool is_trivially_castable<TinyMatrix<M, N, T>> = is_trivially_castable<T>;
+template <size_t M, size_t N, typename T>
+inline constexpr bool is_trivially_castable<const TinyMatrix<M, N, T>> = is_trivially_castable<T>;
// Traits is_false
@@ -105,8 +105,8 @@ inline constexpr bool is_tiny_vector_v<TinyVector<N, T>> = true;
template <typename T>
inline constexpr bool is_tiny_matrix_v = false;
-template <size_t N, typename T>
-inline constexpr bool is_tiny_matrix_v<TinyMatrix<N, T>> = true;
+template <size_t M, size_t N, typename T>
+inline constexpr bool is_tiny_matrix_v<TinyMatrix<M, N, T>> = true;
// Trais is ItemValue
diff --git a/tests/test_ASTNodeFunctionExpressionBuilder.cpp b/tests/test_ASTNodeFunctionExpressionBuilder.cpp
index b37dd10869bb5665c9292b22e363af09be4b9ead..4bf775a679490306efecc2617547cb424c6a256b 100644
--- a/tests/test_ASTNodeFunctionExpressionBuilder.cpp
+++ b/tests/test_ASTNodeFunctionExpressionBuilder.cpp
@@ -655,7 +655,7 @@ f(1);
std::string_view result = R"(
(root:ASTNodeListProcessor)
- `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<1ul, double>, ZeroType>)
+ `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<1ul, 1ul, double>, ZeroType>)
+-(language::name:f:NameProcessor)
`-(language::integer:1:ValueProcessor)
)";
@@ -672,7 +672,7 @@ f(1);
std::string_view result = R"(
(root:ASTNodeListProcessor)
- `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<2ul, double>, ZeroType>)
+ `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<2ul, 2ul, double>, ZeroType>)
+-(language::name:f:NameProcessor)
`-(language::integer:1:ValueProcessor)
)";
@@ -689,7 +689,7 @@ f(1);
std::string_view result = R"(
(root:ASTNodeListProcessor)
- `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<3ul, double>, ZeroType>)
+ `-(language::function_evaluation:FunctionExpressionProcessor<TinyMatrix<3ul, 3ul, double>, ZeroType>)
+-(language::name:f:NameProcessor)
`-(language::integer:1:ValueProcessor)
)";
diff --git a/tests/test_ASTNodeListAffectationExpressionBuilder.cpp b/tests/test_ASTNodeListAffectationExpressionBuilder.cpp
index 6eab14fe5162b64a1d883cd54d8b9a70e20a5d66..8d495053e0d0e09b67c103a4f387ecf6afc018bf 100644
--- a/tests/test_ASTNodeListAffectationExpressionBuilder.cpp
+++ b/tests/test_ASTNodeListAffectationExpressionBuilder.cpp
@@ -239,17 +239,17 @@ let (x1,x2,x3,x) : R^1x1*R^2x2*R^3x3*R,
std::string_view result = R"(
(root:ASTNodeListProcessor)
- +-(language::eq_op:AffectationProcessor<language::eq_op, TinyMatrix<1ul, double>, long>)
+ +-(language::eq_op:AffectationProcessor<language::eq_op, TinyMatrix<1ul, 1ul, double>, long>)
| +-(language::name:a:NameProcessor)
| `-(language::integer:0:ValueProcessor)
- +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<2ul, double> >)
+ +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<2ul, 2ul, double> >)
| +-(language::name:b:NameProcessor)
| `-(language::expression_list:ASTNodeExpressionListProcessor)
| +-(language::integer:1:ValueProcessor)
| +-(language::integer:2:ValueProcessor)
| +-(language::integer:3:ValueProcessor)
| `-(language::integer:4:ValueProcessor)
- +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<3ul, double> >)
+ +-(language::eq_op:AffectationToTinyMatrixFromListProcessor<language::eq_op, TinyMatrix<3ul, 3ul, double> >)
| +-(language::name:c:NameProcessor)
| `-(language::expression_list:ASTNodeExpressionListProcessor)
| +-(language::integer:9:ValueProcessor)
diff --git a/tests/test_AffectationToTupleProcessor.cpp b/tests/test_AffectationToTupleProcessor.cpp
index 650e60ca61bb6bc1d8b2e239c9e4c67713c3212a..d832ac2739cb6b8417e6dfd3818e11252639b7b4 100644
--- a/tests/test_AffectationToTupleProcessor.cpp
+++ b/tests/test_AffectationToTupleProcessor.cpp
@@ -87,7 +87,7 @@ let s :(R^1); s = 1.3;
const std::string A_string = []() -> std::string {
std::ostringstream os;
- os << TinyMatrix<2, double>{1, 2, 3, 4};
+ os << TinyMatrix<2>{1, 2, 3, 4};
return os.str();
}();
@@ -162,7 +162,7 @@ let t :(R^1); t = (x,2);
const std::string A_string = []() -> std::string {
std::ostringstream os;
- os << TinyMatrix<2, double>{1, 2, 3, 4};
+ os << TinyMatrix<2>{1, 2, 3, 4};
return os.str();
}();
@@ -216,7 +216,7 @@ let s :(string); s = x;
const std::string A_string = []() -> std::string {
std::ostringstream os;
- os << TinyMatrix<3, double>{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ os << TinyMatrix<3>{1, 2, 3, 4, 5, 6, 7, 8, 9};
return os.str();
}();
diff --git a/tests/test_ArrayUtils.cpp b/tests/test_ArrayUtils.cpp
index 79df610d909873a84328b59403ec861269b4f8b3..eb976b40df67c640bbe1a376671706835ca4c2e0 100644
--- a/tests/test_ArrayUtils.cpp
+++ b/tests/test_ArrayUtils.cpp
@@ -64,7 +64,7 @@ TEST_CASE("ArrayUtils", "[utils]")
SECTION("TinyMatrix Sum")
{
- using N22 = TinyMatrix<2, int>;
+ using N22 = TinyMatrix<2, 2, int>;
Array<N22> b(10);
b[0] = {13, 2, 0, 1};
b[1] = {1, 3, 6, 3};
diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index 1fefcdc6ea1e50e4c48fdff2be95d737a118de5a..8b7e48ea8f999b7d3d166a2f0b6eefc2a363bc75 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -302,9 +302,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
SECTION("matrix")
{
- const TinyMatrix<3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
- const TinyMatrix<3, size_t> zero{ZeroType{}};
- DiscreteFunctionP0<Dimension, TinyMatrix<3, size_t>> f{mesh};
+ const TinyMatrix<3, 3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ const TinyMatrix<3, 3, size_t> zero{ZeroType{}};
+ DiscreteFunctionP0<Dimension, TinyMatrix<3, 3, size_t>> f{mesh};
f.fill(value);
REQUIRE(all_values_equal(f, value));
@@ -318,7 +318,7 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
copy_to(g, f);
g.fill(zero);
- DiscreteFunctionP0<Dimension, const TinyMatrix<3, size_t>> h = copy(f);
+ DiscreteFunctionP0<Dimension, const TinyMatrix<3, 3, size_t>> h = copy(f);
REQUIRE(all_values_equal(f, value));
REQUIRE(all_values_equal(g, zero));
@@ -397,9 +397,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
SECTION("matrix")
{
- const TinyMatrix<3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
- const TinyMatrix<3, size_t> zero{ZeroType{}};
- DiscreteFunctionP0<Dimension, TinyMatrix<3, size_t>> f{mesh};
+ const TinyMatrix<3, 3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ const TinyMatrix<3, 3, size_t> zero{ZeroType{}};
+ DiscreteFunctionP0<Dimension, TinyMatrix<3, 3, size_t>> f{mesh};
f.fill(value);
REQUIRE(all_values_equal(f, value));
@@ -413,7 +413,7 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
copy_to(g, f);
g.fill(zero);
- DiscreteFunctionP0<Dimension, const TinyMatrix<3, size_t>> h = copy(f);
+ DiscreteFunctionP0<Dimension, const TinyMatrix<3, 3, size_t>> h = copy(f);
REQUIRE(all_values_equal(f, value));
REQUIRE(all_values_equal(g, zero));
@@ -492,9 +492,9 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
SECTION("matrix")
{
- const TinyMatrix<3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
- const TinyMatrix<3, size_t> zero{ZeroType{}};
- DiscreteFunctionP0<Dimension, TinyMatrix<3, size_t>> f{mesh};
+ const TinyMatrix<3, 3, size_t> value{1, 2, 3, 4, 5, 6, 7, 8, 9};
+ const TinyMatrix<3, 3, size_t> zero{ZeroType{}};
+ DiscreteFunctionP0<Dimension, TinyMatrix<3, 3, size_t>> f{mesh};
f.fill(value);
REQUIRE(all_values_equal(f, value));
@@ -508,7 +508,7 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
copy_to(g, f);
g.fill(zero);
- DiscreteFunctionP0<Dimension, const TinyMatrix<3, size_t>> h = copy(f);
+ DiscreteFunctionP0<Dimension, const TinyMatrix<3, 3, size_t>> h = copy(f);
REQUIRE(all_values_equal(f, value));
REQUIRE(all_values_equal(g, zero));
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index f50d2abc61f26bdbebbac4ba14ce884b352ec0aa..800f8497ece648bf81e75b3cf41b3e6adb1e76a2 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -12,39 +12,41 @@
#include <sstream>
// Instantiate to ensure full coverage is performed
-template class TinyMatrix<1, int>;
-template class TinyMatrix<2, int>;
-template class TinyMatrix<3, int>;
-template class TinyMatrix<4, double>;
+template class TinyMatrix<1, 1, int>;
+template class TinyMatrix<2, 2, int>;
+template class TinyMatrix<3, 4, int>;
+template class TinyMatrix<4, 4, double>;
// clazy:excludeall=non-pod-global-static
TEST_CASE("TinyMatrix", "[algebra]")
{
- TinyMatrix<3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
- REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(1, 0) == 4) and (A(1, 1) == 5) and
- (A(1, 2) == 6) and (A(2, 0) == 7) and (A(2, 1) == 8) and (A(2, 2) == 9)));
+ TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
+ REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(0, 3) == 4) and //
+ (A(1, 0) == 5) and (A(1, 1) == 6) and (A(1, 2) == 7) and (A(1, 3) == 8) and //
+ (A(2, 0) == 9) and (A(2, 1) == 10) and (A(2, 2) == 11) and (A(2, 3) == 12)));
- TinyMatrix<3, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7);
+ TinyMatrix<3, 4, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7, 1, 3, 6);
SECTION("checking for opposed matrix")
{
- const TinyMatrix<3, int> minus_A = -A;
- REQUIRE(((minus_A(0, 0) == -1) and (minus_A(0, 1) == -2) and (minus_A(0, 2) == -3) and (minus_A(1, 0) == -4) and
- (minus_A(1, 1) == -5) and (minus_A(1, 2) == -6) and (minus_A(2, 0) == -7) and (minus_A(2, 1) == -8) and
- (minus_A(2, 2) == -9)));
+ const TinyMatrix minus_A = -A;
+ REQUIRE(
+ ((minus_A(0, 0) == -1) and (minus_A(0, 1) == -2) and (minus_A(0, 2) == -3) and (minus_A(0, 3) == -4) and //
+ (minus_A(1, 0) == -5) and (minus_A(1, 1) == -6) and (minus_A(1, 2) == -7) and (minus_A(1, 3) == -8) and //
+ (minus_A(2, 0) == -9) and (minus_A(2, 1) == -10) and (minus_A(2, 2) == -11) and (minus_A(2, 3) == -12)));
}
SECTION("checking for equality and difference tests")
{
- const TinyMatrix<3, int> copy_A = A;
- REQUIRE(((copy_A(0, 0) == 1) and (copy_A(0, 1) == 2) and (copy_A(0, 2) == 3) and (copy_A(1, 0) == 4) and
- (copy_A(1, 1) == 5) and (copy_A(1, 2) == 6) and (copy_A(2, 0) == 7) and (copy_A(2, 1) == 8) and
- (copy_A(2, 2) == 9)));
+ const TinyMatrix copy_A = A;
+ REQUIRE(((copy_A(0, 0) == 1) and (copy_A(0, 1) == 2) and (copy_A(0, 2) == 3) and (copy_A(0, 3) == 4) and //
+ (copy_A(1, 0) == 5) and (copy_A(1, 1) == 6) and (copy_A(1, 2) == 7) and (copy_A(1, 3) == 8) and //
+ (copy_A(2, 0) == 9) and (copy_A(2, 1) == 10) and (copy_A(2, 2) == 11) and (copy_A(2, 3) == 12)));
REQUIRE(copy_A == A);
REQUIRE_FALSE(copy_A != A);
- TinyMatrix<3, int> affected_A;
+ TinyMatrix<3, 4, int> affected_A;
affected_A = A;
REQUIRE(affected_A == A);
REQUIRE_FALSE(affected_A != A);
@@ -55,101 +57,109 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for scalar left product")
{
- const int a = 2;
- const TinyMatrix<3, int> aA = a * A;
+ const int a = 2;
+ const TinyMatrix<3, 4, int> aA = a * A;
- REQUIRE(aA == TinyMatrix<3, int>(2, 4, 6, 8, 10, 12, 14, 16, 18));
+ REQUIRE(aA == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
}
SECTION("checking for scalar seft product")
{
- const int a = 2;
- TinyMatrix<3, int> copy_A = A;
+ const int a = 2;
- REQUIRE((copy_A *= a) == TinyMatrix<3, int>(2, 4, 6, 8, 10, 12, 14, 16, 18));
+ TinyMatrix copy_A = A;
+
+ REQUIRE((copy_A *= a) == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
}
SECTION("checking for null matrix management")
{
- TinyMatrix<3, int> Z = zero;
- REQUIRE(Z == TinyMatrix<3, int>(0, 0, 0, 0, 0, 0, 0, 0, 0));
+ TinyMatrix<4, 3, int> Z = zero;
+ REQUIRE(Z == TinyMatrix<4, 3, int>(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
- TinyMatrix<3, int> affected_Z;
+ TinyMatrix<4, 3, int> affected_Z;
affected_Z = zero;
REQUIRE(affected_Z == Z);
}
SECTION("checking for identity management")
{
- TinyMatrix<3, int> I = identity;
- REQUIRE(I == TinyMatrix<3, int>(1, 0, 0, 0, 1, 0, 0, 0, 1));
+ TinyMatrix<3, 3, int> I = identity;
+ REQUIRE(I == TinyMatrix<3, 3, int>(1, 0, 0, 0, 1, 0, 0, 0, 1));
- TinyMatrix<3, int> affected_I;
+ TinyMatrix<3, 3, int> affected_I;
affected_I = identity;
REQUIRE(affected_I == I);
}
SECTION("checking for matrices sum")
{
- REQUIRE(A + B == TinyMatrix<3, int>(7, 7, 6, 12, 39, 12, 42, 14, 16));
+ REQUIRE(A + B == TinyMatrix<3, 4, int>(7, 7, 6, 12, 39, 12, 42, 14, 16, 11, 14, 18));
- TinyMatrix<3, int> ApB = A;
+ TinyMatrix<3, 4, int> ApB = A;
ApB += B;
REQUIRE(ApB == A + B);
- TinyMatrix<3, int> Ap2B = A + 2 * B;
+ TinyMatrix<3, 4, int> Ap2B = A + 2 * B;
REQUIRE(Ap2B == ApB + B);
}
SECTION("checking for matrices difference ")
{
- REQUIRE(A - B == TinyMatrix<3, int>(-5, -3, 0, -4, -29, 0, -28, 2, 2));
+ REQUIRE(A - B == TinyMatrix<3, 4, int>(-5, -3, 0, -4, -29, 0, -28, 2, 2, 9, 8, 6));
- TinyMatrix<3, int> AmB = A;
+ TinyMatrix<3, 4, int> AmB = A;
AmB -= B;
REQUIRE(AmB == A - B);
- TinyMatrix<3, int> Am2B = A - 2 * B;
+ TinyMatrix<3, 4, int> Am2B = A - 2 * B;
REQUIRE(Am2B == AmB - B);
}
SECTION("checking for matrices product")
{
- REQUIRE(A * B == TinyMatrix<3, int>(127, 91, 36, 274, 226, 84, 421, 361, 132));
+ TinyMatrix<4, 2, int> C{3, -2, 2, 6, -2, 5, 7, 2};
+ REQUIRE(A * C == TinyMatrix<3, 2, int>(29, 33, 69, 77, 109, 121));
+
+ TinyMatrix<2, 3, int> D{-3, 2, 3, 2, 5, -7};
+ REQUIRE(D * A == TinyMatrix<2, 4, int>(34, 36, 38, 40, -36, -36, -36, -36));
}
SECTION("checking for matrix-vector product")
{
- REQUIRE(A * TinyVector<3, int>(2, -3, 5) == TinyVector<3, int>(11, 23, 35));
+ REQUIRE(A * TinyVector<4, int>(2, -3, 5, 2) == TinyVector<3, int>(19, 43, 67));
}
SECTION("checking for tensor product")
{
- const TinyVector<3, int> u(1, 3, 7);
+ const TinyVector<4, int> u(1, 3, 7, 5);
const TinyVector<3, int> v(6, 2, -3);
- REQUIRE(tensorProduct(u, v) == TinyMatrix<3, int>(6, 2, -3, 18, 6, -9, 42, 14, -21));
+ REQUIRE(tensorProduct(u, v) == TinyMatrix<4, 3, int>(6, 2, -3, 18, 6, -9, 42, 14, -21, 30, 10, -15));
}
SECTION("checking for minor calculation")
{
- TinyMatrix<3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
- REQUIRE(getMinor(A, 0, 0) == TinyMatrix<2, int>(5, 6, 8, 9));
- REQUIRE(getMinor(A, 1, 0) == TinyMatrix<2, int>(2, 3, 8, 9));
- REQUIRE(getMinor(A, 2, 0) == TinyMatrix<2, int>(2, 3, 5, 6));
-
- REQUIRE(getMinor(A, 0, 1) == TinyMatrix<2, int>(4, 6, 7, 9));
- REQUIRE(getMinor(A, 1, 1) == TinyMatrix<2, int>(1, 3, 7, 9));
- REQUIRE(getMinor(A, 2, 1) == TinyMatrix<2, int>(1, 3, 4, 6));
-
- REQUIRE(getMinor(A, 0, 2) == TinyMatrix<2, int>(4, 5, 7, 8));
- REQUIRE(getMinor(A, 1, 2) == TinyMatrix<2, int>(1, 2, 7, 8));
- REQUIRE(getMinor(A, 2, 2) == TinyMatrix<2, int>(1, 2, 4, 5));
+ REQUIRE(getMinor(A, 0, 0) == TinyMatrix<2, 3, int>(6, 7, 8, 10, 11, 12));
+ REQUIRE(getMinor(A, 1, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 10, 11, 12));
+ REQUIRE(getMinor(A, 2, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 6, 7, 8));
+
+ REQUIRE(getMinor(A, 0, 1) == TinyMatrix<2, 3, int>(5, 7, 8, 9, 11, 12));
+ REQUIRE(getMinor(A, 1, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 9, 11, 12));
+ REQUIRE(getMinor(A, 2, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 5, 7, 8));
+
+ REQUIRE(getMinor(A, 0, 2) == TinyMatrix<2, 3, int>(5, 6, 8, 9, 10, 12));
+ REQUIRE(getMinor(A, 1, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 9, 10, 12));
+ REQUIRE(getMinor(A, 2, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 5, 6, 8));
+
+ REQUIRE(getMinor(A, 0, 3) == TinyMatrix<2, 3, int>(5, 6, 7, 9, 10, 11));
+ REQUIRE(getMinor(A, 1, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 9, 10, 11));
+ REQUIRE(getMinor(A, 2, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 5, 6, 7));
}
SECTION("checking for cofactors")
{
- TinyMatrix<3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
+ TinyMatrix<3, 3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
REQUIRE(cofactor(A, 0, 0) == (5 * 9 - 8 * 6));
REQUIRE(cofactor(A, 1, 0) == -(2 * 9 - 8 * 3));
REQUIRE(cofactor(A, 2, 0) == (2 * 6 - 5 * 3));
@@ -165,26 +175,25 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for determinant calculations")
{
- REQUIRE(det(TinyMatrix<1, int>(6)) == 6);
- REQUIRE(det(TinyMatrix<2, int>(3, 1, -3, 6)) == 21);
- REQUIRE(det(TinyMatrix<3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1);
- REQUIRE(det(B) == -1444);
- REQUIRE(det(TinyMatrix<4, double>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
+ REQUIRE(det(TinyMatrix<1, 1, int>(6)) == 6);
+ REQUIRE(det(TinyMatrix<2, 2, int>(3, 1, -3, 6)) == 21);
+ REQUIRE(det(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1);
+ REQUIRE(det(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == -1444);
+ REQUIRE(det(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
Catch::Approx(6661.455).epsilon(1E-14));
-
- REQUIRE(det(TinyMatrix<4, double>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
+ REQUIRE(det(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
}
SECTION("checking for inverse calculations")
{
{
- const TinyMatrix<1, double> A1(2);
+ const TinyMatrix<1> A1(2);
REQUIRE(inverse(A1)(0, 0) == Catch::Approx(0.5).epsilon(1E-14));
}
{
- const TinyMatrix<2, double> A2(2, 3, 4, 1);
- const TinyMatrix<2, double> I = inverse(A2) * A2;
+ const TinyMatrix<2> A2(2, 3, 4, 1);
+ const TinyMatrix<2> I = inverse(A2) * A2;
REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
@@ -193,8 +202,8 @@ TEST_CASE("TinyMatrix", "[algebra]")
}
{
- const TinyMatrix<3, double> A3(2, 3, 1, 4, -1, 5, -2, 3, 4);
- const TinyMatrix<3, double> I = inverse(A3) * A3;
+ const TinyMatrix<3> A3(2, 3, 1, 4, -1, 5, -2, 3, 4);
+ const TinyMatrix<3> I = inverse(A3) * A3;
REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
@@ -223,6 +232,11 @@ TEST_CASE("TinyMatrix", "[algebra]")
REQUIRE(TinyMatrix<3>{}.numberOfColumns() == 3);
REQUIRE(TinyMatrix<3>{}.dimension() == 9);
REQUIRE(TinyMatrix<3>{}.numberOfValues() == 9);
+
+ REQUIRE(TinyMatrix<3, 4>{}.numberOfRows() == 3);
+ REQUIRE(TinyMatrix<3, 4>{}.numberOfColumns() == 4);
+ REQUIRE(TinyMatrix<3, 4>{}.dimension() == 12);
+ REQUIRE(TinyMatrix<3, 4>{}.numberOfValues() == 12);
}
SECTION("is square matrix")
@@ -230,26 +244,28 @@ TEST_CASE("TinyMatrix", "[algebra]")
REQUIRE(TinyMatrix<1>{}.isSquare());
REQUIRE(TinyMatrix<2>{}.isSquare());
REQUIRE(TinyMatrix<3>{}.isSquare());
- REQUIRE(TinyMatrix<4>{}.isSquare());
+
+ REQUIRE_FALSE(TinyMatrix<3, 4>{}.isSquare());
}
SECTION("checking for matrices output")
{
- REQUIRE(Catch::Detail::stringify(A) == "[(1,2,3)(4,5,6)(7,8,9)]");
- REQUIRE(Catch::Detail::stringify(TinyMatrix<1, int>(7)) == "[(7)]");
+ REQUIRE(Catch::Detail::stringify(A) == "[(1,2,3,4)(5,6,7,8)(9,10,11,12)]");
+ REQUIRE(Catch::Detail::stringify(TinyMatrix<1, 1, int>(7)) == "[(7)]");
}
#ifndef NDEBUG
SECTION("output with signaling NaN")
{
- TinyMatrix<2> A;
+ TinyMatrix<2, 3> A;
A(0, 0) = 1;
+ A(0, 2) = 3;
A(1, 0) = 2;
std::ostringstream A_ost;
A_ost << A;
std::ostringstream ref_ost;
- ref_ost << "[(1,nan)(2,nan)]";
+ ref_ost << "[(1,nan,3)(2,nan,nan)]";
REQUIRE(A_ost.str() == ref_ost.str());
}
@@ -257,19 +273,19 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for bounds violation")
{
REQUIRE_THROWS_AS(A(3, 0), AssertError);
- REQUIRE_THROWS_AS(A(0, 3), AssertError);
+ REQUIRE_THROWS_AS(A(0, 4), AssertError);
REQUIRE_THROWS_AS(getMinor(A, 3, 0), AssertError);
- REQUIRE_THROWS_AS(getMinor(A, 0, 3), AssertError);
+ REQUIRE_THROWS_AS(getMinor(A, 0, 4), AssertError);
- const TinyMatrix<3, int>& constA = A;
+ const TinyMatrix<3, 4, int>& constA = A;
REQUIRE_THROWS_AS(constA(3, 0), AssertError);
- REQUIRE_THROWS_AS(constA(0, 3), AssertError);
+ REQUIRE_THROWS_AS(constA(0, 4), AssertError);
}
SECTION("checking for nan initialization")
{
- TinyMatrix<3, double> B;
+ TinyMatrix<3, 4, double> B;
for (size_t i = 0; i < B.numberOfRows(); ++i) {
for (size_t j = 0; j < B.numberOfColumns(); ++j) {
@@ -280,7 +296,7 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for bad initialization")
{
- TinyMatrix<3, int> B;
+ TinyMatrix<3, 4, int> B;
for (size_t i = 0; i < B.numberOfRows(); ++i) {
for (size_t j = 0; j < B.numberOfColumns(); ++j) {