gsAdditiveOp.h

Go to the documentation of this file.

#pragma once

#include <gsSolver/gsLinearOperator.h>

namespace gismo
{


template<class T>
class gsAdditiveOp GISMO_FINAL : public gsLinearOperator<T>
{
    typedef typename gsLinearOperator<T>::Ptr   OpPtr;
    typedef std::vector<OpPtr>                  OpContainer;
    typedef gsSparseMatrix<T,RowMajor>          Transfer;
    typedef memory::shared_ptr<Transfer>        TransferPtr;
    typedef std::vector<Transfer>               TransferContainer;
    typedef std::vector<TransferPtr>            TransferPtrContainer;

public:

    typedef memory::shared_ptr<gsAdditiveOp> Ptr;

    typedef memory::unique_ptr<gsAdditiveOp> uPtr;

    gsAdditiveOp() : m_transfers(), m_ops() {}

    gsAdditiveOp(TransferContainer transfers, OpContainer ops)
    : m_transfers(), m_ops(give(ops))
    {
        const size_t sz = transfers.size();
        m_transfers.reserve(sz);
        for (size_t i=0; i<sz; ++i)
            m_transfers.push_back( transfers[i].moveToPtr() );
#ifndef NDEBUG
        GISMO_ASSERT( m_transfers.size() == m_ops.size(), "Sizes do not agree" );
        for (size_t i=0; i<sz; ++i)
        {
            GISMO_ASSERT ( m_transfers[i]->rows()==m_transfers[0]->rows()
                       && m_transfers[i]->cols() == m_ops[i]->rows()
                       && m_ops[i]->cols() == m_ops[i]->rows(),
                       "Dimensions of the operators do not fit." );
        }
#endif
    }
    
    gsAdditiveOp(TransferPtrContainer transfers, OpContainer ops)
    : m_transfers(give(transfers)), m_ops(give(ops))
    {
#ifndef NDEBUG
        GISMO_ASSERT( m_transfers.size() == m_ops.size(), "Sizes do not agree" );
        const size_t sz = m_transfers.size();
        for (size_t i=0; i<sz; ++i)
        {
            GISMO_ASSERT ( m_transfers[i]->rows()==m_transfers[0]->rows()
                       && m_transfers[i]->cols() == m_ops[i]->rows()
                       && m_ops[i]->cols() == m_ops[i]->rows(),
                       "Dimensions of the operators do not fit." );
        }
#endif
    }

    static uPtr make()
    { return uPtr( new gsAdditiveOp() ); }

    static uPtr make(TransferContainer transfers, OpContainer ops)
    { return uPtr( new gsAdditiveOp( give(transfers), give(ops) ) ); }

    static uPtr make(TransferPtrContainer transfers, OpContainer ops)
    { return uPtr( new gsAdditiveOp( give(transfers), give(ops) ) ); }


    void addOperator(Transfer transfer, OpPtr op)
    {
        m_transfers.push_back(transfer.moveToPtr());
        m_ops.push_back(give(op));
        GISMO_ASSERT ( m_transfers.back()->rows()==m_transfers[0]->rows()
                       && m_transfers.back()->cols() == m_ops.back()->rows()
                       && m_ops.back()->cols() == m_ops.back()->rows(),
                       "Dimensions of the operators do not fit." );
    }
    
    void addOperator(TransferPtr transfer, OpPtr op)
    {
        m_transfers.push_back(give(transfer));
        m_ops.push_back(give(op));
        GISMO_ASSERT ( m_transfers.back()->rows()==m_transfers[0]->rows()
                       && m_transfers.back()->cols() == m_ops.back()->rows()
                       && m_ops.back()->cols() == m_ops.back()->rows(),
                       "Dimensions of the operators do not fit." );
    }

    void apply(const gsMatrix<T>& input, gsMatrix<T>& x) const;

    index_t rows() const
    {
        GISMO_ASSERT( !m_transfers.empty(), "gsAdditiveOp::rows does not work for 0 operators." );
        return m_transfers[0]->rows();
    }

    index_t cols() const
    {
        GISMO_ASSERT( !m_transfers.empty(), "gsAdditiveOp::cols does not work for 0 operators." );
        return m_transfers[0]->rows();
    }

protected:
    TransferPtrContainer m_transfers;   
    OpContainer m_ops;                  

};

} // namespace gismo

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