gsINSAssembler.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsCore/gsField.h>
#include <gsIncompressibleFlow/src/gsFlowAssemblerBase.h>
#include <gsIncompressibleFlow/src/gsINSVisitors.h>
#include <gsIncompressibleFlow/src/gsFlowUtils.h>
 
namespace gismo
{
 
template<class T, int MatOrder>
class gsINSAssembler: public gsFlowAssemblerBase<T, MatOrder>
{
 
public:
    typedef gsFlowAssemblerBase<T, MatOrder> Base;
 
 
protected: // *** Class members ***
 
    T m_viscosity;
    gsSparseMatrix<T, MatOrder> m_baseMatrix, m_matrix;
    gsMatrix<T> m_baseRhs, m_rhs;
 
    index_t m_udofs;
    index_t m_pdofs;
    index_t m_pshift;
    index_t m_nnzPerRowU, m_nnzPerRowP;
    gsINSVisitorUUlin<T, MatOrder> m_visitorUUlin;
    gsINSVisitorUUnonlin<T, MatOrder> m_visitorUUnonlin;
    gsINSVisitorPU_withUPrhs<T, MatOrder> m_visitorUP;
    gsINSVisitorRhsU<T, MatOrder> m_visitorF;
    gsINSVisitorRhsP<T, MatOrder> m_visitorG;
    gsSparseMatrix<T, MatOrder> m_blockUUlin_comp, m_blockUUnonlin_comp, m_blockUP;
    gsSparseMatrix<T, MatOrder> m_blockUUlin_whole, m_blockUUnonlin_whole;
    gsMatrix<T> m_rhsUlin, m_rhsUnonlin, m_rhsBtB, m_rhsFG;
    gsField<T>  m_currentVelField, m_currentPresField;
 
    bool m_isMassMatReady;
    std::vector< gsSparseMatrix<T, MatOrder> > m_massMatBlocks;
 
    // PCD members
    // std::vector<index_t> m_presInIDs, m_presOutIDs, m_presWallIDs
    // std::vector< gsSparseMatrix<T, MatOrder> > m_pcdBlocks(3); // [laplaceP, robinBCblock, convectionP]
 
 
protected: // *** Base class members ***
 
    using Base::m_params;
    using Base::m_dofs;
    using Base::m_tarDim;
    using Base::m_dofMappers;
    using Base::m_ddof;
    using Base::m_solution;
    using Base::m_isBaseReady;
    using Base::m_isSystemReady;
 
public: // *** Base class member functions ***
 
    using Base::getBases;
    using Base::getPatches;
    using Base::getAssemblerOptions;
    using Base::getBCs;
 
public: // *** Constructor/destructor ***
 
    gsINSAssembler(const gsFlowSolverParams<T>& params):
    Base(params)
    { }
 
    virtual ~gsINSAssembler()
    { }
 
protected: // *** Member functions ***
 
    void initMembers();
 
    virtual void updateSizes();
 
    virtual void updateDofMappers();
 
    virtual void updateCurrentSolField(const gsMatrix<T> & solVector, bool updateSol);
 
    virtual void assembleLinearPart();
 
    virtual void assembleNonlinearPart();
 
    void fillGlobalMat_UU(gsSparseMatrix<T, MatOrder>& globalMat, const gsSparseMatrix<T, MatOrder>& sourceMat);
 
    void fillGlobalMat_UP(gsSparseMatrix<T, MatOrder>& globalMat, const gsSparseMatrix<T, MatOrder>& sourceMat);
 
    void fillGlobalMat_PU(gsSparseMatrix<T, MatOrder>& globalMat, const gsSparseMatrix<T, MatOrder>& sourceMat);
 
    void fillGlobalMat_PP(gsSparseMatrix<T, MatOrder>& globalMat, const gsSparseMatrix<T, MatOrder>& sourceMat);
 
    virtual void fillBaseSystem();
 
    virtual void fillSystem();
 
 
    // --- PCD member functions ---
 
    // void findPressureBoundaryIDs();
 
    // void findPressureBoundaryPartIDs(std::vector<std::pair<int, boxSide> > bndPart, std::vector<index_t>& idVector);
 
 
public: // *** Member functions ***
 
    virtual void initialize();
 
    virtual void update(const gsMatrix<T> & solVector, bool updateSol = true);
 
    virtual void markDofsAsEliminatedZeros(const std::vector< gsMatrix< index_t > > & boundaryDofs, const index_t unk);
 
    virtual void fillStokesSystem(gsSparseMatrix<T, MatOrder>& stokesMat, gsMatrix<T>& stokesRhs);
 
    virtual gsField<T> constructSolution(const gsMatrix<T>& solVector, index_t unk) const;
 
    T computeFlowRate(index_t patch, boxSide side, gsMatrix<T> solution) const;
 
 
public: // *** Getters/setters ***
 
    T getViscosity() const { return m_viscosity; }
 
    virtual const gsSparseMatrix<T, MatOrder>& matrix() const
    {
        GISMO_ASSERT(m_isSystemReady, "Matrix not ready, is fillGlobalSyst in the solver params true? If so, update() must be called first.");
        return m_matrix;
    }
 
    virtual const gsMatrix<T>& rhs() const
    {
        GISMO_ASSERT(m_isSystemReady, "Rhs not ready, is fillGlobalSyst in the solver params true? If so, update() must be called first.");
        return m_rhs;
    }
 
    const gsMatrix<T>& getSolution() const { return m_solution; }
 
