From b2abd51d24ab0c8794e3bfacd719656d108b2b08 Mon Sep 17 00:00:00 2001
From: Uli Sack <usack@math.fu-berlin.de>
Date: Mon, 23 Jul 2012 13:53:13 +0000
Subject: [PATCH] introduces two new methods interlace and deinterlace to
 'interlace' resp. 'deinterlace' blockvectorsi. Implemented so far as
 interlace: \f$\left((R^k)^m\right)^n\longrightarrow(R^k\cdot n)^m\f$, i.e.
 for example       / / /a11\ \ \               / /a11\ \      | |  \a12/  | | 
            | | a12 | |      | |  /a21\  | |             | | b11 | |      | | 
 \a22/  | |             | | b12 | |      | |  /a31\  | |             | | c11 |
 |      |  \ \a32/ /  |             | | c12 | |      |             |          
   | | d11 | |      |  / /b11\ \  |             |  \d12/  |      | |  \b12/  |
 |             |         |      | |  /b21\  | |             |  /a21\  |      |
 |  \b22/  | |             | | a22 | |      | |  /b31\  | |             | |
 b21 | |      |  \ \b32/ /  |  interlace  | | b22 | |      |             |
 ----------> | | c11 | |      |  / /c11\ \  |             | | c22 | |      | |
  \c12/  | |             | | d21 | |      | |  /c21\  | |             |  \d22/
  |      | |  \c22/  | |             |         |      | |  /c31\  | |         
    |  /a31\  |      |  \ \c32/ /  |             | | a32 | |      |           
  |             | | b31 | |      |  / /d11\ \  |             | | b32 | |     
 | |  \d12/  | |             | | c31 | |      | |  /d21\  | |             | |
 c32 | |      | |  \d22/  | |             | | d31 | |      | |  /d31\  | |    
         |  \d32/  |       \ \ \d32/ / /               \       /

[[Imported from SVN: r6541]]
---
 dune/solvers/common/genericvectortools.hh | 50 +++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/dune/solvers/common/genericvectortools.hh b/dune/solvers/common/genericvectortools.hh
index 30a9a954..98b1b8f2 100644
--- a/dune/solvers/common/genericvectortools.hh
+++ b/dune/solvers/common/genericvectortools.hh
@@ -9,8 +9,10 @@
 #include<bitset>
 
 #include <dune/common/fvector.hh>
+#include <dune/common/static_assert.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
 
 #include <dune/solvers/operators/sumoperator.hh>
 
@@ -124,6 +126,54 @@ struct GenericVector
                 *it = 0;
     }
 
+    //! weave two vector blocks into each other e.g. [[u1...un][w1...wn]] --> [[u1 w1]...[un wn]]
+    template <class LVectorType, class RVectorType>
+    static void interlace(const LVectorType& lvec, RVectorType& rvec)
+    {
+        DUNE_THROW(Dune::NotImplemented,"GenericVector::interlace not implemented for given VectorTypes");
+    }
+
+    template <class LVectorBlock, class RVectorBlock>
+    static void interlace(const Dune::BlockVector<Dune::BlockVector<LVectorBlock> >& lvec, Dune::BlockVector<RVectorBlock>& rvec)
+    {
+        dune_static_assert(RVectorBlock::dimension % LVectorBlock::dimension == 0, "Block sizes don't match.");
+
+        rvec.resize(lvec[0].size());
+        rvec = 0.0;
+
+        size_t offset = LVectorBlock::dimension;
+
+        for (size_t i=0; i<lvec[0].size(); ++i)
+            for (size_t k=0; k<lvec.size(); ++k)
+                for (size_t j=0; j<offset; ++j)
+                    rvec[i][j+k*offset] = lvec[k][i][j];
+    }
+
+    //! unweave two vectors previously interlaced e.g. [[u1 w1]...[un wn]] --> [[u1...un][w1...wn]]
+    template <class LVectorType, class RVectorType>
+    static void deinterlace(const LVectorType& lvec, RVectorType& rvec)
+    {
+        DUNE_THROW(Dune::NotImplemented,"GenericVector::deinterlace not implemented for given VectorTypes");
+    }
+
+    template <class LVectorBlock, class RVectorBlock>
+    static void deinterlace(const Dune::BlockVector<LVectorBlock>& lvec, Dune::BlockVector<Dune::BlockVector<RVectorBlock> >& rvec)
+    {
+        dune_static_assert(LVectorBlock::dimension % RVectorBlock::dimension == 0, "Block sizes don't match.");
+
+        rvec.resize(LVectorBlock::dimension / RVectorBlock::dimension);
+        for (size_t k=0; k<rvec.size(); ++k)
+            rvec[k].resize(lvec.size());
+        rvec = 0.0;
+
+        size_t offset = RVectorBlock::dimension;
+
+        for (size_t i=0; i<lvec.size(); ++i)
+            for (size_t k=0; k<rvec.size(); ++k)
+                for (size_t j=0; j<offset; ++j)
+                    rvec[k][i][j] = lvec[i][j+k*offset];
+    }
+
 };
 
 #endif
-- 
GitLab