diff --git a/doc/userdoc.org b/doc/userdoc.org
index 320a9c183270985d0d3ad17fbda9ebd562cad2c9..c338582d59ef05c6b2241c09b615a515e56b2fe0 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -2589,7 +2589,8 @@ The output is
#+RESULTS: dot-examples
A set of ~det~ functions is defined to get the determinant of matrices
-$\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$.
+of $\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and
+$\mathbb{R}^{3\times3}$.
#+NAME: det-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
@@ -2601,20 +2602,47 @@ $\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$.
The output is
#+RESULTS: det-examples
+The ~trace~ functions compute the trace of matrices of
+$\mathbb{R}^{1\times1}$, $\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$.
+#+NAME: trace-examples
+#+BEGIN_SRC pugs :exports both :results output
+ import math;
+ cout << "trace([[1.2]]) = " << trace([[1.2]]) << "\n";
+ cout << "trace([[1,2],[3,4]]) = " << trace([[1,2],[3,4]]) << "\n";
+ cout << "trace([[1,2,3],[4,5,6],[7,8,9]]) = "
+ << trace([[1,2,3],[4,5,6],[7,8,9]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: trace-examples
+
Also, one can compute inverses of $\mathbb{R}^{1\times1}$,
$\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ matrices using the
~inverse~ function set.
#+NAME: inverse-examples
#+BEGIN_SRC pugs :exports both :results output
import math;
- cout << "det([[1.2]]) = " << inverse([[1.2]]) << "\n";
- cout << "det([[1,2],[3,4]]) = " << inverse([[1,2],[3,4]]) << "\n";
- cout << "det([[3,2,1],[5,6,4],[7,8,9]]) = "
- << det([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
+ cout << "inverse([[1.2]]) = " << inverse([[1.2]]) << "\n";
+ cout << "inverse([[1,2],[3,4]]) = " << inverse([[1,2],[3,4]]) << "\n";
+ cout << "inverse([[3,2,1],[5,6,4],[7,8,9]]) = "
+ << inverse([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
#+END_SRC
The output is
#+RESULTS: inverse-examples
+Transpose of matrices of $\mathbb{R}^{1\times1}$,
+$\mathbb{R}^{2\times2}$ and $\mathbb{R}^{3\times3}$ are obtained using the
+~transpose~ functions.
+#+NAME: transpose-examples
+#+BEGIN_SRC pugs :exports both :results output
+ import math;
+ cout << "transpose([[1.2]]) = " << transpose([[1.2]]) << "\n";
+ cout << "transpose([[1,2],[3,4]]) = " << transpose([[1,2],[3,4]]) << "\n";
+ cout << "transpose([[3,2,1],[5,6,4],[7,8,9]]) = "
+ << transpose([[3,2,1],[5,6,4],[7,8,9]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: transpose-examples
+
#+BEGIN_note
Observe that the use of a proper rounding or truncation function is
the right way to convert a real value to an integer one. Available
@@ -3246,18 +3274,20 @@ Here is the list of these functions
- ~tan: Vh -> Vh~
- ~tanh: Vh -> Vh~
-The ~det~ functions are defined for $\mathbb{P}_0(\mathbb{R^{}}^{1\times1})$,
+A few functions are defined for $\mathbb{P}_0(\mathbb{R^{}}^{1\times1})$,
$\mathbb{P}_0(\mathbb{R^{}}^{2\times2})$ and
$\mathbb{P}_0(\mathbb{R^{}}^{3\times3})$ data and the return value is a
$\mathbb{P}_0(\mathbb{R})$ function. The value is simply the
application of the function to the cell values.
- ~det: Vh -> Vh~
+- ~trace: Vh -> Vh~
-The ~inverse~ functions are defined for $\mathbb{P}_0(\mathbb{R}^{d\times d})$
+Also, functions are defined for $\mathbb{P}_0(\mathbb{R}^{d\times d})$
data and the return value is a $\mathbb{P}_0(\mathbb{R}^{d\times d})$
function, with $d\in\{1,2,3\}$. The value is simply the application of
the function to the cell values.
- ~inverse: Vh -> Vh~
+- ~transpose: Vh -> Vh~
****** ~Vh*Vh -> Vh~
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index e34f98ed8c25df7a77ba90bbd7cd05305b5eb115..1717b48c78deaedc7cabd9f60feb215e29a80ebe 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -21,16 +21,21 @@ class [[nodiscard]] TinyMatrix
using data_type = T;
private:
+ static_assert((M > 0), "TinyMatrix number of rows must be strictly positive");
+ static_assert((N > 0), "TinyMatrix number of columns 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)
+ constexpr size_t
+ _index(size_t i, size_t j) const noexcept // LCOV_EXCL_LINE (due to forced inline)
{
return i * N + j;
}
template <typename... Args>
- PUGS_FORCEINLINE constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
+ PUGS_FORCEINLINE constexpr void
+ _unpackVariadicInput(const T& t, Args&&... args) noexcept
{
m_values[M * N - 1 - sizeof...(args)] = t;
if constexpr (sizeof...(args) > 0) {
@@ -40,13 +45,15 @@ class [[nodiscard]] TinyMatrix
public:
PUGS_INLINE
- constexpr bool isSquare() const noexcept
+ constexpr bool
+ isSquare() const noexcept
{
return M == N;
}
PUGS_INLINE
- constexpr friend TinyMatrix<N, M, T> transpose(const TinyMatrix& A)
+ constexpr friend TinyMatrix<N, M, T>
+ transpose(const TinyMatrix& A)
{
TinyMatrix<N, M, T> tA;
for (size_t i = 0; i < M; ++i) {
@@ -58,31 +65,36 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr size_t dimension() const
+ constexpr size_t
+ dimension() const
{
return M * N;
}
PUGS_INLINE
- constexpr size_t numberOfValues() const
+ constexpr size_t
+ numberOfValues() const
{
return this->dimension();
}
PUGS_INLINE
- constexpr size_t numberOfRows() const
+ constexpr size_t
+ numberOfRows() const
{
return M;
}
PUGS_INLINE
- constexpr size_t numberOfColumns() const
+ constexpr size_t
+ numberOfColumns() const
{
return N;
}
PUGS_INLINE
- constexpr TinyMatrix operator-() const
+ constexpr TinyMatrix
+ operator-() const
{
TinyMatrix opposite;
for (size_t i = 0; i < M * N; ++i) {
@@ -92,20 +104,23 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr friend TinyMatrix operator*(const T& t, const TinyMatrix& A)
+ constexpr friend TinyMatrix
+ operator*(const T& t, const TinyMatrix& A)
{
TinyMatrix B = A;
return B *= t;
}
PUGS_INLINE
- constexpr friend TinyMatrix operator*(const T& t, TinyMatrix&& A)
+ constexpr friend TinyMatrix
+ operator*(const T& t, TinyMatrix&& A)
{
return std::move(A *= t);
}
PUGS_INLINE
- constexpr TinyMatrix& operator*=(const T& t)
+ constexpr TinyMatrix&
+ operator*=(const T& t)
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] *= t;
@@ -114,7 +129,8 @@ class [[nodiscard]] TinyMatrix
}
template <size_t P>
- PUGS_INLINE constexpr TinyMatrix<M, P, T> operator*(const TinyMatrix<N, P, T>& B) const
+ PUGS_INLINE constexpr TinyMatrix<M, P, T>
+ operator*(const TinyMatrix<N, P, T>& B) const
{
const TinyMatrix& A = *this;
TinyMatrix<M, P, T> AB;
@@ -131,7 +147,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyVector<M, T> operator*(const TinyVector<N, T>& x) const
+ constexpr TinyVector<M, T>
+ operator*(const TinyVector<N, T>& x) const
{
const TinyMatrix& A = *this;
TinyVector<M, T> Ax;
@@ -146,7 +163,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
+ constexpr friend std::ostream&
+ operator<<(std::ostream& os, const TinyMatrix& A)
{
os << '[';
for (size_t i = 0; i < M; ++i) {
@@ -165,7 +183,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr bool operator==(const TinyMatrix& A) const
+ constexpr bool
+ operator==(const TinyMatrix& A) const
{
for (size_t i = 0; i < M * N; ++i) {
if (m_values[i] != A.m_values[i])
@@ -175,13 +194,15 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr bool operator!=(const TinyMatrix& A) const
+ constexpr bool
+ operator!=(const TinyMatrix& A) const
{
return not this->operator==(A);
}
PUGS_INLINE
- constexpr TinyMatrix operator+(const TinyMatrix& A) const
+ constexpr TinyMatrix
+ operator+(const TinyMatrix& A) const
{
TinyMatrix sum;
for (size_t i = 0; i < M * N; ++i) {
@@ -191,14 +212,16 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix operator+(TinyMatrix&& A) const
+ constexpr TinyMatrix
+ operator+(TinyMatrix&& A) const
{
A += *this;
return std::move(A);
}
PUGS_INLINE
- constexpr TinyMatrix operator-(const TinyMatrix& A) const
+ constexpr TinyMatrix
+ operator-(const TinyMatrix& A) const
{
TinyMatrix difference;
for (size_t i = 0; i < M * N; ++i) {
@@ -208,7 +231,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix operator-(TinyMatrix&& A) const
+ constexpr TinyMatrix
+ operator-(TinyMatrix&& A) const
{
for (size_t i = 0; i < M * N; ++i) {
A.m_values[i] = m_values[i] - A.m_values[i];
@@ -217,7 +241,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix& operator+=(const TinyMatrix& A)
+ constexpr TinyMatrix&
+ operator+=(const TinyMatrix& A)
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
@@ -226,7 +251,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr void operator+=(const volatile TinyMatrix& A) volatile
+ constexpr void
+ operator+=(const volatile TinyMatrix& A) volatile
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
@@ -234,7 +260,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix& operator-=(const TinyMatrix& A)
+ constexpr TinyMatrix&
+ operator-=(const TinyMatrix& A)
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] -= A.m_values[i];
@@ -243,21 +270,24 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr T& operator()(size_t i, size_t j) noexcept(NO_ASSERT)
+ constexpr T&
+ operator()(size_t i, size_t j) noexcept(NO_ASSERT)
{
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)
+ constexpr const T&
+ operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
{
Assert((i < M) and (j < N));
return m_values[_index(i, j)];
}
PUGS_INLINE
- constexpr TinyMatrix& operator=(ZeroType) noexcept
+ 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 < M * N; ++i) {
@@ -267,7 +297,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix& operator=(IdentityType) noexcept
+ 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 < M; ++i) {
@@ -329,7 +360,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix(const TinyMatrix&) noexcept = default;
PUGS_INLINE
- TinyMatrix(TinyMatrix && A) noexcept = default;
+ TinyMatrix(TinyMatrix&& A) noexcept = default;
PUGS_INLINE
~TinyMatrix() = default;
@@ -450,6 +481,19 @@ getMinor(const TinyMatrix<M, N, T>& A, size_t I, size_t J)
return m;
}
+template <size_t N, typename T>
+PUGS_INLINE T
+trace(const TinyMatrix<N, N, T>& A)
+{
+ static_assert(std::is_arithmetic<T>::value, "trace is not defined for non-arithmetic types");
+
+ T t = A(0, 0);
+ for (size_t i = 1; i < N; ++i) {
+ t += A(i, i);
+ }
+ return t;
+}
+
template <size_t N, typename T>
PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index 5f2de33d8502fa84edd809f2dfbfa51c91c12ef4..4a1e7b35b76a4d9b8b05ebb1355f3cdca4f751d9 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -77,6 +77,12 @@ MathModule::MathModule()
this->_addBuiltinFunction("det", std::function([](const TinyMatrix<3>& A) -> double { return det(A); }));
+ this->_addBuiltinFunction("trace", std::function([](const TinyMatrix<1> A) -> double { return trace(A); }));
+
+ this->_addBuiltinFunction("trace", std::function([](const TinyMatrix<2>& A) -> double { return trace(A); }));
+
+ this->_addBuiltinFunction("trace", std::function([](const TinyMatrix<3>& A) -> double { return trace(A); }));
+
this->_addBuiltinFunction("inverse",
std::function([](const TinyMatrix<1> A) -> TinyMatrix<1> { return inverse(A); }));
@@ -85,6 +91,15 @@ MathModule::MathModule()
this->_addBuiltinFunction("inverse",
std::function([](const TinyMatrix<3>& A) -> TinyMatrix<3> { return inverse(A); }));
+
+ this->_addBuiltinFunction("transpose",
+ std::function([](const TinyMatrix<1> A) -> TinyMatrix<1> { return transpose(A); }));
+
+ this->_addBuiltinFunction("transpose",
+ std::function([](const TinyMatrix<2>& A) -> TinyMatrix<2> { return transpose(A); }));
+
+ this->_addBuiltinFunction("transpose",
+ std::function([](const TinyMatrix<3>& A) -> TinyMatrix<3> { return transpose(A); }));
}
void
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index e4ac44a3328d40415a4cbdca52600e224f6885cc..67622f4a7c2ebbc7284f7305be98aeb3beb8768f 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -554,6 +554,57 @@ det(const std::shared_ptr<const IDiscreteFunction>& A)
}
}
+template <typename MatrixType>
+std::shared_ptr<const IDiscreteFunction>
+trace(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+ switch (A->mesh()->dimension()) {
+ case 1: {
+ using DiscreteFunctionType = DiscreteFunctionP0<1, MatrixType>;
+ return std::make_shared<const DiscreteFunctionP0<1, double>>(trace(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 2: {
+ using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+ return std::make_shared<const DiscreteFunctionP0<2, double>>(trace(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 3: {
+ using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+ return std::make_shared<const DiscreteFunctionP0<3, double>>(trace(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+trace(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+ if (A->dataType() == ASTNodeDataType::matrix_t and A->descriptor().type() == DiscreteFunctionType::P0) {
+ if (A->dataType().numberOfRows() != A->dataType().numberOfColumns()) {
+ throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+ }
+ switch (A->dataType().numberOfRows()) {
+ case 1: {
+ return trace<TinyMatrix<1>>(A);
+ }
+ case 2: {
+ return trace<TinyMatrix<2>>(A);
+ }
+ case 3: {
+ return trace<TinyMatrix<3>>(A);
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("invalid matrix type");
+ }
+ // LCOV_EXCL_STOP
+ }
+ } else {
+ throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+ }
+}
+
template <typename MatrixType>
std::shared_ptr<const IDiscreteFunction>
inverse(const std::shared_ptr<const IDiscreteFunction>& A)
@@ -605,6 +656,57 @@ inverse(const std::shared_ptr<const IDiscreteFunction>& A)
}
}
+template <typename MatrixType>
+std::shared_ptr<const IDiscreteFunction>
+transpose(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+ switch (A->mesh()->dimension()) {
+ case 1: {
+ using DiscreteFunctionType = DiscreteFunctionP0<1, MatrixType>;
+ return std::make_shared<const DiscreteFunctionType>(transpose(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 2: {
+ using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+ return std::make_shared<const DiscreteFunctionType>(transpose(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 3: {
+ using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+ return std::make_shared<const DiscreteFunctionType>(transpose(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+transpose(const std::shared_ptr<const IDiscreteFunction>& A)
+{
+ if (A->dataType() == ASTNodeDataType::matrix_t and A->descriptor().type() == DiscreteFunctionType::P0) {
+ if (A->dataType().numberOfRows() != A->dataType().numberOfColumns()) {
+ throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+ }
+ switch (A->dataType().numberOfRows()) {
+ case 1: {
+ return transpose<TinyMatrix<1>>(A);
+ }
+ case 2: {
+ return transpose<TinyMatrix<2>>(A);
+ }
+ case 3: {
+ return transpose<TinyMatrix<3>>(A);
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("invalid matrix type");
+ }
+ // LCOV_EXCL_STOP
+ }
+ } else {
+ throw NormalError(EmbeddedIDiscreteFunctionUtils::invalidOperandType(A));
+ }
+}
+
double
min(const std::shared_ptr<const IDiscreteFunction>& f)
{
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
index 1908dd30cf16c6cc7fa3c3281d42429b0c661580..24994c743428d98befb8d34cc18c2cbbf12a2ba4 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -67,8 +67,12 @@ std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&,
std::shared_ptr<const IDiscreteFunction> det(const std::shared_ptr<const IDiscreteFunction>&);
+std::shared_ptr<const IDiscreteFunction> trace(const std::shared_ptr<const IDiscreteFunction>&);
+
std::shared_ptr<const IDiscreteFunction> inverse(const std::shared_ptr<const IDiscreteFunction>&);
+std::shared_ptr<const IDiscreteFunction> transpose(const std::shared_ptr<const IDiscreteFunction>&);
+
std::shared_ptr<const IDiscreteFunction> sum_of_Vh_components(const std::shared_ptr<const IDiscreteFunction>&);
std::shared_ptr<const IDiscreteFunction> vectorize(
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 6e1ef6a0d28a3436543fa090bee813586874bbee..d1bd72c3394346db10898cd83f1c609974f17697 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -536,6 +536,18 @@ class DiscreteFunctionP0 : public IDiscreteFunction
return result;
}
+ PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+ trace(const DiscreteFunctionP0& A)
+ {
+ static_assert(is_tiny_matrix_v<DataType>);
+ Assert(A.m_cell_values.isBuilt());
+ DiscreteFunctionP0<Dimension, double> result{A.m_mesh};
+ parallel_for(
+ A.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = trace(A[cell_id]); });
+
+ return result;
+ }
+
PUGS_INLINE friend DiscreteFunctionP0
inverse(const DiscreteFunctionP0& A)
{
@@ -549,6 +561,19 @@ class DiscreteFunctionP0 : public IDiscreteFunction
return result;
}
+ PUGS_INLINE friend DiscreteFunctionP0
+ transpose(const DiscreteFunctionP0& A)
+ {
+ static_assert(is_tiny_matrix_v<DataType>);
+ static_assert(DataType::NumberOfRows == DataType::NumberOfColumns, "cannot compute inverse of non square matrices");
+ Assert(A.m_cell_values.isBuilt());
+ DiscreteFunctionP0 result{A.m_mesh};
+ parallel_for(
+ A.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = transpose(A[cell_id]); });
+
+ return result;
+ }
+
PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
dot(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
{
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index 7899ea66d6ca7bb0f1fcf79d38b524cc07fed1e5..0168633cea1d2d2e764b6b835093d1f17504024d 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -396,6 +396,44 @@ let d:R, d = det(A);
-2.6, 5.2, -3.5}));
}
+ { // trace
+ tested_function_set.insert("trace:R^1x1");
+ std::string_view data = R"(
+import math;
+let A:R^1x1, A = -2;
+let t:R, t = trace(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "t", trace(TinyMatrix<1>{-2}));
+ }
+
+ { // trace
+ tested_function_set.insert("trace:R^2x2");
+ std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+ [ 3.2, 2.3]];
+let t:R, t = trace(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "t",
+ trace(TinyMatrix<2>{-2.5, 3.2, //
+ +3.2, 2.3}));
+ }
+
+ { // trace
+ tested_function_set.insert("trace:R^3x3");
+ std::string_view data = R"(
+import math;
+let A:R^3x3, A = [[-2.5, 2.9,-1.3],
+ [ 3.2, 2.3, 2.7],
+ [-2.6, 5.2,-3.5]];
+let t:R, t = trace(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "t",
+ trace(TinyMatrix<3>{-2.5, 2.9, -1.3, //
+ +3.2, 2.3, +2.7, //
+ -2.6, 5.2, -3.5}));
+ }
+
{ // inverse
tested_function_set.insert("inverse:R^1x1");
std::string_view data = R"(
@@ -434,6 +472,44 @@ let invA:R^3x3, invA = inverse(A);
-2.6, 5.2, -3.5}));
}
+ { // transpose
+ tested_function_set.insert("transpose:R^1x1");
+ std::string_view data = R"(
+import math;
+let A:R^1x1, A = -2;
+let tA:R^1x1, tA = transpose(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "tA", transpose(TinyMatrix<1>{-2}));
+ }
+
+ { // transpose
+ tested_function_set.insert("transpose:R^2x2");
+ std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+ [ 3.2, 2.3]];
+let tA:R^2x2, tA = transpose(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "tA",
+ transpose(TinyMatrix<2>{-2.5, 3.2, //
+ +3.2, 2.3}));
+ }
+
+ { // transpose
+ tested_function_set.insert("transpose:R^3x3");
+ std::string_view data = R"(
+import math;
+let A:R^3x3, A = [[-2.5, 2.9,-1.3],
+ [ 3.2, 2.3, 2.7],
+ [-2.6, 5.2,-3.5]];
+let tA:R^3x3, tA = transpose(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "tA",
+ transpose(TinyMatrix<3>{-2.5, 2.9, -1.3, //
+ +3.2, 2.3, +2.7, //
+ -2.6, 5.2, -3.5}));
+ }
+
MathModule math_module;
bool missing_test = false;
diff --git a/tests/test_DiscreteFunctionP0.cpp b/tests/test_DiscreteFunctionP0.cpp
index 6c420b8443b11b40d7931a7d628e6374d3a83c49..48c107e4ad2f07391a49270a04c9a3dc20d2db4a 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -3450,6 +3450,28 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
}
}
+ SECTION("trace(Ah)")
+ {
+ DiscreteFunctionP0<Dimension, TinyMatrix<2>> Ah{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ Ah[cell_id] = TinyMatrix<2>{x + 1, 2 * x - 3, //
+ -0.2 * x - 1, 2 + x};
+ });
+
+ {
+ DiscreteFunctionP0 result = trace(Ah);
+ bool is_same = true;
+ parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+ if (result[cell_id] != trace(Ah[cell_id])) {
+ is_same = false;
+ }
+ });
+ REQUIRE(is_same);
+ }
+ }
+
SECTION("inverse(Ah)")
{
DiscreteFunctionP0<Dimension, TinyMatrix<2>> Ah{mesh};
@@ -3472,6 +3494,28 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
}
}
+ SECTION("transpose(Ah)")
+ {
+ DiscreteFunctionP0<Dimension, TinyMatrix<2>> Ah{mesh};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double x = xj[cell_id][0];
+ Ah[cell_id] = TinyMatrix<2>{x + 1, 2 * x - 3, //
+ -0.2 * x - 1, 2 + x};
+ });
+
+ {
+ DiscreteFunctionP0 result = transpose(Ah);
+ bool is_same = true;
+ parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+ if (result[cell_id] != transpose(Ah[cell_id])) {
+ is_same = false;
+ }
+ });
+ REQUIRE(is_same);
+ }
+ }
+
SECTION("scalar sum")
{
const CellValue<const double> cell_value = positive_function.cellValues();
diff --git a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
index 0fcc3c148a5da83d1f43663bb52c1eda05e02fd9..6d00782f96379548ebb66927044782def35612e8 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -603,6 +603,71 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
}
+ SECTION("trace Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = trace(p_R1x1_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R1x1_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = trace(p_R2x2_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R2x2_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = trace(p_R3x3_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R3x3_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(trace(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(trace(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(trace(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(trace(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("inverse Vh -> Vh")
{
{
@@ -686,6 +751,89 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
}
+ SECTION("transpose Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = transpose(p_R1x1_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+ auto values = p_R1x1_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id](0, 0) != transpose(values[cell_id])(0, 0)) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = transpose(p_R2x2_u);
+ constexpr size_t nb_row = 2;
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+ auto values = p_R2x2_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ const TinyMatrix invA = transpose(values[cell_id]);
+ const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+ for (size_t i = 0; i < nb_row; ++i) {
+ for (size_t j = 0; j < nb_row; ++j) {
+ if (fu_i(i, j) != invA(i, j)) {
+ is_same = false;
+ }
+ }
+ }
+ if (not is_same) {
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = transpose(p_R3x3_u);
+ constexpr size_t nb_row = 3;
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+ auto values = p_R3x3_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ const TinyMatrix invA = transpose(values[cell_id]);
+ const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+ for (size_t i = 0; i < nb_row; ++i) {
+ for (size_t j = 0; j < nb_row; ++j) {
+ if (fu_i(i, j) != invA(i, j)) {
+ is_same = false;
+ }
+ }
+ }
+ if (not is_same) {
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(transpose(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(transpose(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(transpose(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(transpose(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("sum_of_Vh Vh -> Vh")
{
{
@@ -1323,6 +1471,71 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
}
+ SECTION("trace Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = trace(p_R1x1_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R1x1_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = trace(p_R2x2_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R2x2_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = trace(p_R3x3_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R3x3_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(trace(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(trace(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(trace(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(trace(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("inverse Vh -> Vh")
{
{
@@ -1406,6 +1619,89 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
}
+ SECTION("transpose Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = transpose(p_R1x1_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+ auto values = p_R1x1_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id](0, 0) != transpose(values[cell_id])(0, 0)) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = transpose(p_R2x2_u);
+ constexpr size_t nb_row = 2;
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+ auto values = p_R2x2_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ const TinyMatrix invA = transpose(values[cell_id]);
+ const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+ for (size_t i = 0; i < nb_row; ++i) {
+ for (size_t j = 0; j < nb_row; ++j) {
+ if (fu_i(i, j) != invA(i, j)) {
+ is_same = false;
+ }
+ }
+ }
+ if (not is_same) {
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = transpose(p_R3x3_u);
+ constexpr size_t nb_row = 3;
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+ auto values = p_R3x3_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ const TinyMatrix invA = transpose(values[cell_id]);
+ const TinyMatrix<nb_row>& fu_i = fu[cell_id];
+ for (size_t i = 0; i < nb_row; ++i) {
+ for (size_t j = 0; j < nb_row; ++j) {
+ if (fu_i(i, j) != invA(i, j)) {
+ is_same = false;
+ }
+ }
+ }
+ if (not is_same) {
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(transpose(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(transpose(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(transpose(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(transpose(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("sum_of_Vh Vh -> Vh")
{
{
@@ -2041,6 +2337,71 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
}
+ SECTION("trace Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = trace(p_R1x1_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R1x1_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = trace(p_R2x2_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R2x2_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = trace(p_R3x3_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_fu);
+
+ auto values = p_R3x3_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != trace(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(trace(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(trace(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(trace(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(trace(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("inverse Vh -> Vh")
{
{
@@ -2108,6 +2469,73 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
}
+ SECTION("transpose Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = transpose(p_R1x1_u);
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*p_fu);
+
+ auto values = p_R1x1_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != transpose(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = transpose(p_R2x2_u);
+ constexpr size_t nb_row = 2;
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+ auto values = p_R2x2_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != transpose(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = transpose(p_R3x3_u);
+ constexpr size_t nb_row = 3;
+
+ REQUIRE(p_fu.use_count() > 0);
+ REQUIRE_NOTHROW(dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu));
+
+ const auto& fu = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<nb_row>>&>(*p_fu);
+
+ auto values = p_R3x3_u->cellValues();
+ bool is_same = true;
+ for (CellId cell_id = 0; cell_id < values.numberOfItems(); ++cell_id) {
+ if (fu[cell_id] != transpose(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(transpose(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(transpose(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(transpose(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(transpose(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("sum_of_Vh Vh -> Vh")
{
{
diff --git a/tests/test_MathModule.cpp b/tests/test_MathModule.cpp
index cea6c79519fe3530f64e3adc0cb1fd5fc6c57274..2cb3481b1e6cde7561733d8ea03decd46153db9d 100644
--- a/tests/test_MathModule.cpp
+++ b/tests/test_MathModule.cpp
@@ -13,7 +13,7 @@ TEST_CASE("MathModule", "[language]")
MathModule math_module;
const auto& name_builtin_function = math_module.getNameBuiltinFunctionMap();
- REQUIRE(name_builtin_function.size() == 36);
+ REQUIRE(name_builtin_function.size() == 42);
SECTION("Z -> N")
{
@@ -509,6 +509,56 @@ TEST_CASE("MathModule", "[language]")
const double result = det(arg);
REQUIRE(std::get<double>(result_variant) == result);
}
+
+ SECTION("trace:R^1x1")
+ {
+ TinyMatrix<1> arg{3};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("trace:R^1x1");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg_variant});
+
+ const double result = trace(arg);
+ REQUIRE(std::get<double>(result_variant) == result);
+ }
+
+ SECTION("trace:R^2x2")
+ {
+ TinyMatrix<2> arg{3, 2, 1, -2};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("trace:R^2x2");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg_variant});
+
+ const double result = trace(arg);
+ REQUIRE(std::get<double>(result_variant) == result);
+ }
+
+ SECTION("trace:R^3x3")
+ {
+ TinyMatrix<3> arg{3, 2, 4, //
+ -1, 4, 1, //
+ -4, -1, 5};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("trace:R^3x3");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg_variant});
+
+ const double result = trace(arg);
+ REQUIRE(std::get<double>(result_variant) == result);
+ }
}
SECTION("R^dxd -> R^dxd")
@@ -562,5 +612,55 @@ TEST_CASE("MathModule", "[language]")
const TinyMatrix<3> result = inverse(arg);
REQUIRE(std::get<TinyMatrix<3>>(result_variant) == result);
}
+
+ SECTION("transpose:R^1x1")
+ {
+ TinyMatrix<1> arg{3};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("transpose:R^1x1");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg_variant});
+
+ const TinyMatrix<1> result = transpose(arg);
+ REQUIRE(std::get<TinyMatrix<1>>(result_variant) == result);
+ }
+
+ SECTION("transpose:R^2x2")
+ {
+ TinyMatrix<2> arg{3, 2, 1, -2};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("transpose:R^2x2");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg_variant});
+
+ const TinyMatrix<2> result = transpose(arg);
+ REQUIRE(std::get<TinyMatrix<2>>(result_variant) == result);
+ }
+
+ SECTION("transpose:R^3x3")
+ {
+ TinyMatrix<3> arg{3, 2, 4, //
+ -1, 4, 1, //
+ -4, -1, 5};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("transpose:R^3x3");
+ REQUIRE(i_function != name_builtin_function.end());
+
+ IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
+ DataVariant result_variant = function_embedder.apply({arg_variant});
+
+ const TinyMatrix<3> result = transpose(arg);
+ REQUIRE(std::get<TinyMatrix<3>>(result_variant) == result);
+ }
}
}
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 2f4c09d4ad9ce1dbb5e3d856d8737d27344459f2..f6f09cf55dd2ad1ffef55c9275b989c203598f7e 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -194,6 +194,17 @@ TEST_CASE("TinyMatrix", "[algebra]")
REQUIRE(det(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
}
+ SECTION("checking for trace calculations")
+ {
+ REQUIRE(trace(TinyMatrix<1, 1, int>(6)) == 6);
+ REQUIRE(trace(TinyMatrix<2, 2, int>(5, 1, -3, 6)) == 5 + 6);
+ REQUIRE(trace(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1 + 2 + 3);
+ REQUIRE(trace(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == 6 + 34 + 7);
+ REQUIRE(trace(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
+ Catch::Approx(1 + 4 + 2 + 17.5).epsilon(1E-14));
+ REQUIRE(trace(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 1 + 0 + 1 + 2);
+ }
+
SECTION("checking for inverse calculations")
{
{