densemultigridtransfer.hh

#ifndef DENSE_MULTIGRID_TRANSFER_HH
#define DENSE_MULTIGRID_TRANSFER_HH

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

#include "multigridtransfer.hh"
#include "genericmultigridtransfer.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 BCRSMatrix of small dense matrices.
 * 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 DenseMultigridTransfer :
    public MultigridTransfer<VectorType,BitVectorType,MatrixType>
{

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

    typedef typename VectorType::field_type field_type;

public:

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

    virtual ~DenseMultigridTransfer() {}


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


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


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


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


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


    /** \brief Galerkin assemble a coarse stiffness matrix
     */
    virtual void galerkinRestrict(const MatrixType& fineMat, MatrixType& coarseMat) const
    {
        GenericMultigridTransfer::galerkinRestrict(matrix_, fineMat, coarseMat);
    }


    /** \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 MatrixType& fineMat, MatrixType& coarseMat) const
    {
        GenericMultigridTransfer::galerkinRestrictSetOccupation(matrix_, fineMat, coarseMat);
    }


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


protected:

    TransferOperatorType matrix_;

};

#endif