compressedmultigridtransfer.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef COMPRESSED_MULTIGRID_TRANSFER_HH
#define COMPRESSED_MULTIGRID_TRANSFER_HH

#include <memory>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/shared_ptr.hh>

#include <dune/matrix-vector/axpy.hh>
#include <dune/solvers/operators/sumoperator.hh>

#include "genericmultigridtransfer.hh"
#include "multigridtransfer.hh"


/** \brief Restriction and prolongation operator for standard multigrid
 *
 * This class implements the standard prolongation and restriction
 * operators for standard multigrid solvers.  Restriction and prolongation
 * of block vectors is provided.  Internally, the operator is stored
 * as a scalar BCRSMatrix.  Therefore, the template parameter VectorType
 * has to comply with the ISTL requirements.

 * \todo Currently only works for first-order Lagrangian elements!
 */
template<
    class VectorType,
    class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
    class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
        typename VectorType::field_type, VectorType::block_type::dimension, VectorType::block_type::dimension> > >
class CompressedMultigridTransfer :
    public virtual MultigridTransfer<VectorType,BitVectorType,MatrixType>
{

    enum {blocksize = VectorType::block_type::dimension};

    typedef typename VectorType::field_type field_type;

public:

    typedef MatrixType OperatorType;

    typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> > TransferOperatorType;

    CompressedMultigridTransfer()
    {
        matrix_ = std::make_shared<TransferOperatorType>();
    }

    CompressedMultigridTransfer(typename std::shared_ptr<TransferOperatorType>& matrix) :
        matrix_(matrix)
    {}

    virtual ~CompressedMultigridTransfer()
    {}


    /** \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, coarseMat, -1);
    }


    /** \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, coarseMat, -1);
    }


    /** \brief Direct access to the operator matrix, if you absolutely want it! */
    const TransferOperatorType& getMatrix() const {
        return *matrix_;
    }

    /** \brief Set matrix! */
    void setMatrix(typename std::shared_ptr<TransferOperatorType>& matrix)
    {
        matrix_ = matrix;
    }

    /** \brief Set matrix! */
    void setMatrix(TransferOperatorType& matrix)
    {
        matrix_ = Dune::stackobject_to_shared_ptr<TransferOperatorType>(matrix);
    }



protected:

    typename std::shared_ptr<TransferOperatorType> matrix_;
};

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_ = std::make_shared<TransferOperatorType>();
    }

    CompressedMultigridTransfer(typename std::shared_ptr<TransferOperatorType>& matrix) :
        matrix_(matrix)
    {}

    virtual ~CompressedMultigridTransfer()
    {}


    /** \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)
                    Dune::MatrixVector::addProduct(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(typename std::shared_ptr<TransferOperatorType>& matrix)
    {
        matrix_ = matrix;
    }

    /** \brief Set matrix! */
    void setMatrix(TransferOperatorType& matrix)
    {
        matrix_ = Dune::stackobject_to_shared_ptr<TransferOperatorType>(matrix);
    }



protected:

    typename std::shared_ptr<TransferOperatorType> matrix_;
};

#endif