gsIetiSystem.h

Go to the documentation of this file.

#pragma once

#include <gsSolver/gsMatrixOp.h>

namespace gismo
{

template< typename T >
class gsIetiSystem
{
    typedef gsLinearOperator<T>               Op;              
    typedef memory::shared_ptr<Op>            OpPtr;           
    typedef gsSparseMatrix<T>                 SparseMatrix;    
    typedef gsMatrixOp<SparseMatrix>          SparseMatrixOp;  
    typedef gsSparseMatrix<T,RowMajor>        JumpMatrix;      
    typedef memory::shared_ptr<JumpMatrix>    JumpMatrixPtr;   
    typedef gsMatrix<T>                       Matrix;          
public:

    void reserve(index_t n);

    void addSubdomain(JumpMatrixPtr jumpMatrix, OpPtr localMatrixOp,
        Matrix localRhs, OpPtr localSolverOp = OpPtr());

    JumpMatrixPtr&       jumpMatrix(index_t k)           { return m_jumpMatrices[k];   }
    const JumpMatrixPtr& jumpMatrix(index_t k) const     { return m_jumpMatrices[k];   }

    OpPtr&               localMatrixOp(index_t k)        { return m_localMatrixOps[k]; }
    const OpPtr&         localMatrixOp(index_t k) const  { return m_localMatrixOps[k]; }

    Matrix&              localRhs(index_t k)             { return m_localRhs[k];       }
    const Matrix&        localRhs(index_t k) const       { return m_localRhs[k];       }

    OpPtr&               localSolverOp(index_t k)        { return m_localSolverOps[k]; }
    const OpPtr&         localSolverOp(index_t k) const  { return m_localSolverOps[k]; }

    index_t nLagrangeMultipliers() const
    {
        GISMO_ASSERT( m_jumpMatrices.size()>0, "gsIetiSystem: Number of Lagrange multipliers "
            "can only be determined if there are jump matrices.");
        return m_jumpMatrices[0]->rows();
    }

    OpPtr schurComplement() const;

    Matrix rhsForSchurComplement() const;

    std::vector<Matrix> constructSolutionFromLagrangeMultipliers(const Matrix& multipliers) const;

    OpPtr saddlePointProblem() const;

    Matrix rhsForSaddlePoint() const;

    std::vector<Matrix> constructSolutionFromSaddlePoint(const Matrix& x) const;


private:
    void setupSparseLUSolvers() const;                

    std::vector<JumpMatrixPtr>  m_jumpMatrices;       
    std::vector<OpPtr>          m_localMatrixOps;     
    std::vector<Matrix>         m_localRhs;           
    mutable std::vector<OpPtr>  m_localSolverOps;     
};

} // namespace gismo

#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsIetiSystem.hpp)
#endif