gsElTimeIntegrator.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsElasticity/gsBaseAssembler.h>
#include <gsElasticity/gsBaseUtils.h>
#include <gsElasticity/gsElasticityAssembler.h>
#include <gsElasticity/gsMassAssembler.h>
 
namespace gismo
{
 
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);
 
    using Base::assemble;
    virtual bool assemble(const gsMatrix<T> & solutionVector,
                          const std::vector<gsMatrix<T> > & fixedDoFs);
 
    virtual int numDofs() const;
 
    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;
 
    using Base::constructSolution;
    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