gsElTimeIntegrator.h

Go to the documentation of this file.

#pragma once

#include <gsElasticity/gsBaseAssembler.h>
#include <gsElasticity/gsBaseUtils.h>

namespace gismo
{

template <class T>
class gsElasticityAssembler;
template <class T>
class gsMassAssembler;

template <class T>
class gsElTimeIntegrator : public gsBaseAssembler<T>
{
public:
    typedef gsBaseAssembler<T> Base;
    gsElTimeIntegrator(gsElasticityAssembler<T> & stiffAssembler_,
                       gsMassAssembler<T> & massAssembler_);

    static gsOptionList defaultOptions();

    void setDisplacementVector(const gsMatrix<T> & displacementVector)
    {
        GISMO_ENSURE(displacementVector.rows() == massAssembler.numDofs(),
                     "Wrong size of the displacement vector: " + util::to_string(displacementVector.rows()) +
                     ". Must be: " + util::to_string(massAssembler.numDofs()));
        solVector.middleRows(0,massAssembler.numDofs()) = displacementVector;
        initialized = false;
    }

    void setVelocityVector(const gsMatrix<T> & velocityVector)
    {
        GISMO_ENSURE(velocityVector.rows() == massAssembler.numDofs(),
                     "Wrong size of the velocity vector: " + util::to_string(velocityVector.rows()) +
                     ". Must be: " + util::to_string(massAssembler.numDofs()));
        velVector = velocityVector;
        initialized = false;
    }

    void makeTimeStep(T timeStep);

    virtual bool assemble(const gsMatrix<T> & solutionVector,
                          const std::vector<gsMatrix<T> > & fixedDoFs);

    virtual int numDofs() const { return stiffAssembler.numDofs(); }

    const gsMatrix<T> & solutionVector() const { return solVector; }

    //const gsMatrix<T> & displacementVector() const { return solVector.middleRows(0,massAssembler.numDofs()); }

    const gsMatrix<T> & velocityVector() const { return velVector; }

    void saveState();

    void recoverState();

    index_t numberIterations() const { return numIters;}

    void constructSolution(gsMultiPatch<T> & displacement) const;

    void constructSolution(gsMultiPatch<T> & displacement, gsMultiPatch<T> & pressure) const;

    gsBaseAssembler<T> & mAssembler() { return massAssembler; }
    gsBaseAssembler<T> & assembler() { return stiffAssembler; }

protected:
    void initialize();

    gsMatrix<T> implicitLinear();
    gsMatrix<T> implicitNonlinear();

    T alpha1() {return 1./m_options.getReal("Beta")/pow(tStep,2); }
    T alpha2() {return 1./m_options.getReal("Beta")/tStep; }
    T alpha3() {return (1-2*m_options.getReal("Beta"))/2/m_options.getReal("Beta"); }
    T alpha4() {return m_options.getReal("Gamma")/m_options.getReal("Beta")/tStep; }
    T alpha5() {return 1 - m_options.getReal("Gamma")/m_options.getReal("Beta"); }
    T alpha6() {return (1-m_options.getReal("Gamma")/m_options.getReal("Beta")/2)*tStep; }

protected:
    gsElasticityAssembler<T> & stiffAssembler;
    gsMassAssembler<T> & massAssembler;
    bool initialized;
    T tStep;
    gsMatrix<T> solVector;
    gsMatrix<T> velVector;
    gsMatrix<T> accVector;
    using Base::m_system;
    using Base::m_options;
    using Base::m_ddof;
    index_t numIters;
    bool hasSavedState;
    gsMatrix<T> solVecSaved;
    gsMatrix<T> velVecSaved;
    gsMatrix<T> accVecSaved;
    std::vector<gsMatrix<T> > ddofsSaved;
    gsMatrix<T> newSolVector, oldVelVector, dispVectorDiff;
    gsSparseMatrix<T> tempMassBlock;
};

}

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