gsMultiGrid.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsCore/gsForwardDeclarations.h>
#include <gsSolver/gsPreconditioner.h>
#include <gsIO/gsOptionList.h>
 
namespace gismo
{
 
template<class T>
class gsMultiGridOp : public gsPreconditionerOp<T>
{
 
public:
 
    typedef memory::shared_ptr<gsMultiGridOp> Ptr;
 
    typedef memory::unique_ptr<gsMultiGridOp> uPtr;
 
    typedef typename gsLinearOperator<T>::Ptr OpPtr;
 
    typedef typename gsPreconditionerOp<T>::Ptr PrecondPtr;
 
    typedef gsPreconditionerOp<T> Base;
 
    typedef gsMatrix<T> Matrix;
 
    typedef gsSparseMatrix<T> SpMatrix;
 
    typedef gsSparseMatrix<T, RowMajor> SpMatrixRowMajor;
 
    typedef memory::shared_ptr<SpMatrix> SpMatrixPtr;
 
    typedef memory::shared_ptr<SpMatrixRowMajor> SpMatrixRowMajorPtr;
 
    gsMultiGridOp(
        SpMatrix fineMatrix,
        std::vector<SpMatrixRowMajor> transferMatrices,
        OpPtr coarseSolver = OpPtr()
    );
 
    gsMultiGridOp(
        SpMatrixPtr fineMatrix,
        std::vector<SpMatrixRowMajorPtr> transferMatrices,
        OpPtr coarseSolver = OpPtr()
    );
 
    gsMultiGridOp(
        const std::vector<OpPtr>& ops,
        const std::vector<OpPtr>& prolongation,
        const std::vector<OpPtr>& restriction,
        OpPtr coarseSolver = OpPtr()
    );
 
    static uPtr make(
        SpMatrix fineMatrix,
        std::vector<SpMatrixRowMajor> transferMatrices,
        OpPtr coarseSolver = OpPtr()
    )
    { return uPtr( new gsMultiGridOp( give(fineMatrix), give(transferMatrices), give(coarseSolver) ) ); }
 
    static uPtr make(
        SpMatrixPtr fineMatrix,
        std::vector<SpMatrixRowMajorPtr> transferMatrices,
        OpPtr coarseSolver = OpPtr()
    )
    { return uPtr( new gsMultiGridOp( give(fineMatrix), give(transferMatrices), give(coarseSolver) ) ); }
 
    static uPtr make(
        const std::vector<OpPtr>& ops,
        const std::vector<OpPtr>& prolongation,
        const std::vector<OpPtr>& restriction,
        OpPtr coarseSolver = OpPtr()
    )
    { return uPtr( new gsMultiGridOp( ops, prolongation, restriction, give(coarseSolver) ) ); }
 
private:
    void init(SpMatrixPtr fineMatrix, std::vector<SpMatrixRowMajorPtr> transferMatrices);
    void initCoarseSolver() const;
public:
 
    void smoothingStep(index_t level, const Matrix& rhs, Matrix& x) const;
 
    void smoothingStep(const Matrix& rhs, Matrix& x) const
    { smoothingStep(finestLevel(), rhs, x); }
 
    void multiGridStep(index_t level, const Matrix& rhs, Matrix& x) const;
 
    void step(const Matrix& rhs, Matrix& x) const override
    { multiGridStep(finestLevel(), rhs, x); }
 
    void stepT(const Matrix& rhs, Matrix& x) const override
    {
        std::swap(m_numPreSmooth,m_numPostSmooth);
        step(rhs,x);
        std::swap(m_numPreSmooth,m_numPostSmooth);
    }
 
    void fullMultiGrid(
        const std::vector<Matrix>& rhs,
        const std::vector<Matrix>& dirichletIntp,
        Matrix& result
    ) const;
 
    void cascadicMultiGrid(
        const std::vector<Matrix>& rhs,
        const std::vector<Matrix>& fixedValues,
        Matrix& result
    ) const;
 
