added specialization for MatrixType=SumOperator<SparseMatrix, LowRankMatrix>

[[Imported from SVN: r5474]]

dune/solvers/transferoperators/compressedmultigridtransfer.hh

@@ -5,6 +5,9 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/bitsetvector.hh>
 
+#include <dune/solvers/operators/sumoperator.hh>
+#include <dune/solvers/common/staticmatrixtools.hh>
+
 #include "genericmultigridtransfer.hh"
 #include "multigridtransfer.hh"
 
@@ -142,6 +145,157 @@ public:
 
 
 
+protected:
+
+    TransferOperatorType* matrix_;
+    bool matrixInternallyAllocated;
+
+};
+
+template<
+    class VectorType,
+    class BitVectorType,
+    class SparseMatrixType,
+    class LowRankMatrixType >
+class CompressedMultigridTransfer<VectorType, BitVectorType, SumOperator<SparseMatrixType, LowRankMatrixType> > :
+    public MultigridTransfer<VectorType, BitVectorType, SumOperator<SparseMatrixType,LowRankMatrixType> >
+{
+
+    enum {blocksize = VectorType::block_type::dimension};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+
+    typedef SumOperator<SparseMatrixType, LowRankMatrixType> OperatorType;
+
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> > TransferOperatorType;
+
+    CompressedMultigridTransfer()
+    {
+        matrix_ = new TransferOperatorType;
+        matrixInternallyAllocated = true;
+    }
+
+    virtual ~CompressedMultigridTransfer()
+    {
+        if (matrixInternallyAllocated)
+            delete matrix_;
+    }
+
+
+    /** \brief Sets up the operator between two P1 spaces
+     *
+     * \param grid The grid
+     * \param cL The coarse grid level
+     * \param fL The fine grid level
+     */
+    template <class GridType>
+    void setup(const GridType& grid, int cL, int fL)
+    {
+        GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(*matrix_, grid, cL, fL);
+    }
+
+
+    /** \brief Restrict a function from the fine onto the coarse grid
+     */
+    virtual void restrict(const VectorType& f, VectorType &t) const
+    {
+        GenericMultigridTransfer::restrict(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Restrict a bitfield from the fine onto the coarse grid
+     */
+    virtual void restrictScalarBitField(const Dune::BitSetVector<1>& f, Dune::BitSetVector<1>& t) const
+    {
+        GenericMultigridTransfer::restrictScalarBitField(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     */
+    virtual void restrict(const BitVectorType& f, BitVectorType& t) const
+    {
+        GenericMultigridTransfer::restrictBitField(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Restrict a vector valued bitfield from the fine onto the coarse grid
+     * Fine bits only influence their father bits
+     */
+    virtual void restrictToFathers(const BitVectorType& f, BitVectorType& t) const
+    {
+        GenericMultigridTransfer::restrictBitFieldToFathers(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Prolong a function from the coarse onto the fine grid
+     */
+    virtual void prolong(const VectorType& f, VectorType &t) const
+    {
+        GenericMultigridTransfer::prolong(*matrix_, f, t, -1);
+    }
+
+
+    /** \brief Galerkin assemble a coarse stiffness matrix
+     */
+    virtual void galerkinRestrict(const OperatorType& fineMat, OperatorType& coarseMat) const
+    {
+        GenericMultigridTransfer::galerkinRestrict(*matrix_, fineMat.sparseMatrix(), coarseMat.sparseMatrix(), -1);
+
+        const typename LowRankMatrixType::LowRankFactorType& lr_fine_factor(fineMat.lowRankMatrix().lowRankFactor());
+        typename LowRankMatrixType::LowRankFactorType& lr_coarse_factor(coarseMat.lowRankMatrix().lowRankFactor());
+
+        lr_coarse_factor = 0;
+
+        for (size_t i=0; i<matrix_->N(); ++i)
+        {
+            typename TransferOperatorType::ColIterator j_it  = (*matrix_)[i].begin();
+            typename TransferOperatorType::ColIterator j_end = (*matrix_)[i].end();
+
+            for (; j_it!=j_end; ++j_it)
+            {
+                size_t j = j_it.index();
+
+                for (size_t k=0; k<lr_fine_factor.N(); ++k)
+                    StaticMatrix::axpy(lr_coarse_factor[k][j], *j_it, lr_fine_factor[k][i]);
+            }
+        }
+    }
+
+
+    /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix   
+    *
+    * Set occupation of Galerkin restricted coarse matrix. Call this one before
+    * galerkinRestrict to ensure all non-zeroes are present
+    * \param fineMat The fine level matrix
+    * \param coarseMat The coarse level matrix
+    */
+    virtual void galerkinRestrictSetOccupation(const OperatorType& fineMat, OperatorType& coarseMat) const
+    {
+        GenericMultigridTransfer::galerkinRestrictSetOccupation(*matrix_, fineMat.sparseMatrix(), coarseMat.sparseMatrix(), -1);
+
+        coarseMat.lowRankMatrix().lowRankFactor().setSize(fineMat.lowRankMatrix().lowRankFactor().N(), matrix_->M());
+    }
+
+
+    /** \brief Direct access to the operator matrix, if you absolutely want it! */
+    const TransferOperatorType& getMatrix() const {
+        return *matrix_;
+    }
+
+    /** \brief Set matrix! */
+    void setMatrix(TransferOperatorType& matrix)
+    {
+        if (matrixInternallyAllocated)
+            delete matrix_;
+        matrixInternallyAllocated = false;
+        matrix_ = &matrix;
+    }
+
+
+
 protected:
 
     TransferOperatorType* matrix_;