diff --git a/src/amdis/common/FieldMatVec.hpp b/src/amdis/common/FieldMatVec.hpp index 6a316262b4e16a8e517d4989c385dd73190378fc..a7d61c19836310ef5d7b737a9215b1ce545b81ad 100644 --- a/src/amdis/common/FieldMatVec.hpp +++ b/src/amdis/common/FieldMatVec.hpp @@ -23,22 +23,168 @@ namespace std namespace Dune { - // some arithmetic operations with FieldVector + namespace MatVec + { + /// Traits to detect fixed size containers like FieldVector and FieldMatrix + /// @{ + template + struct IsMatrix : std::false_type {}; + + template + struct IsMatrix> : std::true_type {}; + + template + struct IsMatrix> : std::true_type {}; + + + template + struct IsVector : std::false_type {}; + + template + struct IsVector> : std::true_type {}; + + + template + struct IsMatVec + : std::integral_constant::value || IsVector::value> {}; + /// @} + + /// Convert the field_Type of container to type S + /// @{ + template + struct MakeMatVec + { + using type = A; + }; + + template + struct MakeMatVec,S> + { + using type = FieldMatrix; + }; + + template + struct MakeMatVec,S> + { + using type = DiagonalMatrix; + }; + + template + struct MakeMatVec,S> + { + using type = FieldVector; + }; + /// @} + + /// Convert pseudo-scalar to real scalar type + /// @{ + template + decltype(auto) simplify(T&& t) + { + return std::forward(t); + } + + template + T simplify(FieldVector const& t) + { + return t[0]; + } + + template + T simplify(FieldMatrix const& t) + { + return t[0][0]; + } + + template + T simplify(DiagonalMatrix const& t) + { + return t.diagonal(0); + } + /// @} + + // returns -a + template + auto negate(A const& a); + + // returns a+b + template + auto plus(A a, B const& b) + { + return a += b; + } + + // returns a-b + template + auto minus(A a, B const& b) + { + return a -= b; + } + + // returns a*b + template ::value || IsNumber::value, int> = 0> + auto multiplies(A const& a, B const& b); + + template + auto multiplies(FieldVector const& a, FieldVector const& b); + + template ::value && IsVector::value, int> = 0> + auto multiplies(Mat const& mat, Vec const& vec); + + template ::value && IsMatrix::value, int> = 0> + auto multiplies(Vec const& vec, Mat const& mat); + + template + auto multiplies(FieldMatrix const& a, FieldMatrix const& b); + + // return a/b + template + auto divides(A a, B const& b) + { + return a /= b; + } + + } // end namespace MatVec + + // some arithmetic operations with FieldVector and FieldMatrix + + template ::value, int> = 0> + auto operator-(A const& a) + { + return MatVec::negate(MatVec::simplify(a)); + } - template - FieldVector operator-(FieldVector v); + template ::value || MatVec::IsMatVec::value, int> = 0> + auto operator+(A const& a, B const& b) + { + return MatVec::plus(MatVec::simplify(a), MatVec::simplify(b)); + } - template ::value, int> = 0> - FieldVector operator*(FieldVector v, S factor); + template ::value || MatVec::IsMatVec::value, int> = 0> + auto operator-(A const& a, B const& b) + { + return MatVec::minus(MatVec::simplify(a), MatVec::simplify(b)); + } - template ::value, int> = 0> - FieldVector operator*(S factor, FieldVector v); + template ::value || MatVec::IsMatVec::value, int> = 0> + auto operator*(A const& a, B const& b) + { + return MatVec::multiplies(MatVec::simplify(a), MatVec::simplify(b)); + } - template ::value, int> = 0> - FieldVector operator/(FieldVector v, S factor); + template ::value || MatVec::IsMatVec::value, int> = 0> + auto operator/(A const& a, B const& b) + { + return MatVec::divides(MatVec::simplify(a), MatVec::simplify(b)); + } // ---------------------------------------------------------------------------- @@ -197,38 +343,8 @@ namespace Dune template FieldMatrix trans(FieldMatrix const& A); - template - FieldMatrix operator-(FieldMatrix A); - - - template ::value, int> = 0> - FieldMatrix operator*(S scalar, FieldMatrix A); - - template ::value, int> = 0> - FieldMatrix operator*(FieldMatrix A, S scalar); - - template ::value, int> = 0> - FieldMatrix operator/(FieldMatrix A, S scalar); - - - template - FieldMatrix operator+(FieldMatrix A, FieldMatrix const& B); - - template - FieldMatrix operator-(FieldMatrix A, FieldMatrix const& B); - - - template - FieldVector operator*(FieldMatrix const& mat, FieldVector const& vec); - // ----------------------------------------------------------------------------- - template - FieldMatrix multiplies(FieldMatrix const& A, FieldMatrix const& B); - template FieldMatrix multiplies_AtB(FieldMatrix const& A, FieldMatrix const& B); @@ -243,20 +359,6 @@ namespace Dune // ----------------------------------------------------------------------------- - template - T operator*(FieldVector lhs, FieldVector rhs); - - template - T operator*(FieldMatrix lhs, FieldMatrix rhs); - - template - T operator*(FieldVector lhs, FieldMatrix rhs); - - template - T operator*(FieldMatrix lhs, FieldVector rhs); - - // ----------------------------------------------------------------------------- - template T const& at(FieldMatrix const& vec, std::size_t i); diff --git a/src/amdis/common/FieldMatVec.inc.hpp b/src/amdis/common/FieldMatVec.inc.hpp index cadb8d3a551bc7a2d6f61d281880585f7b64f964..8e57a3ea506ba0386e70145be9426d01d8ca3604 100644 --- a/src/amdis/common/FieldMatVec.inc.hpp +++ b/src/amdis/common/FieldMatVec.inc.hpp @@ -11,34 +11,81 @@ namespace Dune { -// some arithmetic operations with FieldVector +namespace MatVec { -template -FieldVector operator-(FieldVector v) -{ - return v *= -1; -} + template + auto negate(A const& a) + { + return multiplies(a, -1); + } -template ::value, int>> -FieldVector operator*(FieldVector v, S factor) -{ - return v *= factor; -} + template ::value || IsNumber::value, int>> + auto multiplies(A const& a, B const& b) + { + using T = std::common_type_t::field_type, typename FieldTraits::field_type>; +#if AMDIS_HAS_CXX_CONSTEXPR_IF + if constexpr(IsNumber::value) { + typename MakeMatVec::type b_{b}; + return b_ *= a; + } else { + typename MakeMatVec::type a_{a}; + return a_ *= b; + } +#else + return Hybrid::ifElse(IsNumber{}, + [&](auto id) { typename MakeMatVec::type b_{b}; return id(b_) *= id(a); }, + [&](auto id) { typename MakeMatVec::type a_{a}; return id(a_) *= id(b); }); +#endif + } -template ::value, int>> -FieldVector operator*(S factor, FieldVector v) -{ - return v *= factor; -} + template + auto multiplies(FieldVector const& a, FieldVector const& b) + { + return a.dot(b); + } -template ::value, int>> -FieldVector operator/(FieldVector v, S factor) -{ - return v /= factor; -} + + template ::value && IsVector::value, int>> + auto multiplies(Mat const& mat, Vec const& vec) + { + static_assert(Mat::cols == Vec::dimension, ""); + using T = std::common_type_t::field_type, typename FieldTraits::field_type>; + FieldVector y; + mat.mv(vec, y); + return y; + } + + + template ::value && IsMatrix::value, int>> + auto multiplies(Vec const& vec, Mat const& mat) + { + static_assert(Mat::rows == Vec::dimension, ""); + using T = std::common_type_t::field_type, typename FieldTraits::field_type>; + FieldVector y; + mat.mtv(vec, y); + return y; + } + + + template + auto multiplies(FieldMatrix const& a, FieldMatrix const& b) + { + FieldMatrix,L,N> C; + + for (int i = 0; i < L; ++i) { + for (int j = 0; j < N; ++j) { + C[i][j] = 0; + for (int k = 0; k < M; ++k) + C[i][j] += a[i][k]*b[k][j]; + } + } + return C; + } + +} // end namespace MatVec // ---------------------------------------------------------------------------- @@ -349,57 +396,6 @@ FieldMatrix trans(FieldMatrix const& A) return At; } -template -FieldMatrix operator-(FieldMatrix A) -{ - return A *= -1; -} - -template ::value, int>> -FieldMatrix operator*(S scalar, FieldMatrix A) -{ - return A *= scalar; -} - -template ::value, int>> -FieldMatrix operator*(FieldMatrix A, S scalar) -{ - return A *= scalar; -} - -template ::value, int>> -FieldMatrix operator/(FieldMatrix A, S scalar) -{ - return A /= scalar; -} - -template -FieldMatrix operator+(FieldMatrix A, FieldMatrix const& B) -{ - return A += B; -} - -template -FieldMatrix operator-(FieldMatrix A, FieldMatrix const& B) -{ - return A -= B; -} - - -template -FieldVector operator*(FieldMatrix const& mat, FieldVector const& vec) -{ - return Dune::FMatrixHelp::mult(mat, vec); -} - -template -FieldMatrix multiplies(FieldMatrix const& A, FieldMatrix const& B) -{ - return A.rightmultiplyany(B); -} template FieldMatrix multiplies_AtB(FieldMatrix const& A, FieldMatrix const& B) @@ -447,32 +443,6 @@ FieldMatrix& multiplies_ABt(FieldMatrix const& A, DiagonalMatri return C; } - -template -T operator*(FieldVector lhs, FieldVector rhs) -{ - return lhs[0] * rhs[0]; -} - -template -T operator*(FieldMatrix lhs, FieldMatrix rhs) -{ - return lhs[0][0] * rhs[0][0]; -} - -template -T operator*(FieldVector lhs, FieldMatrix rhs) -{ - return lhs[0] * rhs[0][0]; -} - -template -T operator*(FieldMatrix lhs, FieldVector rhs) -{ - return lhs[0][0] * rhs[0]; -} - - template T const& at(FieldMatrix const& vec, std::size_t i) { diff --git a/test/FieldMatVecTest.cpp b/test/FieldMatVecTest.cpp index 9f3768389b75c2da00afeec4c0e4b4d914d6cf57..53ea56a4695e1e0f73ed1f46ae4f459fe49efe0e 100644 --- a/test/FieldMatVecTest.cpp +++ b/test/FieldMatVecTest.cpp @@ -101,7 +101,7 @@ void test2() AMDIS_TEST_EQ( det(A), 1.0 ); AMDIS_TEST_EQ( det(B), 0.0 ); - AMDIS_TEST_EQ( det(multiplies(C,D)), det(C)*det(D) ); + AMDIS_TEST_EQ( det(C*D), det(C)*det(D) ); AMDIS_TEST_EQ( det(trans(C)), det(C) ); AMDIS_TEST_EQ( det(2.0*D), Math::pow<3>(2.0)*det(D) ); @@ -201,6 +201,15 @@ void test3() AMDIS_TEST_EQ( vi1*mi1, 1 ); AMDIS_TEST_EQ( mi1*mi1, 1 ); AMDIS_TEST_EQ( mi1*vi1, 1 ); + + AMDIS_TEST_EQ( vd1*vi1, 1.0 ); + AMDIS_TEST_EQ( vd1*mi1, 1.0 ); + AMDIS_TEST_EQ( vi1*vd1, 1.0 ); + AMDIS_TEST_EQ( vi1*md1, 1.0 ); + AMDIS_TEST_EQ( md1*vi1, 1.0 ); + AMDIS_TEST_EQ( md1*mi1, 1.0 ); + AMDIS_TEST_EQ( mi1*vd1, 1.0 ); + AMDIS_TEST_EQ( mi1*md1, 1.0 ); }