    void restrictVector(index_t lf, const Matrix& fine, Matrix& coarse) const;
 
    void prolongVector(index_t lc, const Matrix& coarse, Matrix& fine) const;
 
    void solveCoarse(const Matrix& rhs, Matrix& result) const
    {
        if (!m_coarseSolver) initCoarseSolver();
        m_coarseSolver->apply(rhs, result);
    }
 
    index_t numLevels() const                         { return m_nLevels;                 }
 
    index_t finestLevel() const                       { return m_nLevels - 1;             }
 
    void setUnderlyingOp(index_t lvl, OpPtr op)
    {
        GISMO_ASSERT ( lvl >= 0 && lvl < m_nLevels, "gsMultiGridOp: The given level "<<lvl<<" is not feasible." );
        m_ops[lvl] = op;
    }
 
    OpPtr underlyingOp(index_t lvl) const             { return m_ops[lvl];               }
 
    OpPtr underlyingOp() const override               { return m_ops[finestLevel()];     }
 
    const SpMatrix& matrix(index_t lvl) const;
 
    const SpMatrix& matrix() const                    { return matrix(finestLevel());     }
 
    index_t nDofs(index_t lvl) const                  { return underlyingOp(lvl)->cols(); }
 
    index_t nDofs() const                             { return nDofs( finestLevel() );    }
 
    // For docs, see base class
    index_t rows() const override                     { return underlyingOp()->rows();    }
 
    // For docs, see base class
    index_t cols() const override                     { return underlyingOp()->cols();    }
 
    void setSmoother(index_t lvl, const PrecondPtr& sm);
 
    PrecondPtr smoother(index_t lvl) const            { return m_smoother[lvl];           }
 
    void setCoarseSolver(const OpPtr& sol)            { m_coarseSolver = sol;             }
 
    OpPtr coarseSolver() const                        { return m_coarseSolver;            }
 
    void setNumPreSmooth(index_t n)                   { m_numPreSmooth = n;               }
 
    index_t setNumPreSmooth() const                   { return m_numPreSmooth;            }
 
    void setNumPostSmooth(index_t n)                  { m_numPostSmooth = n;              }
 
    index_t setNumPostSmooth() const                  { return m_numPostSmooth;           }
 
    void setSymmSmooth(bool s)                        { m_symmSmooth = s;                 }
 
    bool symmSmooth() const                           { return m_symmSmooth;              }
 
    void setNumCycles(index_t n)
    {
        m_numCycles.setConstant(m_nLevels-1,n);
        // The direct solver on coarsest level is only invoked once
        m_numCycles[0] = 1;
    }
 
    void setNumCycles(index_t lvl, index_t n)
    {
        GISMO_ASSERT(lvl>=0&&lvl<m_nLevels-1, "gsMultiGrid::setNumCycles: Givel level out of bounds.");
        m_numCycles[lvl] = n;
    }
 
    index_t numCycles(index_t lvl) const              { return m_numCycles[lvl];           }
 
    void setCoarseGridCorrectionDamping(T damping)    { m_damping = damping;               }
 
    T coarseGridCorrectionDamping() const             { return m_damping;                  }
 
    static gsOptionList defaultOptions();
 
    void setOptions(const gsOptionList& opt) override;
 
private:
 
    index_t m_nLevels;
 
    std::vector<OpPtr> m_ops;
 
    std::vector<PrecondPtr> m_smoother;
 
    std::vector<OpPtr> m_prolong;
 
    std::vector<OpPtr> m_restrict;
 
    mutable OpPtr m_coarseSolver;
 
    mutable index_t m_numPreSmooth;
 
    mutable index_t m_numPostSmooth;
 
    bool m_symmSmooth;
 
    gsVector<index_t> m_numCycles;
 
    T m_damping;
 
}; // class gsMultiGridOp
 
}  // namespace gismo
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsMultiGrid.hpp)
#endif