From 3e51ffef04a9468da8906e6f308f345a6b127abb Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Wed, 1 Feb 2023 01:08:56 +0100
Subject: [PATCH] Add det and inverse functions to the language
- These apply to square matrices and discrete functions of square
matrices.
- Tests have been added
- Documentation is also updated
---
doc/userdoc.org | 49 +-
.../modules/MathFunctionRegisterForVh.cpp | 14 +
src/language/modules/MathModule.cpp | 15 +
...EmbeddedIDiscreteFunctionMathFunctions.cpp | 102 +++++
...EmbeddedIDiscreteFunctionMathFunctions.hpp | 4 +
src/scheme/DiscreteFunctionP0.hpp | 26 ++
tests/test_BuiltinFunctionProcessor.cpp | 76 ++++
tests/test_DiscreteFunctionP0.cpp | 44 ++
...EmbeddedIDiscreteFunctionMathFunctions.cpp | 428 ++++++++++++++++++
tests/test_MathModule.cpp | 108 ++++-
10 files changed, 861 insertions(+), 5 deletions(-)
diff --git a/doc/userdoc.org b/doc/userdoc.org
index a4d0ce4a9..320a9c183 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -2588,6 +2588,33 @@ $\mathbb{R}^1$, $\mathbb{R}^2$ and $\mathbb{R}^3$.
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}$.
+#+NAME: det-examples
+#+BEGIN_SRC pugs :exports both :results output
+ import math;
+ cout << "det([[1.2]]) = " << det([[1.2]]) << "\n";
+ cout << "det([[1,2],[3,4]]) = " << det([[1,2],[3,4]]) << "\n";
+ cout << "det([[1,2,3],[4,5,6],[7,8,9]]) = "
+ << det([[1,2,3],[4,5,6],[7,8,9]]) << "\n";
+#+END_SRC
+The output is
+#+RESULTS: det-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";
+#+END_SRC
+The output is
+#+RESULTS: inverse-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
@@ -3196,11 +3223,12 @@ description.
****** ~Vh -> Vh~
-These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
-return value is also a $\mathbb{P}_0(\mathbb{R})$ function. The value
-is simply the application of the function to the cell values.
+The majority of these functions are defined for
+$\mathbb{P}_0(\mathbb{R})$ data and the return value is also a
+$\mathbb{P}_0(\mathbb{R})$ function. The value is simply the
+application of the function to the cell values.
-Here is the list of the functions
+Here is the list of these functions
- ~abs: Vh -> Vh~
- ~acos: Vh -> Vh~
- ~acosh: Vh -> Vh~
@@ -3218,6 +3246,19 @@ Here is the list of the functions
- ~tan: Vh -> Vh~
- ~tanh: Vh -> Vh~
+The ~det~ 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~
+
+The ~inverse~ 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~
+
****** ~Vh*Vh -> Vh~
These functions are defined for $\mathbb{P}_0(\mathbb{R})$ data and the
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index fcf3b31a7..d0468777f 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -209,6 +209,20 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
));
+ scheme_module._addBuiltinFunction("det", std::function(
+
+ [](std::shared_ptr<const IDiscreteFunction> A)
+ -> std::shared_ptr<const IDiscreteFunction> { return det(A); }
+
+ ));
+
+ scheme_module._addBuiltinFunction("inverse", std::function(
+
+ [](std::shared_ptr<const IDiscreteFunction> A)
+ -> std::shared_ptr<const IDiscreteFunction> { return inverse(A); }
+
+ ));
+
scheme_module._addBuiltinFunction("min", std::function(
[](std::shared_ptr<const IDiscreteFunction> a) -> double { return min(a); }
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index 4e8e10d15..5f2de33d8 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -70,6 +70,21 @@ MathModule::MathModule()
this->_addBuiltinFunction("dot", std::function([](const TinyVector<3>& x, const TinyVector<3>& y) -> double {
return dot(x, y);
}));
+
+ this->_addBuiltinFunction("det", std::function([](const TinyMatrix<1> A) -> double { return det(A); }));
+
+ this->_addBuiltinFunction("det", std::function([](const TinyMatrix<2>& A) -> double { return det(A); }));
+
+ this->_addBuiltinFunction("det", std::function([](const TinyMatrix<3>& A) -> double { return det(A); }));
+
+ this->_addBuiltinFunction("inverse",
+ std::function([](const TinyMatrix<1> A) -> TinyMatrix<1> { return inverse(A); }));
+
+ this->_addBuiltinFunction("inverse",
+ std::function([](const TinyMatrix<2>& A) -> TinyMatrix<2> { return inverse(A); }));
+
+ this->_addBuiltinFunction("inverse",
+ std::function([](const TinyMatrix<3>& A) -> TinyMatrix<3> { return inverse(A); }));
}
void
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
index f83bd4185..e4ac44a33 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -503,6 +503,108 @@ template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<2>&,
template std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<3>&,
const std::shared_ptr<const IDiscreteFunction>&);
+template <typename MatrixType>
+std::shared_ptr<const IDiscreteFunction>
+det(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>>(det(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 2: {
+ using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+ return std::make_shared<const DiscreteFunctionP0<2, double>>(det(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 3: {
+ using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+ return std::make_shared<const DiscreteFunctionP0<3, double>>(det(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+det(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 det<TinyMatrix<1>>(A);
+ }
+ case 2: {
+ return det<TinyMatrix<2>>(A);
+ }
+ case 3: {
+ return det<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)
+{
+ switch (A->mesh()->dimension()) {
+ case 1: {
+ using DiscreteFunctionType = DiscreteFunctionP0<1, MatrixType>;
+ return std::make_shared<const DiscreteFunctionType>(inverse(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 2: {
+ using DiscreteFunctionType = DiscreteFunctionP0<2, MatrixType>;
+ return std::make_shared<const DiscreteFunctionType>(inverse(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ case 3: {
+ using DiscreteFunctionType = DiscreteFunctionP0<3, MatrixType>;
+ return std::make_shared<const DiscreteFunctionType>(inverse(dynamic_cast<const DiscreteFunctionType&>(*A)));
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+inverse(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 inverse<TinyMatrix<1>>(A);
+ }
+ case 2: {
+ return inverse<TinyMatrix<2>>(A);
+ }
+ case 3: {
+ return inverse<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 2369429a6..1908dd30c 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionMathFunctions.hpp
@@ -65,6 +65,10 @@ template <size_t VectorDimension>
std::shared_ptr<const IDiscreteFunction> dot(const TinyVector<VectorDimension>&,
const std::shared_ptr<const IDiscreteFunction>&);
+std::shared_ptr<const IDiscreteFunction> det(const std::shared_ptr<const IDiscreteFunction>&);
+
+std::shared_ptr<const IDiscreteFunction> inverse(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 bb6b484a6..6e1ef6a0d 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -524,11 +524,37 @@ class DiscreteFunctionP0 : public IDiscreteFunction
return result;
}
+ PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
+ det(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] = det(A[cell_id]); });
+
+ return result;
+ }
+
+ PUGS_INLINE friend DiscreteFunctionP0
+ inverse(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] = inverse(A[cell_id]); });
+
+ return result;
+ }
+
PUGS_INLINE friend DiscreteFunctionP0<Dimension, double>
dot(const DiscreteFunctionP0& f, const DiscreteFunctionP0& g)
{
static_assert(is_tiny_vector_v<DataType>);
Assert(f.m_cell_values.isBuilt() and g.m_cell_values.isBuilt());
+ Assert(f.m_mesh == g.m_mesh);
DiscreteFunctionP0<Dimension, double> result{f.m_mesh};
parallel_for(
f.m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { result[cell_id] = dot(f[cell_id], g[cell_id]); });
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index 60a4bd491..7899ea66d 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -358,6 +358,82 @@ let s:R, s = dot(x,y);
CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", dot(TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
}
+ { // det
+ tested_function_set.insert("det:R^1x1");
+ std::string_view data = R"(
+import math;
+let A:R^1x1, A = -2;
+let d:R, d = det(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "d", det(TinyMatrix<1>{-2}));
+ }
+
+ { // det
+ tested_function_set.insert("det:R^2x2");
+ std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+ [ 3.2, 2.3]];
+let d:R, d = det(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "d",
+ det(TinyMatrix<2>{-2.5, 3.2, //
+ +3.2, 2.3}));
+ }
+
+ { // det
+ tested_function_set.insert("det: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 d:R, d = det(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "d",
+ det(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"(
+import math;
+let A:R^1x1, A = -2;
+let invA:R^1x1, invA = inverse(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "invA", inverse(TinyMatrix<1>{-2}));
+ }
+
+ { // inverse
+ tested_function_set.insert("inverse:R^2x2");
+ std::string_view data = R"(
+import math;
+let A:R^2x2, A = [[-2.5, 3.2],
+ [ 3.2, 2.3]];
+let invA:R^2x2, invA = inverse(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "invA",
+ inverse(TinyMatrix<2>{-2.5, 3.2, //
+ +3.2, 2.3}));
+ }
+
+ { // inverse
+ tested_function_set.insert("inverse: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 invA:R^3x3, invA = inverse(A);
+)";
+ CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "invA",
+ inverse(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 9d83de536..6c420b844 100644
--- a/tests/test_DiscreteFunctionP0.cpp
+++ b/tests/test_DiscreteFunctionP0.cpp
@@ -3428,6 +3428,50 @@ TEST_CASE("DiscreteFunctionP0", "[scheme]")
CHECK_STD_BINARY_MATH_FUNCTION_WITH_LHS_VALUE(u, vh, dot);
}
+ SECTION("det(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 = det(Ah);
+ bool is_same = true;
+ parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+ if (result[cell_id] != det(Ah[cell_id])) {
+ is_same = false;
+ }
+ });
+ REQUIRE(is_same);
+ }
+ }
+
+ SECTION("inverse(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 = inverse(Ah);
+ bool is_same = true;
+ parallel_for(Ah.cellValues().numberOfItems(), [&](const CellId cell_id) {
+ if (result[cell_id] != inverse(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 9508afc69..0fcc3c148 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -538,6 +538,154 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("det Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(det(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(det(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(det(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
+ SECTION("inverse Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = inverse(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) != inverse(values[cell_id])(0, 0)) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = inverse(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 = inverse(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 = inverse(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 = inverse(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(inverse(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(inverse(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(inverse(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("sum_of_Vh Vh -> Vh")
{
{
@@ -1110,6 +1258,154 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("det Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(det(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(det(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(det(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
+ SECTION("inverse Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = inverse(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) != inverse(values[cell_id])(0, 0)) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = inverse(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 = inverse(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 = inverse(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 = inverse(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(inverse(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(inverse(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(inverse(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(inverse(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
SECTION("sum_of_Vh Vh -> Vh")
{
{
@@ -1680,6 +1976,138 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
REQUIRE_THROWS_WITH(dot(p_Vector3_u, p_Vector2_w), "error: operands have different dimension");
}
+ SECTION("det Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = det(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] != det(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(det(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(det(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(det(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(det(p_R3_u), "error: invalid operand type Vh(P0:R^3)");
+ }
+
+ SECTION("inverse Vh -> Vh")
+ {
+ {
+ std::shared_ptr p_fu = inverse(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] != inverse(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = inverse(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] != inverse(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ {
+ std::shared_ptr p_fu = inverse(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] != inverse(values[cell_id])) {
+ is_same = false;
+ break;
+ }
+ }
+ REQUIRE(is_same);
+ }
+
+ REQUIRE_THROWS_WITH(inverse(p_u), "error: invalid operand type Vh(P0:R)");
+ REQUIRE_THROWS_WITH(inverse(p_R1_u), "error: invalid operand type Vh(P0:R^1)");
+ REQUIRE_THROWS_WITH(inverse(p_R2_u), "error: invalid operand type Vh(P0:R^2)");
+ REQUIRE_THROWS_WITH(inverse(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 99966d4ba..cea6c7951 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() == 30);
+ REQUIRE(name_builtin_function.size() == 36);
SECTION("Z -> N")
{
@@ -457,4 +457,110 @@ TEST_CASE("MathModule", "[language]")
REQUIRE(std::get<double>(result_variant) == result);
}
}
+
+ SECTION("R^dxd -> double")
+ {
+ SECTION("det:R^1x1")
+ {
+ TinyMatrix<1> arg{3};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("det: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 = det(arg);
+ REQUIRE(std::get<double>(result_variant) == result);
+ }
+
+ SECTION("det:R^2x2")
+ {
+ TinyMatrix<2> arg{3, 2, 1, -2};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("det: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 = det(arg);
+ REQUIRE(std::get<double>(result_variant) == result);
+ }
+
+ SECTION("det: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("det: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 = det(arg);
+ REQUIRE(std::get<double>(result_variant) == result);
+ }
+ }
+
+ SECTION("R^dxd -> R^dxd")
+ {
+ SECTION("inverse:R^1x1")
+ {
+ TinyMatrix<1> arg{3};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("inverse: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 = inverse(arg);
+ REQUIRE(std::get<TinyMatrix<1>>(result_variant) == result);
+ }
+
+ SECTION("inverse:R^2x2")
+ {
+ TinyMatrix<2> arg{3, 2, 1, -2};
+
+ DataVariant arg_variant = arg;
+
+ auto i_function = name_builtin_function.find("inverse: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 = inverse(arg);
+ REQUIRE(std::get<TinyMatrix<2>>(result_variant) == result);
+ }
+
+ SECTION("inverse: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("inverse: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 = inverse(arg);
+ REQUIRE(std::get<TinyMatrix<3>>(result_variant) == result);
+ }
+ }
}
--
GitLab