Commit 66f67313 authored by Praetorius, Simon's avatar Praetorius, Simon

Redesign of FieldMatrix + FieldVector arithmetic operations

parent 736ed6d7
...@@ -23,22 +23,168 @@ namespace std ...@@ -23,22 +23,168 @@ namespace std
namespace Dune namespace Dune
{ {
// some arithmetic operations with FieldVector namespace MatVec
{
/// Traits to detect fixed size containers like FieldVector and FieldMatrix
/// @{
template <class T>
struct IsMatrix : std::false_type {};
template <class T, int M, int N>
struct IsMatrix<FieldMatrix<T,M,N>> : std::true_type {};
template <class T, int N>
struct IsMatrix<DiagonalMatrix<T,N>> : std::true_type {};
template <class T>
struct IsVector : std::false_type {};
template <class T, int N>
struct IsVector<FieldVector<T,N>> : std::true_type {};
template <class T>
struct IsMatVec
: std::integral_constant<bool, IsMatrix<T>::value || IsVector<T>::value> {};
/// @}
/// Convert the field_Type of container to type S
/// @{
template <class A, class S>
struct MakeMatVec
{
using type = A;
};
template <class T, int M, int N, class S>
struct MakeMatVec<FieldMatrix<T,M,N>,S>
{
using type = FieldMatrix<S,M,N>;
};
template <class T, int N, class S>
struct MakeMatVec<DiagonalMatrix<T,N>,S>
{
using type = DiagonalMatrix<S,N>;
};
template <class T, int N, class S>
struct MakeMatVec<FieldVector<T,N>,S>
{
using type = FieldVector<S,N>;
};
/// @}
/// Convert pseudo-scalar to real scalar type
/// @{
template <class T>
decltype(auto) simplify(T&& t)
{
return std::forward<T>(t);
}
template <class T>
T simplify(FieldVector<T,1> const& t)
{
return t[0];
}
template <class T>
T simplify(FieldMatrix<T,1,1> const& t)
{
return t[0][0];
}
template <class T>
T simplify(DiagonalMatrix<T,1> const& t)
{
return t.diagonal(0);
}
/// @}
// returns -a
template <class A>
auto negate(A const& a);
// returns a+b
template <class A, class B>
auto plus(A a, B const& b)
{
return a += b;
}
// returns a-b
template <class A, class B>
auto minus(A a, B const& b)
{
return a -= b;
}
// returns a*b
template <class A, class B,
std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int> = 0>
auto multiplies(A const& a, B const& b);
template <class T, int N, class S>
auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b);
template <class Mat, class Vec,
std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int> = 0>
auto multiplies(Mat const& mat, Vec const& vec);
template <class Vec, class Mat,
std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int> = 0>
auto multiplies(Vec const& vec, Mat const& mat);
template <class T, int L, int M, int N, class S>
auto multiplies(FieldMatrix<T,L,M> const& a, FieldMatrix<S,M,N> const& b);
// return a/b
template <class A, class B>
auto divides(A a, B const& b)
{
return a /= b;
}
} // end namespace MatVec
// some arithmetic operations with FieldVector and FieldMatrix
template <class A,
std::enable_if_t<MatVec::IsMatVec<A>::value, int> = 0>
auto operator-(A const& a)
{
return MatVec::negate(MatVec::simplify(a));
}
template <class T, int N> template <class A, class B,
FieldVector<T,N> operator-(FieldVector<T,N> v); std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
auto operator+(A const& a, B const& b)
{
return MatVec::plus(MatVec::simplify(a), MatVec::simplify(b));
}
template <class T, int N, class S, template <class A, class B,
std::enable_if_t<std::is_convertible<S,T>::value, int> = 0> std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
FieldVector<T,N> operator*(FieldVector<T,N> v, S factor); auto operator-(A const& a, B const& b)
{
return MatVec::minus(MatVec::simplify(a), MatVec::simplify(b));
}
template <class T, int N, class S, template <class A, class B,
std::enable_if_t<std::is_convertible<S,T>::value, int> = 0> std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
FieldVector<T,N> operator*(S factor, FieldVector<T,N> v); auto operator*(A const& a, B const& b)
{
return MatVec::multiplies(MatVec::simplify(a), MatVec::simplify(b));
}
template <class T, int N, class S, template <class A, class B,
std::enable_if_t<std::is_convertible<S,T>::value, int> = 0> std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
FieldVector<T,N> operator/(FieldVector<T,N> v, S factor); auto operator/(A const& a, B const& b)
{
return MatVec::divides(MatVec::simplify(a), MatVec::simplify(b));
}
// ---------------------------------------------------------------------------- // ----------------------------------------------------------------------------
...@@ -197,38 +343,8 @@ namespace Dune ...@@ -197,38 +343,8 @@ namespace Dune
template <class T, int M, int N> template <class T, int M, int N>
FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A); FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
template <class T, int M, int N>
FieldMatrix<T,M,N> operator-(FieldMatrix<T, M, N> A);
template <class T, int M, int N, class S,
std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
FieldMatrix<T,M,N> operator*(S scalar, FieldMatrix<T, M, N> A);
template <class T, int M, int N, class S,
std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, S scalar);
template <class T, int M, int N, class S,
std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, S scalar);
template <class T, int M, int N>
FieldMatrix<T,M,N> operator+(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B);
template <class T, int M, int N>
FieldMatrix<T,M,N> operator-(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B);
template <class T, int N, int M>
FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec);
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies(FieldMatrix<T, M, L> const& A, FieldMatrix<T, L, N> const& B);
template <class T, int M, int N, int L> template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B); FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B);
...@@ -243,20 +359,6 @@ namespace Dune ...@@ -243,20 +359,6 @@ namespace Dune
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
template <class T>
T operator*(FieldVector<T, 1> lhs, FieldVector<T, 1> rhs);
template <class T>
T operator*(FieldMatrix<T, 1, 1> lhs, FieldMatrix<T, 1, 1> rhs);
template <class T>
T operator*(FieldVector<T, 1> lhs, FieldMatrix<T, 1, 1> rhs);
template <class T>
T operator*(FieldMatrix<T, 1, 1> lhs, FieldVector<T, 1> rhs);
// -----------------------------------------------------------------------------
template <class T, int N> template <class T, int N>
T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i); T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i);
......
...@@ -11,34 +11,81 @@ ...@@ -11,34 +11,81 @@
namespace Dune { namespace Dune {
// some arithmetic operations with FieldVector namespace MatVec {
template <class T, int N> template <class A>
FieldVector<T,N> operator-(FieldVector<T,N> v) auto negate(A const& a)
{ {
return v *= -1; return multiplies(a, -1);
} }
template <class T, int N, class S, template <class A, class B,
std::enable_if_t<std::is_convertible<S,T>::value, int>> std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int>>
FieldVector<T,N> operator*(FieldVector<T,N> v, S factor) auto multiplies(A const& a, B const& b)
{ {
return v *= factor; using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
} #if AMDIS_HAS_CXX_CONSTEXPR_IF
if constexpr(IsNumber<A>::value) {
typename MakeMatVec<B,T>::type b_{b};
return b_ *= a;
} else {
typename MakeMatVec<A,T>::type a_{a};
return a_ *= b;
}
#else
return Hybrid::ifElse(IsNumber<A>{},
[&](auto id) { typename MakeMatVec<B,T>::type b_{b}; return id(b_) *= id(a); },
[&](auto id) { typename MakeMatVec<A,T>::type a_{a}; return id(a_) *= id(b); });
#endif
}
template <class T, int N, class S, template <class T, int N, class S>
std::enable_if_t<std::is_convertible<S,T>::value, int>> auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b)
FieldVector<T,N> operator*(S factor, FieldVector<T,N> v) {
{ return a.dot(b);
return v *= factor; }
}
template <class T, int N, class S,
std::enable_if_t<std::is_convertible<S,T>::value, int>> template <class Mat, class Vec,
FieldVector<T,N> operator/(FieldVector<T,N> v, S factor) std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int>>
{ auto multiplies(Mat const& mat, Vec const& vec)
return v /= factor; {
} static_assert(Mat::cols == Vec::dimension, "");
using T = std::common_type_t<typename FieldTraits<Vec>::field_type, typename FieldTraits<Mat>::field_type>;
FieldVector<T,Mat::rows> y;
mat.mv(vec, y);
return y;
}
template <class Vec, class Mat,
std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int>>
auto multiplies(Vec const& vec, Mat const& mat)
{
static_assert(Mat::rows == Vec::dimension, "");
using T = std::common_type_t<typename FieldTraits<Vec>::field_type, typename FieldTraits<Mat>::field_type>;
FieldVector<T,Mat::cols> y;
mat.mtv(vec, y);
return y;
}
template <class T, int L, int M, int N, class S>
auto multiplies(FieldMatrix<T,L,M> const& a, FieldMatrix<S,M,N> const& b)
{
FieldMatrix<std::common_type_t<T,S>,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<T,N,M> trans(FieldMatrix<T, M, N> const& A) ...@@ -349,57 +396,6 @@ FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
return At; return At;
} }
template <class T, int M, int N>
FieldMatrix<T,M,N> operator-(FieldMatrix<T,M,N> A)
{
return A *= -1;
}
template <class T, int M, int N, class S,
std::enable_if_t<std::is_convertible<S,T>::value, int>>
FieldMatrix<T,M,N> operator*(S scalar, FieldMatrix<T, M, N> A)
{
return A *= scalar;
}
template <class T, int M, int N, class S,
std::enable_if_t<std::is_convertible<S,T>::value, int>>
FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, S scalar)
{
return A *= scalar;
}
template <class T, int M, int N, class S,
std::enable_if_t<std::is_convertible<S,T>::value, int>>
FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, S scalar)
{
return A /= scalar;
}
template <class T, int M, int N>
FieldMatrix<T,M,N> operator+(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
{
return A += B;
}
template <class T, int M, int N>
FieldMatrix<T,M,N> operator-(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
{
return A -= B;
}
template <class T, int N, int M>
FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec)
{
return Dune::FMatrixHelp::mult(mat, vec);
}
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies(FieldMatrix<T, M, L> const& A, FieldMatrix<T, L, N> const& B)
{
return A.rightmultiplyany(B);
}
template <class T, int M, int N, int L> template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B) FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B)
...@@ -447,32 +443,6 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, DiagonalMatri ...@@ -447,32 +443,6 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, DiagonalMatri
return C; return C;
} }
template <class T>
T operator*(FieldVector<T,1> lhs, FieldVector<T,1> rhs)
{
return lhs[0] * rhs[0];
}
template <class T>
T operator*(FieldMatrix<T,1,1> lhs, FieldMatrix<T,1,1> rhs)
{
return lhs[0][0] * rhs[0][0];
}
template <class T>
T operator*(FieldVector<T,1> lhs, FieldMatrix<T,1,1> rhs)
{
return lhs[0] * rhs[0][0];
}
template <class T>
T operator*(FieldMatrix<T,1,1> lhs, FieldVector<T,1> rhs)
{
return lhs[0][0] * rhs[0];
}
template <class T, int N> template <class T, int N>
T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i) T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i)
{ {
......
...@@ -101,7 +101,7 @@ void test2() ...@@ -101,7 +101,7 @@ void test2()
AMDIS_TEST_EQ( det(A), 1.0 ); AMDIS_TEST_EQ( det(A), 1.0 );
AMDIS_TEST_EQ( det(B), 0.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(trans(C)), det(C) );
AMDIS_TEST_EQ( det(2.0*D), Math::pow<3>(2.0)*det(D) ); AMDIS_TEST_EQ( det(2.0*D), Math::pow<3>(2.0)*det(D) );
...@@ -201,6 +201,15 @@ void test3() ...@@ -201,6 +201,15 @@ void test3()
AMDIS_TEST_EQ( vi1*mi1, 1 ); AMDIS_TEST_EQ( vi1*mi1, 1 );
AMDIS_TEST_EQ( mi1*mi1, 1 ); AMDIS_TEST_EQ( mi1*mi1, 1 );
AMDIS_TEST_EQ( mi1*vi1, 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 );
} }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment