gsALMBase.h

Go to the documentation of this file.

 
#pragma once
#include <gsCore/gsLinearAlgebra.h>
 
#ifdef gsSpectra_ENABLED
#include <gsSpectra/gsSpectra.h>
#endif
#include <gsIO/gsOptionList.h>
#include <gsStructuralAnalysis/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h>
 
namespace gismo
{
 
template <class T>
class gsALMBase
{
protected:
 
    typedef typename gsStructuralAnalysisOps<T>::ALResidual_t    ALResidual_t;
    typedef typename gsStructuralAnalysisOps<T>::Jacobian_t      Jacobian_t;
    typedef typename gsStructuralAnalysisOps<T>::dJacobian_t     dJacobian_t;
 
public:
 
    virtual ~gsALMBase() {};
 
    gsALMBase(  const Jacobian_t   & Jacobian,
                const ALResidual_t & ALResidual,
                const gsVector<T>  & Force )
    : m_residualFun(ALResidual),
      m_forcing(Force)
    {
        m_jacobian  = Jacobian;
        m_djacobian = [this](gsVector<T> const & x, gsVector<T> const & /*dx*/, gsSparseMatrix<T> & m) -> bool
        {
            return m_jacobian(x,m);
        };
 
        // initialize variables
        m_numIterations = 0;
        this->defaultOptions();
        this->setLength(1e-2);
        m_converged = false;
 
        // initialize errors
        m_basisResidualF = 0.0;
        m_basisResidualU = 0.0;
 
        m_status = gsStatus::NotStarted;
    }
 
    gsALMBase(  const dJacobian_t &  dJacobian,
                const ALResidual_t & Residual,
                const gsVector<T>  & Force )
    : m_residualFun(Residual),
      m_forcing(Force)
    {
        m_djacobian  = dJacobian;
 
        // initialize variables
        m_numIterations = 0;
        this->defaultOptions();
        this->setLength(1e-2);
        m_converged = false;
 
        // initialize errors
        m_basisResidualF = 0.0;
        m_basisResidualU = 0.0;
    }
 
// General functions
public:
 
    virtual gsStatus status() { return m_status; }
 
    // Returns the number of DoFs
    virtual index_t numDofs() {return m_forcing.size();}
 
    // Returns the current length
    virtual T getLength() {return m_arcLength; }
 
    virtual gsStatus step();
 
    virtual void initialize(bool stability = true)
    {
        m_initialized = true;
        this -> initMethods();
        this -> init(stability);
    }
 
    virtual void setLength(T length)
    {
      m_options.setReal("Length",length);
      m_arcLength = m_arcLength_prev = m_arcLength_ori = m_options.getReal("Length");
    }
 
    virtual void setLength(T length, bool adaptive)
    {
      this->setLength(length);
      m_adaptiveLength = adaptive;
      m_desiredIterations = 10; // number of desired iterations defaults to 10
    }
 
    virtual void setLength(T length, bool adaptive, index_t iterations)
    {
      this->setLength(length);
      m_adaptiveLength = adaptive;
      m_desiredIterations = iterations;
    }
    virtual void setLength(T length, index_t iterations)
    {
      this->setLength(length);
      m_adaptiveLength = true;
      m_desiredIterations = iterations;
    }
 
    // Output
    virtual bool converged() const {return m_converged;}
 
    virtual index_t numIterations() const { return m_numIterations;}
 
    virtual T tolerance() const {return m_tolerance;}
 
    virtual T residue()   const {return m_residue;}
 
    virtual T indicator() const {return m_indicator;}
    virtual void setIndicator(T indicator) {m_indicator = indicator; m_stability = this->stability();}
 
    virtual const gsVector<T> & solutionU() const {return  m_U;}
    virtual const gsVector<T> & solutionDU() const {return  m_DeltaU;}
    virtual T  solutionL() const {return  m_L;}
    virtual T  solutionDL() const {return  m_DeltaL;}
    virtual const gsVector<T> & solutionV() const {return  m_V;}
 
    virtual bool isStable() const {return m_stability;}
 
    // /// Returns the value of the deterimant of the jacobian
    // virtual T determinant() const
    // {
    //     return m_jacobian(m_U).toDense().determinant();
    // }
 
    virtual void resetStep() {m_DeltaUold.setZero(); m_DeltaLold = 0;}
 
    // Set initial guess for solution
    virtual void setInitialGuess(const gsVector<T> & Uguess, const T & Lguess) {m_Uguess = Uguess; m_Lguess = Lguess;}
    virtual void setPrevious(const gsVector<T> & Uprev, const T & Lprev)
    {
        m_Uprev = Uprev;
        m_Lprev = Lprev;
        m_DeltaUold = m_U - m_Uprev;
        m_DeltaLold = m_L - m_Lprev;
    }
 
    virtual void setSolution(const gsVector<T> & U, const T & L) {m_L = L; m_U = U; }// m_DeltaUold.setZero(); m_DeltaLold = 0;}
 
    virtual void setSolutionStep(const gsVector<T> & DU, const T & DL) {m_DeltaUold = DU; m_DeltaLold = DL;}// m_DeltaUold.setZero(); m_DeltaLold = 0;}
 
    virtual gsOptionList & options() {return m_options;};
 
    virtual void setOptions(gsOptionList options) {m_options.update(options,gsOptionList::addIfUnknown); this->getOptions(); };
 
    virtual void options_into(gsOptionList options) {options = m_options;};
 
    virtual void applyOptions() {this->getOptions(); }
 
    virtual T distance(const gsVector<T>& /*DeltaU*/, const T /*DeltaL*/) const
    {
        GISMO_NO_IMPLEMENTATION;
    }
 
// ------------------------------------------------------------------------------------------------------------
// ---------------------------------------Singular point methods-----------------------------------------------
// ------------------------------------------------------------------------------------------------------------
public:
 
 
    virtual gsStatus computeSingularPoint(const gsVector<T> & U, const T & L, bool switchBranch=false, bool jacobian=false, bool testPoint=true);
 
    virtual bool isBifurcation(bool jacobian = false);
 
    virtual bool stabilityChange() const;
 
    virtual gsStatus computeStability(bool jacobian=true, T shift = -1e2);
 
    virtual index_t stability() const;
 
    virtual void switchBranch();
 
    virtual T reduceLength(T fac = 0.5);
 
    virtual T resetLength();
 
protected:
 
    virtual void _computeStability(const gsVector<T> & x, bool jacobian=true, T shift = -1e2);
 
    virtual void _computeSingularPoint(const gsVector<T> & U, const T & L, bool switchBranch=false, bool jacobian=false, bool testPoint=true);
 
    virtual bool _testSingularPoint(bool jacobian=false);
 
    virtual void _extendedSystemIteration();
 
    virtual index_t _bisectionObjectiveFunction(const gsVector<T> & x, bool jacobian=true);
    virtual T       _bisectionTerminationFunction(const gsVector<T> & x, bool jacobian=true);
 
    virtual bool _extendedSystemSolve(const gsVector<T> & U, const T L, const T tol);
 
    virtual bool _bisectionSolve(const gsVector<T> & U, const T L, const T tol);
 
    virtual void _initOutputExtended();
 
    virtual void _stepOutputExtended();
 
// ------------------------------------------------------------------------------------------------------------
// ---------------------------------------Computations---------------------------------------------------------
// ------------------------------------------------------------------------------------------------------------
protected:
    virtual void _step();
 
    virtual void defaultOptions();
 
    virtual void getOptions();
 
    virtual void init(bool stability);
 
    virtual void factorizeMatrix(const gsSparseMatrix<T> & M);
 
    virtual gsVector<T> solveSystem(const gsVector<T> & F);
 
    virtual gsVector<T> computeResidual(const gsVector<T> & U, const T & L);
    virtual void computeResidual();
 
    virtual void computeResidualNorms();
 
    virtual gsSparseMatrix<T> _computeJacobian(const gsVector<T> & U, const gsVector<T> & dU);
    virtual gsSparseMatrix<T> computeJacobian(const gsVector<T> & U, const gsVector<T> & dU);
    virtual gsSparseMatrix<T> computeJacobian(const gsVector<T> & U);
    virtual gsSparseMatrix<T> computeJacobian();
 
