Merge branch 'feature/more_fieldmatvec_operations' into 'develop'

Added some missing FieldVector and FieldMatrix operations

See merge request !62

src/amdis/common/FieldMatVec.hpp

@@ -6,31 +6,40 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
+namespace std
+{
+  template <class T, int N>
+  struct common_type<Dune::FieldVector<T,N>, T>
+  {
+    using type = T;
+  };
+
+  template <class T, int N, int M>
+  struct common_type<Dune::FieldMatrix<T,N,M>, T>
+  {
+    using type = T;
+  };
+}
+
 namespace Dune
 {
   // some arithmetic operations with FieldVector
 
+  template <class T, int N>
+  FieldVector<T,N> operator-(FieldVector<T,N> v);
+
   template <class T, int N, class S,
-    std::enable_if_t<std::is_arithmetic<S>::value,int> = 0 >
+    std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
   FieldVector<T,N> operator*(FieldVector<T,N> v, S factor);
 
-  template <class S, class T, int N,
-    std::enable_if_t<std::is_arithmetic<S>::value,int> = 0 >
+  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 T, int N, class S,
-    std::enable_if_t<std::is_arithmetic<S>::value,int> = 0 >
+    std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
   FieldVector<T,N> operator/(FieldVector<T,N> v, S factor);
 
-  template <class T>
-  FieldVector<T,1> operator*(FieldVector<T,1> const& v, FieldVector<T,1> const& w);
-
-  template <class T, int N>
-  FieldVector<T,N> operator*(FieldVector<T,1> const& factor, FieldVector<T,N> v);
-
-  template <class T, int N>
-  FieldVector<T,N> operator*(FieldVector<T,N> v, FieldVector<T,1> const& factor);
-
   // ----------------------------------------------------------------------------
 
   /// Cross-product a 2d-vector = orthogonal vector
@@ -45,10 +54,6 @@ namespace Dune
   template <class T, class S, int N>
   auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2);
 
-  template <class T, int N, int M,
-    std::enable_if_t<( N!=1 && M!=1 ),int> = 0>
-  auto operator*(FieldVector<T,N> const& v, FieldVector<T,M> const& w);
-
   template <class T, class S, int N>
   auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2);
 
@@ -194,15 +199,15 @@ namespace Dune
 
 
   template <class T, int M, int N, class S,
-    std::enable_if_t<std::is_arithmetic<S>::value,int> = 0 >
+    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_arithmetic<S>::value,int> = 0 >
+    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_arithmetic<S>::value,int> = 0 >
+    std::enable_if_t<std::is_convertible<S,T>::value, int> = 0>
   FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, S scalar);
 
 
@@ -216,19 +221,11 @@ namespace Dune
   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 N, int M>
-  FieldMatrix<T,N,M> operator*(FieldMatrix<T,N,M> mat, FieldVector<T,1> const& scalar);
-
-  template <class T, int N>
-  FieldMatrix<T,N,1> operator*(FieldMatrix<T,N,1> mat, FieldVector<T,1> const& scalar);
-
-
+  // -----------------------------------------------------------------------------
 
   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);
 
@@ -241,6 +238,35 @@ namespace Dune
   template <class T, int M, int N>
   FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A,  DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C);
 
+  // -----------------------------------------------------------------------------
+
+  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);
+
+  template <class T, int M>
+  T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i);
+
+  // necessary specialization to resolve ambiguous calls
+  template <class T>
+  T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i);
+
+  template <class T, int N>
+  T const& at(FieldVector<T,N> const& vec, std::size_t i);
+
 } // end namespace Dune
 
 namespace AMDiS

src/amdis/common/FieldMatVec.inc.hpp

@@ -13,45 +13,33 @@ namespace Dune {
 
 // some arithmetic operations with FieldVector
 
+template <class T, int N>
+FieldVector<T,N> operator-(FieldVector<T,N> v)
+{
+  return v *= -1;
+}
+
 template <class T, int N, class S,
-  std::enable_if_t<std::is_arithmetic<S>::value,int> >
+  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 S, class T, int N,
-  std::enable_if_t<std::is_arithmetic<S>::value,int> >
+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,
-  std::enable_if_t<std::is_arithmetic<S>::value,int> >
+  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 T>
-FieldVector<T,1> operator*(FieldVector<T,1> const& v, FieldVector<T,1> const& w)
-{
-  return {v[0] * w[0]};
-}
-
-template <class T, int N>
-FieldVector<T,N> operator*(FieldVector<T,1> const& factor, FieldVector<T,N> v)
-{
-  return v *= factor[0];
-}
-
-template <class T, int N>
-FieldVector<T,N> operator*(FieldVector<T,N> v, FieldVector<T,1> const& factor)
-{
-  return v *= factor[0];
-}
-
 // ----------------------------------------------------------------------------
 
 /// Cross-product a 2d-vector = orthogonal vector
@@ -77,14 +65,6 @@ auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2)
   return vec1.dot(vec2);
 }
 
-template <class T, int N, int M,
-  std::enable_if_t<( N!=1 && M!=1 ),int> >
-auto operator*(FieldVector<T,N> const& v, FieldVector<T,M> const& w)
-{
-  static_assert(M == N, "Requires vectors of the same type!");
-  return v.dot(w);
-}
-
 template <class T, class S, int N>
 auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2)
 {
@@ -371,21 +351,21 @@ FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
 
 
 template <class T, int M, int N, class S,
-  std::enable_if_t<std::is_arithmetic<S>::value,int> >
+  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_arithmetic<S>::value,int> >
+  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_arithmetic<S>::value,int> >
+  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;
@@ -410,28 +390,12 @@ FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const
   return Dune::FMatrixHelp::mult(mat, vec);
 }
 
