Merge branch 'issue/mixed_type_fieldmatvec' into 'develop'

Redesign of FieldMatrix + FieldVector arithmetic operations

See merge request spraetor/dune-amdis!83

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 <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>
-  FieldVector<T,N> operator-(FieldVector<T,N> v);
+  template <class A, class B,
+    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,
-    std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
-  FieldVector<T,N> operator*(FieldVector<T,N> v, S factor);
+  template <class A, class B,
+    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
+  auto operator-(A const& a, B const& b)
+  {
+    return MatVec::minus(MatVec::simplify(a), MatVec::simplify(b));
+  }
 
-  template <class T, int N, class S,
-    std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
-  FieldVector<T,N> operator*(S factor, FieldVector<T,N> v);
+  template <class A, class B,
+    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
+  auto operator*(A const& a, B const& b)
+  {
+    return MatVec::multiplies(MatVec::simplify(a), MatVec::simplify(b));
+  }
 
-  template <class T, int N, class S,
-    std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
-  FieldVector<T,N> operator/(FieldVector<T,N> v, S factor);
+  template <class A, class B,
+    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::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 <class T, int M, int N>
   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>
   FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A,  FieldMatrix<T, N, L> const& B);
 
@@ -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>
   T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i);
 

src/amdis/common/FieldMatVec.inc.hpp

@@ -11,34 +11,81 @@
 
 namespace Dune {
 
-// some arithmetic operations with FieldVector
+namespace MatVec {
 
-template <class T, int N>
-FieldVector<T,N> operator-(FieldVector<T,N> v)
-{
-  return v *= -1;
-}
+  template <class A>
+  auto negate(A const& a)
+  {
+    return multiplies(a, -1);
+  }
 
-template <class T, int N, class S,
-  std::enable_if_t<std::is_convertible<S,T>::value, int>>
-FieldVector<T,N> operator*(FieldVector<T,N> v, S factor)
-{
-  return v *= factor;
-}
+  template <class A, class B,
+    std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int>>
+  auto multiplies(A const& a, B const& b)
+  {
+    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,
-  std::enable_if_t<std::is_convertible<S,T>::value, int>>
-FieldVector<T,N> operator*(S factor, FieldVector<T,N> v)
-{
-  return v *= factor;
-}
+  template <class T, int N, class S>
+  auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b)
+  {
+    return a.dot(b);
+  }
 
-template <class T, int N, class S,
-  std::enable_if_t<std::is_convertible<S,T>::value, int>>
-FieldVector<T,N> operator/(FieldVector<T,N> v, S factor)
-{
-  return v /= factor;
-}
+
+  template <class Mat, class Vec,
+    std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int>>
+  auto multiplies(Mat const& mat, Vec const& vec)
+  {
+    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)
   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>
 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
   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>
 T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i)
 {

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 );
 }