gsHeatEquation.h

Go to the documentation of this file.

#pragma once
#include <gsAssembler/gsVisitorMass.h>

namespace gismo
{

template <class T>
class gsHeatEquation : public gsAssembler<T>
{
public:
    typedef gsAssembler<T> Base;

public:

    explicit gsHeatEquation(gsAssembler<T> & stationary)
    :  Base(stationary),  // note: unnecessary sliced copy here
       m_stationary(&stationary), m_theta(0.5)
    {
        m_options.addReal("theta",
        "Theta parameter determining the time integration scheme [0..1]", m_theta);
    }

public:

    void setTheta(const T th)
    {
        GISMO_ASSERT(th<=1 && th>=0, "Invalid value");
        m_theta= th;
        m_options.setReal("theta", m_theta);
    }
    
    void assemble()
    {
        // Grab theta once and for all
        m_theta = m_options.getReal("theta");
        
        // Assemble the stationary problem
        m_stationary->assemble();

        //copy the Dirichlet values, to enable calling
        // the construct solution functions
        Base::m_ddof = m_stationary->allFixedDofs();
        
        // Assemble mass matrix
        assembleMass();

        GISMO_ASSERT( m_stationary->matrix().rows() == m_mass.rows(),
                      "Something went terribly wrong.");
    }

    void nextTimeStep(const gsMatrix<T> & curSolution, const T Dt);
    
    void nextTimeStep(const gsSparseMatrix<T> & sysMatrix,
                      const gsSparseMatrix<T> & massMatrix,
                      const gsMatrix<T> & rhs0,
                      const gsMatrix<T> & rhs1,                                  
                      const gsMatrix<T> & curSolution,
                      const T Dt);
    
    void nextTimeStepFixedRhs(const gsSparseMatrix<T> & sysMatrix,
                              const gsSparseMatrix<T> & massMatrix,
                              const gsMatrix<T> & rhs,
                              const gsMatrix<T> & curSolution,
                              const T Dt);

    const gsSparseMatrix<T> & mass() const { return m_mass; }
    const gsSparseMatrix<T> & stationaryMatrix() const { return m_stationary->matrix(); }
    const gsSparseMatrix<T> & stationaryRhs() const { return m_stationary->rhs(); }
    
    void assembleMass();

protected:

    using Base::m_options;

    using Base::m_system;

    gsAssembler<T> * m_stationary;
    
    gsSparseMatrix<T> m_mass;
    
    T m_theta;
    
    using Base::m_pde_ptr;
    using Base::m_bases;

};// end class definition


} // namespace gismo


namespace gismo
{

template<class T>
void gsHeatEquation<T>::nextTimeStep(const gsMatrix<T> & curSolution, const T Dt)
{
    nextTimeStepFixedRhs(m_stationary->matrix(), m_mass,
                         m_stationary->rhs(), curSolution, Dt);
}

/*
template<class T>
void gsHeatEquation<T>::nextTimeStep(const gsMatrix<T> & curSolution, const T Dt)
{
    // Keep previous time
    m_prevRhs.swap(m_stationary->rhs());
    // Get current time
    m_time += Dt;
    m_stationary->rhs(m_time);

    // note: m_stationary->matrix() assumed constant in time
        
    nextTimeStep(m_stationary->matrix(), m_mass,
                 m_prevRhs, m_stationary->rhs(), curSolution, Dt);
    m_prevRhs.swap(m_stationary->rhs());
}
*/

template<class T>
void gsHeatEquation<T>::nextTimeStep(const gsSparseMatrix<T> & sysMatrix,
                                     const gsSparseMatrix<T> & massMatrix,
                                     const gsMatrix<T> & rhs0,
                                     const gsMatrix<T> & rhs1,
                                     const gsMatrix<T> & curSolution,
                                     const T Dt)
{
    GISMO_ASSERT( curSolution.rows() == massMatrix.cols(),
                  "Wrong size in current solution vector.");

    const T c1 = Dt * m_theta;
    m_system.matrix() = massMatrix + c1 * sysMatrix;

    const T c2 = Dt * (1.0 - m_theta);
    m_system.rhs().noalias() = c1 * rhs1 + c2 * rhs0 + (massMatrix - c2 * sysMatrix) * curSolution;
}

template<class T>
void gsHeatEquation<T>::nextTimeStepFixedRhs(const gsSparseMatrix<T> & sysMatrix,
                                             const gsSparseMatrix<T> & massMatrix,
                                             const gsMatrix<T> & rhs,
                                             const gsMatrix<T> & curSolution,
                                             const T Dt)
{
    GISMO_ASSERT( curSolution.rows() == massMatrix.cols(),
                  "Wrong size in current solution vector.");

    const T c1 = Dt * m_theta;
    m_system.matrix() = massMatrix + c1 * sysMatrix;

    const T c2 = Dt * (1.0 - m_theta);
    m_system.rhs().noalias() = Dt * rhs + (massMatrix - c2 * sysMatrix) * curSolution;
}


template<class T>
void gsHeatEquation<T>::assembleMass()
{
    // the mass visitor used to need empty variable
    // Base::m_ddof.resize(1, gsMatrix<T>() );

    const index_t m_dofs = this->numDofs();
    if ( 0  == m_dofs ) // Are there any interior dofs ?
    {
        gsWarn<<"No interior DoFs for mass compuation.\n";
        return;
    }

    // Pre-allocate non-zero elements for each column of the
    // sparse matrix
    m_system.reserve(m_bases[0], m_options, 1);// zero rhs's

    // Assemble mass integrals
    gsVisitorMass<T> mass;
    for (size_t np=0; np < m_pde_ptr->domain().nPatches(); ++np )
    {
        //Assemble mass matrix for this patch
        this->apply(mass, np);
    }

    // Assembly is done, compress the matrix
    this->finalize();

    // Store the mass matrix once and for all
    m_system.matrix().swap(m_mass);
}

} // namespace gismo