-template <class T, int N, int M>
-FieldMatrix<T,N,M> operator*(FieldMatrix<T,N,M> mat, FieldVector<T,1> const& scalar)
-{
-  return mat *= scalar[0];
-}
-
-template <class T, int N>
-FieldMatrix<T,N,1> operator*(FieldMatrix<T,N,1> mat, FieldVector<T,1> const& scalar)
-{
-  return mat *= scalar[0];
-}
-
-
-
 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)
 {
@@ -451,15 +415,7 @@ template <class T, int M, int N, int L>
 FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T, N, L> const& B)
 {
   FieldMatrix<T,M,N> C;
-
-  for (int m = 0; m < M; ++m) {
-    for (int n = 0; n < N; ++n) {
-      C[m][n] = 0;
-      for (int l = 0; l < L; ++l)
-        C[m][n] += A[m][l] * B[n][l];
-    }
-  }
-  return C;
+  return multiplies_ABt(A,B,C);
 }
 
 template <class T, int M, int N, int L>
@@ -486,6 +442,56 @@ 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)
+{
+  return vec[i][0];
+}
+
+template <class T, int M>
+T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i)
+{
+  return vec[0][i];
+}
+
+template <class T>
+T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i)
+{
+  return vec[0][i];
+}
+
+template <class T, int N>
+T const& at(FieldVector<T,N> const& vec, std::size_t i)
+{
+  return vec[i];
+}
+
 } // end namespace AMDiS
 
 #endif

test/FieldMatVecTest.cpp

@@ -128,12 +128,88 @@ void test2()
   AMDIS_TEST_EQ( sol, a );
 }
 
+// test of scalar wrapper FieldVector<T,1> and FieldMatrix<T,1,1>
+void test3()
+{
+  using V = FieldVector<double, 3>;
+  using M = FieldMatrix<double, 3, 3>;
+
+  V a{1.0, 2.0, 3.0};
+  M A{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} };
+
+  using VD = FieldVector<double, 1>;
+  using MD = FieldMatrix<double, 1, 1>;
+
+  VD vd1 = 1.0;
+  MD md1 = 1.0;
+
+  using VI = FieldVector<int, 1>;
+  using MI = FieldMatrix<int, 1, 1>;
+
+  VI vi1 = 1;
+  MI mi1 = 1;
+
+  // scale a vector
+  AMDIS_TEST_EQ( 1*a, a );
+  AMDIS_TEST_EQ( 1.0*a, a );
+  AMDIS_TEST_EQ( a*1, a );
+  AMDIS_TEST_EQ( a*1.0, a );
+  AMDIS_TEST_EQ( vd1*a, a );
+  AMDIS_TEST_EQ( a*vd1, a );
+  AMDIS_TEST_EQ( vi1*a, a );
+  AMDIS_TEST_EQ( a*vi1, a );
+  AMDIS_TEST_EQ( md1*a, a );
+  AMDIS_TEST_EQ( a*md1, a );
+  AMDIS_TEST_EQ( mi1*a, a );
+  AMDIS_TEST_EQ( a*mi1, a );
+  AMDIS_TEST_EQ( a/1, a );
+  AMDIS_TEST_EQ( a/1.0, a );
+  AMDIS_TEST_EQ( a/vd1, a );
+  AMDIS_TEST_EQ( a/vi1, a );
+  AMDIS_TEST_EQ( a/md1, a );
+  AMDIS_TEST_EQ( a/mi1, a );
+
+  // scale a matrix
+  AMDIS_TEST_EQ( 1*A, A );
+  AMDIS_TEST_EQ( 1.0*A, A );
+  AMDIS_TEST_EQ( A*1, A );
+  AMDIS_TEST_EQ( A*1.0, A );
+  AMDIS_TEST_EQ( vd1*A, A );
+  AMDIS_TEST_EQ( A*vd1, A );
+  AMDIS_TEST_EQ( vi1*A, A );
+  AMDIS_TEST_EQ( A*vi1, A );
+  AMDIS_TEST_EQ( md1*A, A );
+  AMDIS_TEST_EQ( A*md1, A );
+  AMDIS_TEST_EQ( mi1*A, A );
+  AMDIS_TEST_EQ( A*mi1, A );
+  AMDIS_TEST_EQ( A/1, A );
+  AMDIS_TEST_EQ( A/1.0, A );
+  AMDIS_TEST_EQ( A/vd1, A );
+  AMDIS_TEST_EQ( A/vi1, A );
+  AMDIS_TEST_EQ( A/md1, A );
+  AMDIS_TEST_EQ( A/mi1, A );
+
+  AMDIS_TEST_EQ( vd1*vd1, 1.0 );
+  AMDIS_TEST_EQ( vd1*md1, 1.0 );
+  AMDIS_TEST_EQ( md1*md1, 1.0 );
+  AMDIS_TEST_EQ( md1*vd1, 1.0 );
+  AMDIS_TEST_EQ( vd1*1.0, 1.0 );
+  AMDIS_TEST_EQ( md1*1.0, 1.0 );
+  AMDIS_TEST_EQ( 1.0*md1, 1.0 );
+  AMDIS_TEST_EQ( 1.0*vd1, 1.0 );
+  AMDIS_TEST_EQ( vi1*vi1, 1 );
+  AMDIS_TEST_EQ( vi1*mi1, 1 );
+  AMDIS_TEST_EQ( mi1*mi1, 1 );
+  AMDIS_TEST_EQ( mi1*vi1, 1 );
+}
+
 
 int main()
 {
   test0();
   test1();
   test2();
+  test3();
 
   return report_errors();
 }