diff --git a/dune/solvers/common/genericvectortools.hh b/dune/solvers/common/genericvectortools.hh index 30a9a954fb31fb3b33c559c4ba8ef9d9556a21a4..98b1b8f22b2874cec47338624a1c89df109cffd4 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