    index_t numDofsUnk(index_t i);
 
    index_t getUdofs() const { return m_udofs; }
 
    index_t getPdofs() const { return m_pdofs; }
 
    index_t getPshift() const { return m_pshift; }
 
    virtual gsSparseMatrix<T, MatOrder> getBlockUU(bool linPartOnly = false);
 
    virtual gsSparseMatrix<T, MatOrder> getBlockUUcompDiag(index_t i = 0)
    { 
        GISMO_ASSERT(i >= 0 && i < m_tarDim, "Component index out of range.");
        return getBlockUU().block(i * m_udofs, i * m_udofs, m_udofs, m_udofs); 
    }
 
 
    const gsSparseMatrix<T, MatOrder>& getBlockUP() const
    { return m_blockUP; }
 
    gsSparseMatrix<T, MatOrder> getBlockPU() const
    { return (-1.0)*gsSparseMatrix<T, MatOrder>(m_blockUP.transpose()); }
 
    virtual gsSparseMatrix<T, MatOrder> getBlockUPcomp(index_t i) const
    { 
        GISMO_ASSERT(i >= 0 && i < m_tarDim, "Component index out of range.");
        return getBlockUP().middleRows(i * m_udofs, m_udofs);
    }
 
    virtual gsSparseMatrix<T, MatOrder> getBlockPUcomp(index_t i) const
    { 
        GISMO_ASSERT(i >= 0 && i < m_tarDim, "Component index out of range.");
        return (-1.0)*gsSparseMatrix<T, MatOrder>(getBlockUPcomp(i).transpose());
    }
 
    virtual gsSparseMatrix<T, MatOrder>& getMassMatrix(index_t unkID)
    { 
        GISMO_ASSERT(unkID == 0 || unkID == 1, "unkID must be 0 (velocity) or 1 (pressure).");
        GISMO_ASSERT(m_isMassMatReady, "Mass matrices not assembled in gsINSAssembler.");
        return m_massMatBlocks[unkID];
    }
 
    virtual const gsSparseMatrix<T, MatOrder>& getMassMatrix(index_t unkID) const
    { 
        GISMO_ASSERT(unkID == 0 || unkID == 1, "unkID must be 0 (velocity) or 1 (pressure).");
        GISMO_ASSERT(m_isMassMatReady, "Mass matrices not assembled in gsINSAssembler.");
        return m_massMatBlocks[unkID];
    }
 
    virtual gsMatrix<T> getRhsU() const;
 
    virtual gsMatrix<T> getRhsUcomp(index_t i) const
    { 
        GISMO_ASSERT(i >= 0 && i < m_tarDim, "Component index out of range.");
        return getRhsU().middleRows(i * m_udofs, m_udofs);
    }
 
    gsMatrix<T> getRhsP() const;
 
}; // gsINSAssembler
 
 
// ===================================================================================================================
 
 
template<class T, int MatOrder>
class gsINSAssemblerSteady: public gsINSAssembler<T, MatOrder>
{
 
public:
    typedef gsINSAssembler<T, MatOrder> Base;
 
 
public: // *** Constructor/destructor ***
 
    gsINSAssemblerSteady(const gsFlowSolverParams<T>& params): 
    Base(params)
    {
        Base::initMembers();
    }
 
    virtual ~gsINSAssemblerSteady()
    { }
 
}; // gsINSAssemblerSteady
 
 
// ===================================================================================================================
 
 
template<class T, int MatOrder>
class gsINSAssemblerUnsteady: public gsINSAssembler<T, MatOrder>
{
 
public:
    typedef gsINSAssembler<T, MatOrder> Base;
 
 
protected: // *** Class members ***
 
    gsINSVisitorUUtimeDiscr<T, MatOrder> m_visitorTimeDiscr;
    gsSparseMatrix<T, MatOrder> m_blockTimeDiscr;
    gsMatrix<T> m_rhsTimeDiscr;
    gsField<T> m_oldTimeVelField;
 
 
protected: // *** Base class members ***
 
    using Base::m_params;
    using Base::m_pshift;
    using Base::m_nnzPerRowU;
    using Base::m_dofMappers;
    using Base::m_solution;
    using Base::m_baseMatrix;
    using Base::m_matrix;
    using Base::m_rhs;
    using Base::m_currentVelField;
 
 
public: // *** Constructor/destructor ***
 
    gsINSAssemblerUnsteady(const gsFlowSolverParams<T>& params): 
    Base(params)
    {
        initMembers();
    }
 
    virtual ~gsINSAssemblerUnsteady()
    { }
 
 
protected: // *** Member functions ***
 
    void initMembers();
 
    virtual void updateSizes();
 
    virtual void updateDofMappers();
 
    virtual void updateCurrentSolField(const gsMatrix<T> & solVector, bool updateSol);
 
    virtual void assembleLinearPart();
 
    virtual void fillSystem();
 
 
public: // *** Member functions ***
 
    virtual void update(const gsMatrix<T> & solVector, bool updateSol = true);
 
 
public: // *** Getters/setters ***
 
    virtual gsSparseMatrix<T, MatOrder> getBlockUU()
    { return Base::getBlockUU() + m_blockTimeDiscr; }
 
 
    virtual gsMatrix<T> getRhsU() const
    { return Base::getRhsU() + m_rhsTimeDiscr; }
 
 
}; //gsINSAssemblerUnsteady
 
} // namespace gismo
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsINSAssembler.hpp)
#endif