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