diff --git a/doc/userdoc.org b/doc/userdoc.org index a4d0ce4a9ac4466ce3fd3c0467a2d3d3ade63965..320a9c183270985d0d3ad17fbda9ebd562cad2c9 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 fcf3b31a73da2c47bb59e27bc169195255effb30..d0468777f4528873115939a69027dcb6844d936f 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 4e8e10d1534f55f8c02cccb2f3f48974fdc944d2..5f2de33d8502fa84edd809f2dfbfa51c91c12ef4 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 f83bd418592c29b5cded8b77bda6c717ff3af804..e4ac44a3328d40415a4cbdca52600e224f6885cc 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 2369429a6267a8fdbc803d75b4213703dd970a48..1908dd30cf16c6cc7fa3c3281d42429b0c661580 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 bb6b484a6bc10bd803e16b8bb54a7cba25a349de..6e1ef6a0d28a3436543fa090bee813586874bbee 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 60a4bd491f392b7cc0c303d94531edcf8634edbc..7899ea66d6ca7bb0f1fcf79d38b524cc07fed1e5 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 9d83de5368f81f1dc46570234c2a7ed19ebc2da4..6c420b8443b11b40d7931a7d628e6374d3a83c49 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 9508afc6987a8bf67a2a090892ba1bee211cd176..0fcc3c148a5da83d1f43663bb52c1eda05e02fd9 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 99966d4bacebb938fce5bc686661b5c53e8d708e..cea6c79519fe3530f64e3adc0cb1fd5fc6c57274 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); + } + } }