multigridstep.hh

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

#include <vector>
#include <map>
#include <dune/common/bitsetvector.hh>

#include <dune/solvers/transferoperators/multigridtransfer.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/common/wrapownshare.hh>

#include "lineariterationstep.hh"

namespace Dune {

  namespace Solvers {

/** \brief A linear multigrid step */
    template<
        class MatrixType,
        class VectorType,
        class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
    class MultigridStep : public LinearIterationStep<MatrixType, VectorType, BitVectorType>
    {

        static const int blocksize = VectorType::block_type::dimension;

    public:
        using LinearStepType = LinearIterationStep<MatrixType, VectorType, BitVectorType>;

        MultigridStep() :
            presmoother_(0),
            postsmoother_(0),
            basesolver_(0),
            ignoreNodesHierarchy_(0),
            preprocessCalled(false)
        {}

        MultigridStep(const MatrixType& mat,
                      VectorType& x,
                      const VectorType& rhs,
                      int mu, int nu1, int nu2,
                      LinearStepType* preSmoother,
                      LinearStepType* postSmoother,
                      Solver* baseSolver,
                      const BitVectorType* ignoreNodes) :
            LinearStepType(mat, x, rhs),
            preprocessCalled(false)
        {
            mu_  = mu;
            nu1_ = nu1;
            nu2_ = nu2;

            setSmoother(preSmoother,postSmoother);

            basesolver_   = baseSolver;

            this->ignoreNodes_ = ignoreNodes;
        }

        MultigridStep(const MatrixType& mat,
                      VectorType& x,
                      const VectorType& rhs) :
            LinearStepType(mat, x, rhs),
            basesolver_(0),
            preprocessCalled(false)
        {}

        virtual ~MultigridStep()
        {
            for (int i=0; i<int(ignoreNodesHierarchy_.size()-1); i++)
            {
                if (ignoreNodesHierarchy_[i])
                    delete(ignoreNodesHierarchy_[i]);
            }
        }

        virtual void setProblem(const MatrixType& mat,
                                VectorType& x,
                                const VectorType& rhs)
        {
            level_ = numLevels()-1;
            LinearIterationStep<MatrixType, VectorType, BitVectorType>::setProblem(mat,x,rhs);

            // Preprocess must be called again, to create the new matrix hierarchy
            preprocessCalled = false;
        }

        void setRhs(const VectorType& rhs)
        {
            level_ = numLevels()-1;
            this->rhs_ = &rhs;
            rhsHierarchy_.back() = rhs;
        }

        template <class DerivedTransferHierarchy>
        void setTransferOperators(const DerivedTransferHierarchy& transfer)
        {
            mgTransfer_.resize(transfer.size());
            for(size_t j=0; j<transfer.size(); ++j)
                mgTransfer_[j] = transfer[j];
        }

        /**
         * \brief Set transfer operator hierarchy from vector of shared_ptr's.
         *
         * Be careful: The Multigrid step does currently not share ownership
         * afterwards. This may change in the future.
         */
        template <class DerivedTransfer>
        void setTransferOperators(const std::vector<typename std::shared_ptr<DerivedTransfer> >& transfer)
        {
            mgTransfer_.resize(transfer.size());
            for(size_t j=0; j<transfer.size(); ++j)
                mgTransfer_[j] = transfer[j].get();
        }

        virtual void iterate();

        virtual void nestedIteration();

        virtual void preprocess();

        virtual void postprocess();

//        virtual const MatrixType* getMatrix();

        /** \brief Return total number of levels of the multigrid hierarchy */
        virtual size_t numLevels() const
        {
            return mgTransfer_.size() + 1;
        }

        virtual int level() const {return level_;}

        /** \brief Sets the number of pre- and postsmoothing steps
            and of coarse corrections.
            \param mu Number of coarse corrections
            \param nu1 Number of presmoothing steps
            \param nu2 Number of postsmoothing steps
        */
        virtual void setMGType(int mu, int nu1, int nu2);

        /** \brief Set the smoother iteration step */
        virtual void setSmoother(LinearStepType* smoother)
        {
                presmootherDefault_ = postsmootherDefault_ = Dune::stackobject_to_shared_ptr(*smoother);

                levelWiseSmoothers_.clear();
        }

        /** \brief Set the smoother iteration step from a smart pointer*/
        virtual void setSmoother(std::shared_ptr<LinearStepType> smoother)
        {
                presmootherDefault_ = postsmootherDefault_ = smoother;

                levelWiseSmoothers_.clear();
        }

        /** \brief Set pre- and post smoothers individually */
        virtual void setSmoother(LinearStepType* preSmoother,
                                 LinearStepType* postSmoother)
        {
                presmootherDefault_  = Dune::stackobject_to_shared_ptr(*preSmoother);
                postsmootherDefault_ = Dune::stackobject_to_shared_ptr(*postSmoother);

                levelWiseSmoothers_.clear();
        }

        /** \brief Set the smoother iteration step for a particular level */
        virtual void setSmoother(LinearStepType* smoother, std::size_t level)
        {
            levelWiseSmoothers_[level] = Dune::stackobject_to_shared_ptr(*smoother);
        }

        /** \brief Set the smoother iteration step for a particular level, from a smart pointer */
        virtual void setSmoother(std::shared_ptr<LinearStepType> smoother, std::size_t level)
        {
            levelWiseSmoothers_[level] = smoother;
        }

        /** \brief Set base solver */
        template <class BaseSolver>
        void setBaseSolver(BaseSolver&& baseSolver)
        {
            basesolver_ = wrap_own_share<Solver>(std::forward<BaseSolver>(baseSolver));
        }

    protected:
        /** \brief The presmoothers, one for each level */
        std::vector<std::shared_ptr<LinearStepType> > presmoother_;

        /** \brief The postsmoothers, one for each level */
        std::vector<std::shared_ptr<LinearStepType> > postsmoother_;

        /** \brief The base solver */
        std::shared_ptr<Solver> basesolver_;

        //! Number of presmoothing steps
        int nu1_;

        //! Number of postsmoothing steps
        int nu2_;

        //! Number of coarse corrections
        int mu_;
    public:
        //! Variable used to store the current level
        int level_;

        //! The linear operators on each level
        std::vector<std::shared_ptr<const MatrixType> > matrixHierarchy_;

    protected:
        //! Flags specifying the dirichlet nodes on each level
        std::vector<BitVectorType*> ignoreNodesHierarchy_;

    public:
        /**  \brief! hierarchy of solution/correction vectors.
          *
          *  on the fine level it contains the iterate, whereas on the coarse levels the corresponding corrections
          */
        std::vector<std::shared_ptr<VectorType> > xHierarchy_;

        std::vector<VectorType> rhsHierarchy_;

        //protected:
        std::vector<MultigridTransfer<VectorType, BitVectorType, MatrixType>* > mgTransfer_;

    protected:

        std::shared_ptr<LinearStepType> presmootherDefault_;
        std::shared_ptr<LinearStepType> postsmootherDefault_;
        typedef std::map<std::size_t, std::shared_ptr<LinearStepType> > SmootherCache;
        SmootherCache levelWiseSmoothers_;

        bool preprocessCalled;
    };

  }   // namespace Solvers

}   // namespace Dune


// For backward compatibility: will be removed eventually
using Dune::Solvers::MultigridStep;

#include "multigridstep.cc"

#endif