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 const 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