gsScaledDirichletPrec.h

Go to the documentation of this file.

#pragma once

#include <gsSolver/gsMatrixOp.h>
#include <gsUtils/gsSortedVector.h>

namespace gismo
{

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

    void reserve( index_t n )
    {
        m_jumpMatrices.reserve(n);
        m_localSchurOps.reserve(n);
        m_localScaling.reserve(n);
    }

    void addSubdomain( JumpMatrixPtr jumpMatrix, OpPtr localSchurOp )
    {
        m_jumpMatrices.push_back(give(jumpMatrix));
        m_localSchurOps.push_back(give(localSchurOp));
        m_localScaling.push_back(Matrix());
    }

    // Adds a new subdomain
    // Documented in comment above
    void addSubdomain( std::pair<JumpMatrix,OpPtr> data )
    { addSubdomain(data.first.moveToPtr(), give(data.second)); }

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

    OpPtr&               localSchurOps(index_t k)        { return m_localSchurOps[k]; }
    const OpPtr&         localSchurOps(index_t k) const  { return m_localSchurOps[k]; }

    Matrix&              localScaling(index_t k)         { return m_localScaling[k];  }
    const Matrix&        localScaling(index_t k) const   { return m_localScaling[k];  }

    static gsSortedVector<index_t> skeletonDofs( const JumpMatrix& jumpMatrix );

    static JumpMatrix restrictJumpMatrix( const JumpMatrix& jumpMatrix, const std::vector<index_t>& dofs );

    struct Blocks { SparseMatrix A00, A01, A10, A11; };

    static Blocks matrixBlocks( const SparseMatrix& localMatrix, const std::vector<index_t>& dofs );

    static OpPtr schurComplement( Blocks matrixBlocks, OpPtr solver );

    static OpPtr schurComplement( const SparseMatrix& localMatrix, const std::vector<index_t>& dofs )
    {
        Blocks blocks = matrixBlocks(localMatrix, dofs);
        return schurComplement( blocks, makeSparseCholeskySolver(blocks.A11) );
    }

    static std::pair<JumpMatrix,OpPtr> restrictToSkeleton(
        const JumpMatrix& jumpMatrix,
        const SparseMatrix& localMatrix,
        const std::vector<index_t>& dofs
    )
    { return std::pair<JumpMatrix,OpPtr>(restrictJumpMatrix(jumpMatrix,dofs),schurComplement(localMatrix,dofs)); }

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

    void setupMultiplicityScaling();

    OpPtr preconditioner() const;

public:
    std::vector<JumpMatrixPtr>  m_jumpMatrices;     
    std::vector<OpPtr>          m_localSchurOps;    
    std::vector<Matrix>         m_localScaling;     
};

} // namespace gismo

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