Use SingleNonZeroRowMatrix implementation that was in ComponentwiseMatrixMap.

dune/matrix-vector/singlenonzerorowmatrix.hh

@@ -8,7 +8,6 @@
 
 #include <dune/fufem/arithmetic.hh>
 #include <dune/fufem/indexedsliceiterator.hh>
-
 /**
  * \brief A static matrix that has only a single nonzero row
  *
@@ -17,39 +16,81 @@
  * all entries accordingly. So it will not reduce memory
  * requirements but allows to implement methods more efficiently.
  */
-template<class K, int ROWS, int COLS>
-class SingleNonZeroRowMatrix :
-    public Dune::FieldMatrix<K, ROWS, COLS>
+template<class K, int ROWS, int COLS> //TODO allow to set storage type as in singlenonzerocolumnmatrix. Update docu accordingly.
+class SingleNonZeroRowMatrix
 {
-    typedef typename Dune::FieldMatrix<K, ROWS, COLS> Base;
+public:
+    typedef typename Dune::FieldMatrix<K, 1, COLS> SingleRowMatrix;
+
+protected:
+    class RowProxy
+    {
+    public:
+        typedef std::size_t size_type;
+        typedef typename SingleRowMatrix::ConstColIterator ConstIterator;
+        typedef ConstIterator const_iterator;
+
+        RowProxy(const SingleRowMatrix* nzRow, size_type rowIndex, size_type nzRowIndex) :
+            nzRow_(nzRow),
+            isNonZero_(rowIndex==nzRowIndex)
+        {}
+
+        ConstIterator begin() const
+        {
+            if (isNonZero_)
+                return (*nzRow_)[0].begin();
+            else
+                return (*nzRow_)[0].end();
+        }
+
+        ConstIterator end() const
+        {
+            return (*nzRow_)[0].end();
+        }
+
+    protected:
+        const SingleRowMatrix* nzRow_;
+        bool isNonZero_;
+    };
 
 public:
+    enum { rows = ROWS, cols = COLS};
 
-    typedef typename Dune::FieldMatrix<K, 1, Base::cols> SingleRowMatrix;
+    typedef RowProxy row_type;
+    typedef row_type const_row_reference;
+    typedef typename std::size_t size_type;
+    typedef typename RowProxy::ConstIterator ConstColIterator;
 
     /**
      * \brief Create from single row matrix and row index
      */
-    SingleNonZeroRowMatrix(const SingleRowMatrix& c, int row) :
-        Base(0.0),
-        row_(row)
+    SingleNonZeroRowMatrix(const SingleRowMatrix& r, size_type nzRowIndex) :
+        nzRow_(r),
+        nzRowIndex_(nzRowIndex)
+    {}
+
+    size_type N() const
+    {
+        return ROWS;
+    }
+
+    size_type M() const
     {
-        for(int i=0; i<Base::cols; ++i)
-            (*this)[row_][i] = c[0][i];
+        return COLS;
     }
 
     template<class X , class Y >
     void umv(const X& x, Y& y) const
     {
-        for(int i=0; i<Base::cols; ++i)
-            y[row_] += (*this)[row_][i] * x[i];
+        for(size_type i=0; i<M(); ++i)
+            y[nzRowIndex_] += nzRow_[0][i] * x[i];
     }
 
     template<class X , class Y >
     void umtv(const X& x, Y& y) const
     {
-        for(int i=0; i<Base::cols; ++i)
-            y[i] = (*this)[row_][i] * x[row_];
+        for(size_type i=0; i<N(); ++i)
+            y[i] = nzRow_[0][i] * x[nzRowIndex_];
     }
 
     template<class X , class Y >
@@ -59,19 +100,19 @@ public:
         umtv(x, y);
     }
 
-    int nonZeroRowIndex() const
+    const_row_reference operator[] (size_type rowIndex) const
     {
-        return row_;
+        return const_row_reference(&nzRow_, rowIndex, nzRowIndex_);
     }
 
-    const Base::row_type& nonZeroRow() const
+    size_type nonZeroRowIndex() const
     {
-        return (*this)[row_];
+        return nzRowIndex_;
     }
 
 protected:
-
-    const int row_;
+    SingleRowMatrix nzRow_;
+    const size_type nzRowIndex_;
 };
 
 
@@ -81,27 +122,14 @@ namespace Arithmetic
     template<class K, int ROWS, int COLS>
     struct MatrixTraits<SingleNonZeroRowMatrix<K, ROWS, COLS> >
     {
-        enum { isMatrix=true };
-        enum { rows=ROWS};
-        enum { cols=COLS};
-        typedef typename SingleNonZeroRowMatrix<K, ROWS, COLS>::row_type::ConstColIterator ConstColIterator;
-
-        static ConstColIterator rowBegin(const typename SingleNonZeroRowMatrix<K, ROWS, COLS>::row_type& m, int row)
-        {
-            if (row == m.nonZeroRowIndex())
-                return m.nonZeroRow().begin();
-            else
-                return m.nonZeroRow().end();
-        }
-
-        static ConstColIterator rowEnd(const typename SingleNonZeroRowMatrix<K, ROWS, COLS>::row_type& m, int row)
-        {
-            return m.nonZeroRow().end();
-        }
+        enum { isMatrix = true };
+        enum { isDense = false };
+        enum { sizeIsStatic = true };  //TODO only one of these should be used
+        enum { hasStaticSize = true }; //TODO only one of these should be used
+        enum { rows = ROWS};
+        enum { cols = COLS};
     };
-
-};
-
+}
 
 
 #endif