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