    virtual void computeLength();
 
    virtual void computeUt();
    virtual void computeUbar();
 
// Purely virtual functions
protected:
    virtual void initMethods() = 0;
    virtual void initiateStep() = 0;
    virtual void iterationFinish() = 0;
 
    virtual void quasiNewtonPredictor() = 0;
    virtual void quasiNewtonIteration() = 0;
 
    virtual void predictor() = 0;
    virtual void predictorGuess() = 0;
    virtual void iteration() = 0;
 
    virtual void initOutput() = 0;
    virtual void stepOutput() = 0;
 
protected:
 
    // Number of degrees of freedom
    index_t m_numDof;
 
    Jacobian_t      m_jacobian;
    dJacobian_t     m_djacobian;
    const ALResidual_t    m_residualFun;
    const gsVector<T>     m_forcing;
 
    mutable typename gsSparseSolver<T>::uPtr m_solver; // Cholesky by default
 
public:
 
 
    struct bifmethod
    {
        enum type
        {
            Nothing = -1,
            Determinant = 0,
            Eigenvalue  = 1,
        };
    };
 
    struct solver
    {
        enum type
        {
            LDLT = 0,
            CG  = 1, // The CG solver is robust for membrane models, where zero-blocks in the matrix might occur.
        };
    };
 
    struct SPfail
    {
        enum type
        {
            Without  = 0,
            With     = 1,
        };
    };
 
  protected:
    mutable gsOptionList m_options;
 
 
    index_t m_numIterations;
 
    index_t m_maxIterations;
 
    index_t m_desiredIterations;
 
    T m_arcLength;
    T m_arcLength_prev;
    T m_arcLength_ori;
    bool m_adaptiveLength;
 
    T m_tolerance;
 
    T m_toleranceF;
 
    T m_toleranceU;
 
    bool m_verbose;
    bool m_initialized;
 
    bool m_quasiNewton;
    index_t m_quasiNewtonInterval;
 
    std::string m_note;
 
    gsStatus m_status;
 
protected:
 
    bool m_converged;
 
    T m_residue;
 
    T m_residueF;
    T m_basisResidualF;
 
    T m_residueU;
    T m_basisResidualU;
 
    T m_residueL;
 
    T m_residueKTPhi;
    T m_basisResidualKTPhi;
 
    T m_indicator;
    index_t m_negatives;
 
    T m_relax;
 
protected:
 
    // Previous update
    gsVector<T> m_DeltaUold;
    real_t m_DeltaLold;
    gsVector<T> m_U, m_Uprev;
    gsVector<T> m_DeltaU;
    gsVector<T> m_deltaUbar;
    gsVector<T> m_deltaUt;
    gsVector<T> m_deltaU;
 
    T m_L, m_Lprev;
    T m_DeltaL;
    T m_deltaL;
    gsVector<T> m_deltaLs;
 
    gsVector<T> m_Uguess;
    T m_Lguess;
 
    gsSparseMatrix<T> m_jacMat;
    T m_detKT;
 
    gsVector<T> m_resVec;
 
    // eigenvector
    gsVector<T> m_V;
    // step eigenvector
    gsVector<T> m_deltaV;
    gsVector<T> m_deltaVbar;
    gsVector<T> m_deltaVt;
    gsVector<T> m_DeltaV;
 
    // Stability indicator
    gsVector<T> m_stabilityVec;
 
    // Integer if point is unstable (-1) or not (+1)
    index_t m_stabilityPrev; //previous step
    index_t m_stability; //current step
    // Method to check if a point is a bifurcation point
    index_t m_bifurcationMethod;
 
    // What to do after computeSingularPoint fails?
    index_t m_SPfail;
 
    // Number of iterations and the tolerance for the singular point test
    index_t m_SPTestIt;
    T m_SPTestTol;
 
    // Singular point computation tolerances
    T m_SPCompTolE; // extended iterations
    T m_SPCompTolB; // bisection method
 
    // Branch switch parameter
    T m_tau;
};
 
 
} // namespace gismo
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsALMBase.hpp)
#endif