From 6efe929728b636dd798fbcb7e1019dc99df22660 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 20 Mar 2019 16:26:50 +0100
Subject: [PATCH] renamed simplify to as_scalar and add as_vector

---
 src/amdis/common/FieldMatVec.hpp     | 80 ++++++++++++++++---------
 src/amdis/common/FieldMatVec.inc.hpp | 87 ++++++++++++++++++++++------
 2 files changed, 123 insertions(+), 44 deletions(-)

diff --git a/src/amdis/common/FieldMatVec.hpp b/src/amdis/common/FieldMatVec.hpp
index 8c4aeb87..f2b8b053 100644
--- a/src/amdis/common/FieldMatVec.hpp
+++ b/src/amdis/common/FieldMatVec.hpp
@@ -82,51 +82,70 @@ namespace Dune
     /// Convert pseudo-scalar to real scalar type
     /// @{
     template <class T>
-    decltype(auto) simplify(T&& t)
+    decltype(auto) as_scalar(T&& t)
     {
       return FWD(t);
     }
 
     template <class T>
-    T simplify(FieldVector<T,1> const& t)
+    T as_scalar(FieldVector<T,1> const& t)
     {
       return t[0];
     }
 
     template <class T>
-    T simplify(FieldMatrix<T,1,1> const& t)
+    T as_scalar(FieldMatrix<T,1,1> const& t)
     {
       return t[0][0];
     }
 
     template <class T>
-    T simplify(DiagonalMatrix<T,1> const& t)
+    T as_scalar(DiagonalMatrix<T,1> const& t)
     {
       return t.diagonal(0);
     }
     /// @}
 
+    /// Convert pseudo-vector to real vector type
+    /// @{
+    template <class T>
+    decltype(auto) as_vector(T&& t)
+    {
+      return FWD(t);
+    }
+
+    template <class T, int N>
+    FieldVector<T,N> const& as_vector(FieldMatrix<T,1,N> const& t)
+    {
+      return t[0];
+    }
+
+    template <class T, int N>
+    FieldVector<T,N>& as_vector(FieldMatrix<T,1,N>& t)
+    {
+      return t[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;
-    }
+    auto plus(A const& a, B const& b);
 
     // returns a-b
     template <class A, class B>
-    auto minus(A a, B const& b)
-    {
-      return a -= b;
-    }
+    auto minus(A const& a, B const& b);
 
     // returns a*b
     template <class A, class B,
-      std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int> = 0>
+      std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int> = 0>
+    auto multiplies(A const& a, B const& 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>
@@ -158,35 +177,35 @@ namespace Dune
     std::enable_if_t<MatVec::IsMatVec<A>::value, int> = 0>
   auto operator-(A const& a)
   {
-    return MatVec::negate(MatVec::simplify(a));
+    return MatVec::negate(MatVec::as_scalar(a));
   }
 
   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));
+    return MatVec::plus(MatVec::as_scalar(a), MatVec::as_scalar(b));
   }
 
   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));
+    return MatVec::minus(MatVec::as_scalar(a), MatVec::as_scalar(b));
   }
 
   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));
+    return MatVec::multiplies(MatVec::as_scalar(a), MatVec::as_scalar(b));
   }
 
   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));
+    return MatVec::divides(MatVec::as_scalar(a), MatVec::as_scalar(b));
   }
 
   // ----------------------------------------------------------------------------
@@ -346,19 +365,28 @@ namespace Dune
   template <class T, int M, int N>
   FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
 
+  template <class T, int N>
+  DiagonalMatrix<T,N> const& trans(DiagonalMatrix<T,N> const& A)
+  {
+    return A;
+  }
+
   // -----------------------------------------------------------------------------
 
-  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);
+  template <class T1, class T2, int M, int N, int L>
+  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A,  FieldMatrix<T2, N, L> const& B);
 
-  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);
+  template <class T1, class T2, int M, int N, int L>
+  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B);
 
-  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);
+  template <class T1, class T2, class T3, int M, int N, int L>
+  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C);
 
-  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 T1, class T2, class T3, int M, int N>
+  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C);
+
+  template <class T1, class T2, class T3, int N>
+  FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);
 
   // -----------------------------------------------------------------------------
 
diff --git a/src/amdis/common/FieldMatVec.inc.hpp b/src/amdis/common/FieldMatVec.inc.hpp
index 4644778c..4e1e0220 100644
--- a/src/amdis/common/FieldMatVec.inc.hpp
+++ b/src/amdis/common/FieldMatVec.inc.hpp
@@ -3,6 +3,8 @@
 #include <algorithm>
 #include <limits>
 
+#include <dune/functions/functionspacebases/flatvectorview.hh>
+
 #include <amdis/common/Math.hpp>
 #include <amdis/operations/Arithmetic.hpp>
 #include <amdis/operations/MaxMin.hpp>
