gsINSPrecondBlocks.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsIncompressibleFlow/src/gsFlowUtils.h>
 
#include <gsSolver/gsLinearOperator.h>
#include <gsSolver/gsMatrixOp.h>
 
namespace gismo
{
 
// === BASE BLOCKS ==== //
 
template <class T>
class gsINSPrecondBlock : public gsLinearOperator<T>
{
 
public:
    typedef gsLinearOperator<T> Base;
 
protected: // *** Class members ***
 
    gsSparseSolver<T>* m_pSolver;
    gsOptionList m_opt;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondBlock> Ptr; 
    typedef memory::unique_ptr<gsINSPrecondBlock> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondBlock()
    {
        m_pSolver = NULL;
    }
 
    gsINSPrecondBlock(const gsOptionList& opt) : gsINSPrecondBlock()
    {
        m_opt = opt;
 
        m_pSolver = createLinSolver();
        setupLinSolver(*m_pSolver);
    }
 
    virtual ~gsINSPrecondBlock()
    {
        if (m_pSolver)
        {
            delete m_pSolver;
            m_pSolver = NULL;
        }
    }
 
public: // *** Member functions ***
 
    gsSparseSolver<T>* createLinSolver();
 
    void setupLinSolver(gsSparseSolver<T>& solver);
 
    virtual std::string getName() const
    { GISMO_NO_IMPLEMENTATION }
 
}; // class gsINSPrecondBlock
 
 
template <class T>
class gsINSPrecondBlockMod : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    std::vector<gsSparseSolver<T>* > m_solvers;
 
protected: // *** Base class members ***
 
    using Base::m_opt;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondBlockMod> Ptr;
    typedef memory::unique_ptr<gsINSPrecondBlockMod> uPtr;
 
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondBlockMod(const gsOptionList& opt) : Base()
    {
        m_opt = opt;
        int dim = m_opt.getInt("dim");
 
        m_solvers.resize(dim);
        
        for (int i = 0; i < dim; i++)
        {
            m_solvers[i] = this->createLinSolver();
            this->setupLinSolver(*m_solvers[i]);
        }
    }
 
    virtual ~gsINSPrecondBlockMod()
    {
        for (int i = 0; i < (int)m_solvers.size(); i++)
        {
            if (m_solvers[i])
            {
                delete m_solvers[i];
                m_solvers[i] = NULL;
            }
        }
    }
 
}; // class gsINSPrecondBlockMod
 
 
// === BLOCK F ==== //
 
template <class T, int MatOrder>
class gsINSPrecondBlockF : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    const gsSparseMatrix<T, MatOrder> m_blockF1;
    int m_size;
 
protected: // *** Base class members ***
 
    using Base::m_pSolver;
    using Base::m_opt;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondBlockF> Ptr;
    typedef memory::unique_ptr<gsINSPrecondBlockF> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondBlockF(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt) :
        Base(opt), m_blockF1(matNS.block(0, 0, opt.getInt("udofs"), opt.getInt("udofs")))
    {
        m_size = opt.getInt("dim") * opt.getInt("udofs");
 
        m_pSolver->compute(m_blockF1);
    }
    
