Merge branch 'feature/fixmatvec_operations' into example/stokes

src/amdis/DirichletBC.inc.hpp

@@ -7,8 +7,7 @@
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/functionspacebases/subspacebasis.hh>
 
-#include <amdis/Output.hpp>
-#include <amdis/linearalgebra/Constraints.hpp>
+#include <amdis/LinearAlgebra.hpp>
 
 namespace AMDiS {
 

src/amdis/common/FieldMatVec.hpp

@@ -82,51 +82,88 @@ 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];
+    }
+    /// @}
+
+
+    /// Convert pseudo-matrix to real matrix type with proper operator[][]
+    /// @{
+    template <class T>
+    decltype(auto) as_matrix(T&& t)
+    {
+      return FWD(t);
+    }
+
+    template <class Mat>
+    class MatrixView;
+
+    template <class T, int N>
+    MatrixView<DiagonalMatrix<T,N>> as_matrix(DiagonalMatrix<T,N> const& mat)
+    {
+      return {mat};
+    }
+
     // 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 +195,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 +383,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);
 
   // -----------------------------------------------------------------------------
 

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>;
@@ -85,6 +126,49 @@ namespace MatVec {
     return C;
   }
 
+  template <class T, int SIZE>
+  class MatrixView<DiagonalMatrix<T,SIZE>>
+  {
+    using Matrix = DiagonalMatrix<T,SIZE>;
+    using size_type = typename Matrix::size_type;
+
+    struct RowView
+    {
+      T operator[](size_type c) const
+      {
+        assert(0 <= c && c < mat_->M());
+        return c == r_ ? mat_->diagonal(r_) : T(0);
+      }
+
+      Matrix const* mat_;
+      size_type r_;
+    };
+
+  public:
+    MatrixView(Matrix const& mat)
+      : mat_(mat)
+    {}
+
+    RowView operator[](size_type r) const
+    {
+      assert(0 <= r && r < mat_.N());
+      return {&mat_,r};
+    }
+
+    size_type N() const
+    {
+      return mat_.N();
+    }
+
+    size_type M() const
+    {
+      return mat_.M();
+    }
+
+  private:
+    Matrix const& mat_;
+  };
+
 } // end namespace MatVec
 
 // ----------------------------------------------------------------------------
@@ -158,14 +242,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 +313,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 +345,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 +377,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 +482,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 +497,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 +517,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 +528,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)
 {

src/amdis/gridfunctions/DiscreteFunction.inc.hpp

@@ -3,6 +3,7 @@
 #include <amdis/common/FieldMatVec.hpp>
 #include <amdis/utility/LocalBasisCache.hpp>
 
+#include <dune/common/ftraits.hh>
 #include <dune/grid/utility/hierarchicsearch.hh>
 
 namespace AMDiS {
@@ -197,7 +198,7 @@ GradientLocalFunction::operator()(Domain const& x) const
   {
     // TODO: may DOFVectorView::Range to FieldVector type if necessary
     using LocalDerivativeTraits
-      = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
+      = Dune::Functions::DefaultDerivativeTraits<VT(Domain)>;
     using GradientBlock = typename LocalDerivativeTraits::Range;
 
     auto&& fe = node.finiteElement();
@@ -211,8 +212,8 @@ GradientLocalFunction::operator()(Domain const& x) const
 
     // Compute the shape function gradients on the real element
     std::vector<GradientBlock> gradients(referenceGradients.size());
-    for (std::size_t i = 0; i < gradients.size(); ++i)
-      multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad
+    for (std::size_t i = 0; i < gradients.size(); ++i) // J^(-T) * D[phi] -> grad^T
+      Dune::MatVec::as_vector(gradients[i]) = jacobian * Dune::MatVec::as_vector(referenceGradients[i]);
 
     // Get range entry associated to this node
     auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));

src/amdis/linearalgebra/mtl/Constraints.hpp

@@ -130,8 +130,7 @@ namespace AMDiS
         ins[t.row][t.col] += t.value;
 
       if (setDiagonal) {
-        using mtl::size;
-        for (std::size_t i = 0; i < size(left); ++i) {
+        for (std::size_t i = 0; i < mtl::size(left); ++i) {
           if (left[i]) {
             ins[i][i] = T(1);
             ins[i][at(left2right,i)] = T(-1);

src/amdis/localoperators/FirstOrderTestPartialTrial.hpp

@@ -69,6 +69,7 @@ namespace AMDiS
 
         // The transposed inverse Jacobian of the map from the reference element to the element
         const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
+        auto&& jacobian_mat = Dune::MatVec::as_matrix(jacobian);
         assert(jacobian.M() == CG::dim);
 
         // The multiplicative factor in the integral transformation formula
@@ -85,9 +86,9 @@ namespace AMDiS
         thread_local std::vector<RangeFieldType> colPartial;
         colPartial.resize(shapeGradients.size());
         for (std::size_t i = 0; i < colPartial.size(); ++i) {
-          colPartial[i] = jacobian[comp_][0] * shapeGradients[i][0][0];
-          for (std::size_t j = 1; j < jacobian.M(); ++j)
-            colPartial[i] += jacobian[comp_][j] * shapeGradients[i][0][j];
+          colPartial[i] = jacobian_mat[comp_][0] * shapeGradients[i][0][0];
+          for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
+            colPartial[i] += jacobian_mat[comp_][j] * shapeGradients[i][0][j];
         }
 
         for (std::size_t j = 0; j < colSize; ++j) {

src/amdis/localoperators/SecondOrderPartialTestPartialTrial.hpp

@@ -67,6 +67,7 @@ namespace AMDiS
 
         // The transposed inverse Jacobian of the map from the reference element to the element
         const auto jacobian = contextGeo.geometry().jacobianInverseTransposed(local);
+        auto&& jacobian_mat = Dune::MatVec::as_matrix(jacobian);
 
         // The multiplicative factor in the integral transformation formula
         const auto factor = contextGeo.integrationElement(quad[iq].position()) * quad[iq].weight();
@@ -80,17 +81,17 @@ namespace AMDiS
         thread_local std::vector<RowFieldType> rowPartial;
         rowPartial.resize(rowShapeGradients.size());
         for (std::size_t i = 0; i < rowPartial.size(); ++i) {
-          rowPartial[i] = jacobian[compTest_][0] * rowShapeGradients[i][0][0];
-          for (std::size_t j = 1; j < jacobian.cols; ++j)
-            rowPartial[i] += jacobian[compTest_][j] * rowShapeGradients[i][0][j];
+          rowPartial[i] = jacobian_mat[compTest_][0] * rowShapeGradients[i][0][0];
+          for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
+            rowPartial[i] += jacobian_mat[compTest_][j] * rowShapeGradients[i][0][j];
         }
 
         thread_local std::vector<ColFieldType> colPartial;
         colPartial.resize(colShapeGradients.size());
         for (std::size_t i = 0; i < colPartial.size(); ++i) {
-          colPartial[i] = jacobian[compTrial_][0] * colShapeGradients[i][0][0];
-          for (std::size_t j = 1; j < jacobian.cols; ++j)
-            colPartial[i] += jacobian[compTrial_][j] * colShapeGradients[i][0][j];
+          colPartial[i] = jacobian_mat[compTrial_][0] * colShapeGradients[i][0][0];
+          for (std::size_t j = 1; j < jacobian_mat.M(); ++j)
+            colPartial[i] += jacobian_mat[compTrial_][j] * colShapeGradients[i][0][j];
         }
 
         for (std::size_t j = 0; j < colSize; ++j) {