wrappers for vector and matrix applied to DiscreteFunction and FOT_partial

src/amdis/common/FieldMatVec.hpp

@@ -127,6 +127,24 @@ namespace Dune
     }
     /// @}
 
+
+    /// 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);

src/amdis/common/FieldMatVec.inc.hpp

@@ -126,6 +126,49 @@ namespace MatVec {
     return C;
   }
 
+  template <class T, int N>
+  class MatrixView<DiagonalMatrix<T,N>>
+  {
+    using Matrix = DiagonalMatrix<T,N>;
+    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(DiagonalMatrix<T,N> 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:
+    DiagonalMatrix<T,N> const& mat_;
+  };
+
 } // end namespace MatVec
 
 // ----------------------------------------------------------------------------

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/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 = 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 = 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) {