gsPrimalSystem.h

Go to the documentation of this file.

#pragma once

#include <gsSolver/gsMatrixOp.h>
#include <gsMatrix/gsVector.h>

namespace gismo
{


template< typename T >
class gsPrimalSystem
{
private:
    typedef gsLinearOperator<T>               Op;              
    typedef memory::shared_ptr<Op>            OpPtr;           
    typedef gsSparseMatrix<T>                 SparseMatrix;    
    typedef gsSparseMatrix<T,RowMajor>        JumpMatrix;      
    typedef gsSparseVector<T>                 SparseVector;    
    typedef gsMatrix<T>                       Matrix;          
public:

    gsPrimalSystem(index_t nPrimalDofs);

    static void incorporateConstraints(
        const std::vector<SparseVector>& primalConstraints,
        bool eliminatePointwiseConstraints,
        const SparseMatrix& localMatrix,
        SparseMatrix& modifiedLocalMatrix,
        SparseMatrix& localEmbedding,
        SparseMatrix& embeddingForBasis,
        Matrix& rhsForBasis
    );

    static gsSparseMatrix<T> primalBasis(
        OpPtr localSaddlePointSolver,
        const SparseMatrix& embeddingForBasis,
        const Matrix& rhsForBasis,
        const std::vector<index_t>& primalDofIndices,
        index_t nPrimalDofs
    );

    void addContribution(
        const JumpMatrix& jumpMatrix,
        const SparseMatrix& localMatrix,
        const Matrix& localRhs,
        SparseMatrix primalBasis,
        OpPtr embedding = OpPtr()
    );

    void handleConstraints(
        const std::vector<SparseVector>& primalConstraints,
        const std::vector<index_t>& primalDofIndices,
        JumpMatrix& jumpMatrix,
        SparseMatrix& localMatrix,
        Matrix& localRhs
    );

    std::vector<Matrix> distributePrimalSolution( std::vector<Matrix> sol );

    JumpMatrix&                           jumpMatrix()        { return m_jumpMatrix;                    }
    const JumpMatrix&                     jumpMatrix() const  { return m_jumpMatrix;                    }

    SparseMatrix&                         localMatrix()       { return m_localMatrix;                   }
    const SparseMatrix&                   localMatrix() const { return m_localMatrix;                   }

    Matrix&                               localRhs()          { return m_localRhs;                      }
    const Matrix&                         localRhs() const    { return m_localRhs;                      }

    index_t nPrimalDofs() const                               { return m_localMatrix.rows();            }

    index_t eliminatePointwiseConstraints() const             { return m_eliminatePointwiseConstraints; }

    void setEliminatePointwiseConstraints(bool v)             { m_eliminatePointwiseConstraints = v;    }

private:
    JumpMatrix                  m_jumpMatrix;   
    SparseMatrix                m_localMatrix;  
    Matrix                      m_localRhs;     
    std::vector<SparseMatrix>   m_primalBases;  
    std::vector<OpPtr>          m_embeddings;   
    bool                        m_eliminatePointwiseConstraints; 
};

} // namespace gismo

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