@@ -19,8 +21,47 @@ namespace MatVec {
     return multiplies(a, -1);
   }
 
+  // returns a+b
+  template <class A, class B>
+  auto plus(A const& a, B const& b)
+  {
+    using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
+    typename MakeMatVec<A,T>::type c{a};
+
+    auto b_ = Dune::Functions::flatVectorView(b);
+    auto c_ = Dune::Functions::flatVectorView(c);
+    assert(int(b_.size()) == int(c_.size()));
+    for(int i = 0; i < int(b_.size()); ++i)
+      c_[i] += b_[i];
+
+    return c;
+  }
+
+  // returns a-b
+  template <class A, class B>
+  auto minus(A const& a, B const& b)
+  {
+    using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
+    typename MakeMatVec<A,T>::type c{a};
+
+    auto b_ = Dune::Functions::flatVectorView(b);
+    auto c_ = Dune::Functions::flatVectorView(c);
+    assert(int(b_.size()) == int(c_.size()));
+    for(int i = 0; i < int(b_.size()); ++i)
+      c_[i] -= b_[i];
+
+    return c;
+  }
+
   template <class A, class B,
-    std::enable_if_t<IsNumber<A>::value || IsNumber<B>::value, int>>
+    std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int>>
+  auto multiplies(A const& a, B const& b)
+  {
+    return a * b;
+  }
+
+  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>;
@@ -158,14 +199,14 @@ T sum(FieldMatrix<T, 1, N> const& x)
 template <class T, int N>
 auto unary_dot(FieldVector<T, N> const& x)
 {
-  auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::sqr(std::abs(b)); };
+  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
   return Impl::accumulate(x, T(0), op);
 }
 
 template <class T, int N>
 auto unary_dot(FieldMatrix<T, 1, N> const& x)
 {
-  auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::sqr(std::abs(b)); };
+  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
   return Impl::accumulate(x, T(0), op);
 }
 
@@ -229,14 +270,14 @@ auto abs_min(FieldMatrix<T, 1, N> const& x)
 template <class T, int N>
 auto one_norm(FieldVector<T, N> const& x)
 {
-  auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
+  auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
   return Impl::accumulate(x, T(0), op);
 }
 
 template <class T, int N>
 auto one_norm(FieldMatrix<T, 1, N> const& x)
 {
-  auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
+  auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
   return Impl::accumulate(x, T(0), op);
 }
 
@@ -261,14 +302,14 @@ auto two_norm(FieldMatrix<T, 1, N> const& x)
 template <int p, class T, int N>
 auto p_norm(FieldVector<T, N> const& x)
 {
-  auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::pow<p>(std::abs(b)); };
+  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
   return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
 }
 
 template <int p, class T, int N>
 auto p_norm(FieldMatrix<T, 1, N> const& x)
 {
-  auto op = [](auto const& a, auto const& b) { return a + AMDiS::Math::pow<p>(std::abs(b)); };
+  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
   return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
 }
 
@@ -293,10 +334,11 @@ auto infty_norm(FieldMatrix<T, 1, N> const& x)
 template <class T, int N>
 T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs)
 {
+  using std::sqrt;
   T result = 0;
   for (int i = 0; i < N; ++i)
     result += AMDiS::Math::sqr(lhs[i] - rhs[i]);
-  return std::sqrt(result);
+  return sqrt(result);
 }
 
 // ----------------------------------------------------------------------------
@@ -397,10 +439,10 @@ FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
 }
 
 
-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)
+template <class T1, class T2, int M, int N, int L>
+FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A,  FieldMatrix<T2, N, L> const& B)
 {
-  FieldMatrix<T,M,N> C;
+  FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
 
   for (int m = 0; m < M; ++m) {
     for (int n = 0; n < N; ++n) {
@@ -412,15 +454,15 @@ FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A,  FieldMatrix<T,
   return C;
 }
 
-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)
+template <class T1, class T2, int M, int N, int L>
+FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B)
 {
-  FieldMatrix<T,M,N> C;
+  FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
   return multiplies_ABt(A,B,C);
 }
 
-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)
+template <class T1, class T2, class T3, int M, int N, int L>
+FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C)
 {
   for (int m = 0; m < M; ++m) {
     for (int n = 0; n < N; ++n) {
@@ -432,8 +474,8 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A,  FieldMatrix<T
   return C;
 }
 
-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 T1, class T2, class T3, int M, int N>
+FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C)
 {
   for (int m = 0; m < M; ++m) {
     for (int n = 0; n < N; ++n) {
@@ -443,6 +485,15 @@ FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A,  DiagonalMatri
   return C;
 }
 
+template <class T1, class T2, class T3, int N>
+FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C)
+{
+  for (int n = 0; n < N; ++n) {
+    C[n] = A[0][n] * B.diagonal(n);
+  }
+  return C;
+}
+
 template <class T, int N>
 T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i)
 {
-- 
GitLab