public: // *** Static functions ***
 
    static uPtr make(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondBlockF<T, MatOrder>(matNS, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const;
 
 
public: // *** Getters/setters ***
 
    int rows() const { return m_size; }
 
    int cols() const { return m_size; }
 
    virtual std::string getName() const
    {
        return "FdiagEqual";
    }
 
}; // class gsINSPrecondBlockF
 
 
// === BLOCK Fwhole ==== //
 
template <class T, int MatOrder>
class gsINSPrecondBlockFwhole : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    const gsSparseMatrix<T, MatOrder> m_blockA;
    int m_size;
 
protected: // *** Base class members ***
 
    using Base::m_pSolver;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondBlockFwhole> Ptr;
    typedef memory::unique_ptr<gsINSPrecondBlockFwhole> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondBlockFwhole(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt) :
        Base(opt), m_blockA(matNS.block(0, 0, opt.getInt("dim")*opt.getInt("udofs"), opt.getInt("dim")*opt.getInt("udofs")))
    {
        m_size = opt.getInt("dim") * opt.getInt("udofs");
 
        m_pSolver->compute(m_blockA);
    }
 
public: // *** Static functions ***
 
    static uPtr make(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondBlockFwhole<T, MatOrder>(matNS, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const;
 
public: // *** Getters/setters ***
 
    int rows() const { return m_size; }
 
    int cols() const { return m_size; }
 
    virtual std::string getName() const
    {
        return "Fwhole";
    }
 
}; // class gsINSPrecondBlockFwhole
 
 
// === BLOCK Fdiag ==== //
 
template <class T, int MatOrder>
class gsINSPrecondBlockFdiag : public gsINSPrecondBlockMod<T>
{
 
public:
    typedef gsINSPrecondBlockMod<T> Base;
 
protected: // *** Class members ***
    std::vector<gsSparseMatrix<T, MatOrder> > m_diagBlocks;
    int m_size;
 
protected: // *** Base class members ***
 
    using Base::m_solvers;
    using Base::Base::m_opt;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondBlockFdiag> Ptr;
    typedef memory::unique_ptr<gsINSPrecondBlockFdiag> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondBlockFdiag(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt) :
        Base(opt)
    {
        int dim = opt.getInt("dim");
        m_size = dim * opt.getInt("udofs");
 
        m_diagBlocks.resize(dim);
 
        for (int i = 0; i < dim; i++)
        {
            int udofs = opt.getInt("udofs");
            m_diagBlocks[i] = matNS.block(i*udofs, i*udofs, udofs, udofs);
            m_solvers[i]->compute(m_diagBlocks[i]);
        }
    }
 
public: // *** Static functions ***
 
    static uPtr make(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondBlockFdiag<T, MatOrder>(matNS, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const;
 
public: // *** Getters/setters ***
 
    int rows() const { return m_size; }
 
    int cols() const { return m_size; }
 
    virtual std::string getName() const
    {
        return "Fdiag";
    }
 
}; // class gsINSPrecondBlockFdiag
 
 
// === BLOCK Fmod ==== //
 
template <class T, int MatOrder>
class gsINSPrecondBlockFmod : public gsINSPrecondBlockMod<T>
{
 
public:
    typedef gsINSPrecondBlockMod<T> Base;
 
protected: // *** Class members ***
 
    const gsSparseMatrix<T, MatOrder>& m_matRef;
    std::vector<gsSparseMatrix<T, MatOrder> > m_diagBlocks;
    int m_size;
 
protected: // *** Base class members ***
 
    using Base::m_solvers;
    using Base::Base::m_opt;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondBlockFmod> Ptr;
    typedef memory::unique_ptr<gsINSPrecondBlockFmod> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondBlockFmod(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt) :
        Base(opt), m_matRef(matNS)
    {
        int dim = opt.getInt("dim");
        int udofs = opt.getInt("udofs");
        m_size = dim * udofs;
 
       m_diagBlocks.resize(dim);
 
        for (int i = 0; i < dim; i++)
        {
            m_diagBlocks[i] = m_matRef.block(i*udofs, i*udofs, udofs, udofs);
            m_solvers[i]->compute(m_diagBlocks[i]);
        }
    }
 
public: // *** Static functions ***
 
    static uPtr make(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondBlockFmod<T, MatOrder>(matNS, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const;
 
public: // *** Getters/setters ***
 
    int rows() const { return m_size; }
 
    int cols() const { return m_size; }
 
    virtual std::string getName() const
    {
        return "Fmod";
    }
 
}; // class gsINSPrecondBlockFmod
 
 
// === BLOCK Bt ==== //
 
template <class T, int MatOrder>
class gsINSPrecondBlockBt : public gsLinearOperator<T>
{
 
public:
    typedef gsLinearOperator<T> Base;
 
protected: // *** Class members ***
 
    gsSparseMatrix<T, MatOrder> m_mat;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondBlockBt> Ptr;
    typedef memory::unique_ptr<gsINSPrecondBlockBt> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondBlockBt(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt)
    {
        int dim = opt.getInt("dim");
        int udofs = opt.getInt("udofs");
        m_mat = matNS.block(0, dim * udofs, dim * udofs, matNS.cols() - dim * udofs);
    }
 
public: // *** Static functions ***
 
    static uPtr make(const gsSparseMatrix<T, MatOrder>& matNS, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondBlockBt<T, MatOrder>(matNS, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const;
 
public: // *** Getters/setters ***
 
    int rows() const { return m_mat.rows(); }
 
    int cols() const { return m_mat.cols(); }
 
}; // class gsINSPrecondBlockBt
 
 
// === Schur LSC ==== //
 
template <class T, int MatOrder>
class gsINSPrecondSchurLSC : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    gsSparseMatrix<T, MatOrder> m_mat1, m_mat2;
 
protected: // *** Base class members ***
 
    using Base::m_pSolver;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondSchurLSC> Ptr;
    typedef memory::unique_ptr<gsINSPrecondSchurLSC> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondSchurLSC(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt) :
        Base(opt)
    {
        // S = - (B * velM^-1 * BT) * (B * velM^-1 * F * velM^-1 * BT)^-1 * (B * velM^-1 * BT)
        // denote S = - mat1 * mat2^-1 * mat1
 
        int dim = opt.getInt("dim");
        int udofs = opt.getInt("udofs");
        int uSize = dim * udofs; // number of all velocity dofs (all components)
        int pdofs = opt.getInt("pdofs"); // number of pressure dofs
 
        const gsSparseMatrix<T, MatOrder>& matNS = mat.at("matNS");
        const gsSparseMatrix<T, MatOrder>& velM = mat.at("matMu");
 
        gsSparseMatrix<T, MatOrder> velMinv(uSize, uSize); // approximation of velocity mass matrix inverse
 
        diagInvMatrix_into(velM, velMinv, uSize / velM.rows(), opt.getSwitch("lumpingM"));
        
        // TODO: could be more efficient if the matrices not involving block F were not updated (they do not change)
        gsSparseMatrix<T, MatOrder> BMinv, MinvBt;
        BMinv = matNS.block(uSize, 0, pdofs, uSize) * velMinv; // B * velM^-1
        MinvBt = velMinv * matNS.block(0, uSize, uSize, pdofs); // velM^-1 * BT
 
        m_mat1 = BMinv * matNS.block(0, uSize, uSize, pdofs); // B * velM^-1 * BT
        m_mat2 = BMinv * matNS.block(0, 0, uSize, uSize) * MinvBt; // B * velM^-1 * F * velM^-1 * BT
 
        m_pSolver->compute(m_mat1);
    }
 
public: // *** Static functions ***
 
    static uPtr make(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondSchurLSC<T, MatOrder>(mat, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T>& input, gsMatrix<T>& x) const;
 
public: // *** Getters/setters ***
 
    int rows() const { return m_mat2.rows(); }
 
    int cols() const { return m_mat2.cols(); }
 
}; // class gsINSPrecondSchurLSC
 
 
// === Schur PCD ==== //
 
template <class T, int MatOrder>
class gsINSPrecondSchurPCD : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    gsSparseMatrix<T, MatOrder> m_matAp, m_matFp, m_matMp;
    gsSparseSolver<T>* m_pMassSolver;
 
protected: // *** Base class members ***
 
    using Base::m_pSolver;
    
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondSchurPCD> Ptr; 
    typedef memory::unique_ptr<gsINSPrecondSchurPCD> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondSchurPCD(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt) :
        Base(opt)
    {
        // S = - Ap * Fp^-1 * Mp
 
        m_matAp = mat.at("matAp");
        m_matFp = mat.at("matFp");
        m_matMp = mat.at("matMp");
        
        m_pMassSolver = this->createLinSolver();
        this->setupLinSolver(*m_pMassSolver);
 
        m_pSolver->compute(m_matAp);
        m_pMassSolver->compute(m_matMp);
    }
 
    ~gsINSPrecondSchurPCD()
    {
        if (m_pMassSolver)
        {
            delete m_pMassSolver;
            m_pMassSolver = NULL;
        }
    }
 
public: // *** Static functions ***
 
    static uPtr make(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondSchurPCD<T, MatOrder>(mat, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T>& input, gsMatrix<T>& x) const;
 
public: // *** Getters/setters ***
 
    int rows() const { return m_matAp.rows(); }
 
    int cols() const { return m_matAp.cols(); }
 
}; // class gsINSPrecondSchurPCD
 
 
// === Schur PCDmod ==== // 
 
template <class T, int MatOrder>
class gsINSPrecondSchurPCDmod : public gsINSPrecondSchurPCD<T, MatOrder>
{
 
public:
    typedef gsINSPrecondSchurPCD<T, MatOrder> Base;
 
protected: // *** Base class members ***
 
    using Base::m_matAp;
    using Base::m_matFp;
    using Base::m_matMp;
    using Base::m_pMassSolver;
    using Base::Base::m_pSolver;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondSchurPCDmod> Ptr;
    typedef memory::unique_ptr<gsINSPrecondSchurPCDmod> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondSchurPCDmod(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt) :
        Base(mat, opt)
    {
        // S = - Mp * Fp^-1 * Ap
    }
 
public: // *** Static functions ***
 
    static uPtr make(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondSchurPCDmod<T, MatOrder>(mat, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T>& input, gsMatrix<T>& x) const;
 
}; // class gsINSPrecondSchurPCDmod
 
 
// === Schur AL ==== // 
 
template <class T, int MatOrder>
class gsINSPrecondSchurAL : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    real_t m_const;
    gsSparseMatrix<T, MatOrder> m_mat;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondSchurAL> Ptr;  
    typedef memory::unique_ptr<gsINSPrecondSchurAL> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondSchurAL(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        // S = - (gamma + nu) Mp
 
        const gsSparseMatrix<T, MatOrder>& presM = mat.at("matMp");
 
        m_const = opt.getReal("visc") + opt.getReal("gamma");
 
        diagInvMatrix_into(presM, m_mat, 1, opt.getSwitch("lumpingM"));
    }
 
public: // *** Static functions ***
 
    static uPtr make(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondSchurAL<T, MatOrder>(mat, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        GISMO_ASSERT(input.rows() == this->rows(), "Wrong input size.");
        x.resize(this->rows(), 1);
 
        x = - m_const * m_mat * input;
    }
 
public: // *** Getters/setters ***
 
    int rows() const { return m_mat.rows(); }
 
    int cols() const { return m_mat.cols(); }
 
}; // class gsINSPrecondSchurAL
 
 
// === Schur SIMPLE ==== // 
 
template <class T, int MatOrder>
class gsINSPrecondSchurSIMPLE : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    gsSparseMatrix<T, MatOrder> m_mat;
 
protected: // *** Base class members ***
 
    using Base::m_pSolver;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondSchurSIMPLE> Ptr;  
    typedef memory::unique_ptr<gsINSPrecondSchurSIMPLE> uPtr;
 
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondSchurSIMPLE(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt) :
        Base(opt)
    {
        // S = - B Dinv Bt
        int dim = opt.getInt("dim");
        int udofs = opt.getInt("udofs");
        int uSize = dim * udofs;
        int pdofs = opt.getInt("pdofs");
 
        const gsSparseMatrix<T, MatOrder>& matNS = mat.at("matNS");
 
        gsSparseMatrix<T, MatOrder> blockA = matNS.block(0, 0, uSize, uSize);
        gsSparseMatrix<T, MatOrder> Dinv(uSize, uSize); // approximation of A inverse
        diagInvMatrix_into(blockA, Dinv, 1, opt.getSwitch("lumpingA"));
 
        m_mat = - matNS.block(uSize, 0, pdofs, uSize) * Dinv * matNS.block(0, uSize, uSize, pdofs);
 
        m_pSolver->compute(m_mat);
    }
 
public: // *** Static functions ***
 
    static uPtr make(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondSchurSIMPLE<T, MatOrder>(mat, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        GISMO_ASSERT(input.rows() == this->rows(), "Wrong input size.");
        x.resize(this->rows(), 1);
 
        x = m_pSolver->solve(input);
    }
 
public: // *** Getters/setters ***
 
    int rows() const { return m_mat.rows(); }
 
    int cols() const { return m_mat.cols(); }
 
}; // class gsINSPrecondSchurSIMPLE
 
 
// === Schur MSIMPLER ==== // 
 
template <class T, int MatOrder>
class gsINSPrecondSchurMSIMPLER : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    gsSparseMatrix<T, MatOrder> m_mat;
 
protected: // *** Base class members ***
 
    using Base::m_pSolver;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondSchurMSIMPLER> Ptr;
    typedef memory::unique_ptr<gsINSPrecondSchurMSIMPLER> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondSchurMSIMPLER(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt) :
        Base(opt)
    {
        // S = - B velMinv Bt
        int dim = opt.getInt("dim");
        int udofs = opt.getInt("udofs");
        int uSize = dim * udofs;
        int pdofs = opt.getInt("pdofs");
 
        const gsSparseMatrix<T, MatOrder>& matNS = mat.at("matNS");
        const gsSparseMatrix<T, MatOrder>& velM = mat.at("matMu");
 
        gsSparseMatrix<T, MatOrder> velMinv(uSize, uSize); // approximation of velocity mass matrix inverse
        diagInvMatrix_into(velM, velMinv, uSize / velM.rows(), opt.getSwitch("lumpingM"));
 
        m_mat = - matNS.block(uSize, 0, pdofs, uSize) * velMinv * matNS.block(0, uSize, uSize, pdofs);
 
        m_pSolver->compute(m_mat);
    }
 
public: // *** Static functions ***
 
    static uPtr make(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondSchurMSIMPLER<T, MatOrder>(mat, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        GISMO_ASSERT(input.rows() == this->rows(), "Wrong input size.");
        x.resize(this->rows(), 1);
 
        x = m_pSolver->solve(input);
    }
 
public: // *** Getters/setters ***
 
    int rows() const { return m_mat.rows(); }
 
    int cols() const { return m_mat.cols(); }
 
}; // class gsINSPrecondSchurMSIMPLER
 
 
// === Schur Stokes ==== // 
 
template <class T, int MatOrder>
class gsINSPrecondSchurStokes : public gsINSPrecondBlock<T>
{
 
public:
    typedef gsINSPrecondBlock<T> Base;
 
protected: // *** Class members ***
 
    gsSparseMatrix<T, MatOrder> m_mat;
    T m_visc;
 
public: // *** Smart pointers ***
 
    typedef memory::shared_ptr<gsINSPrecondSchurStokes> Ptr;
    typedef memory::unique_ptr<gsINSPrecondSchurStokes> uPtr;
 
public: // *** Constructor/destructor ***
 
    gsINSPrecondSchurStokes(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        // S = diag(Mp)
 
        m_visc = opt.getReal("visc");
 
        const gsSparseMatrix<T, MatOrder>& presM = mat.at("matMp");
        diagInvMatrix_into(presM, m_mat, 1, opt.getSwitch("lumpingM"));
    }
 
public: // *** Static functions ***
 
    static uPtr make(const std::map<std::string, gsSparseMatrix<T, MatOrder> >& mat, const gsOptionList& opt)
    {
        return memory::make_unique(new gsINSPrecondSchurStokes<T, MatOrder>(mat, opt));
    }
 
public: // *** Member functions ***
 
    void apply(const gsMatrix<T> & input, gsMatrix<T> & x) const
    {
        GISMO_ASSERT(input.rows() == this->rows(), "Wrong input size.");
        x.resize(this->rows(), 1);
 
        x = - m_visc * m_mat * input;
    }
 
public: // *** Getters/setters ***
 
    int rows() const { return m_mat.rows(); }
 
    int cols() const { return m_mat.cols(); }
 
}; // class gsINSPrecondSchurStokes
 
 
} // namespace gismo
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsINSPrecondBlocks.hpp)
#endif