gsStaticNewton.h

Go to the documentation of this file.

#include <typeinfo>

#include <gsStructuralAnalysis/src/gsStaticSolvers/gsStaticBase.h>
#pragma once


namespace gismo
{

template <class T>
class gsStaticNewton : public gsStaticBase<T>
{
protected:

    typedef gsStaticBase<T> Base;

    typedef typename Base::Residual_t    Residual_t;
    typedef typename Base::ALResidual_t  ALResidual_t;
    typedef typename Base::Jacobian_t    Jacobian_t;
    typedef typename Base::dJacobian_t   dJacobian_t;

public:

    gsStaticNewton(     const gsSparseMatrix<T> &linear,
                        const gsVector<T> &force        )
    :
        m_linear(linear),
        m_force(force),
        m_nonlinear(nullptr),
        m_residualFun(nullptr)
    {
        this->_init();
    }

    gsStaticNewton( const gsSparseMatrix<T> &linear,
                    const gsVector<T> &force,
                    const Jacobian_t &nonlinear,
                    const Residual_t &residual      )
    :
        m_linear(linear),
        m_force(force),
        m_nonlinear(nonlinear),
        m_dnonlinear(nullptr),
        m_residualFun(residual),
        m_ALresidualFun(nullptr)
    {
        m_dnonlinear = [this](gsVector<T> const & x, gsVector<T> const & dx, gsSparseMatrix<T> & m) -> bool
        {
            return m_nonlinear(x,m);
        };

        this->_init();
    }

    gsStaticNewton( const gsSparseMatrix<T> &linear,
                    const gsVector<T>  &force,
                    const Jacobian_t   &nonlinear,
                    const ALResidual_t &ALResidual  )
    :
        m_linear(linear),
        m_force(force),
        m_nonlinear(nonlinear),
        m_dnonlinear(nullptr),
        m_residualFun(nullptr),
        m_ALresidualFun(ALResidual)
    {
        m_L = 1.0;
        m_residualFun = [this](gsVector<T> const & x, gsVector<T> & result) -> bool
        {
            return m_ALresidualFun(x,m_L,result);
        };

        m_dnonlinear = [this](gsVector<T> const & x, gsVector<T> const & dx, gsSparseMatrix<T> & m) -> bool
        {
            return m_nonlinear(x,m);
        };

        this->_init();
    }

    gsStaticNewton(   const gsSparseMatrix<T> &linear,
                      const gsVector<T> &force,
                      const dJacobian_t &dnonlinear,
                      const Residual_t &residual
                      )
    :
        m_linear(linear),
        m_force(force),
        m_nonlinear(nullptr),
        m_dnonlinear(dnonlinear),
        m_residualFun(residual),
        m_ALresidualFun(nullptr)
    {
        this->_init();
    }

public:

    gsStatus solve() override
    {
        if (m_NL)
            return solveNonlinear();
        else
            return solveLinear();
    }
    gsStatus solveLinear();
    gsStatus solveLinear(gsVector<T> & solution)
    {
        gsStatus status = this->solveLinear();
        solution = m_U;
        return status;
    }
    gsStatus solveNonlinear(gsVector<T> & solution)
    {
        gsStatus status = this->solveNonlinear();
        solution = m_U;
        return status;
    }
    gsStatus solveNonlinear();

    void initOutput() override;

    void stepOutput(index_t k) override;

    void defaultOptions() override;

    void reset() override;

    void getOptions() override;

    using Base::indicator;
    using Base::stabilityVec;

    T indicator()
    {
        gsSparseMatrix<T> m;
        GISMO_ENSURE(m_nonlinear(m_U, m),"Assembly failed");
        return indicator(m);
    }

    gsVector<T> stabilityVec()
    {
        gsSparseMatrix<T> m;
        GISMO_ENSURE(m_nonlinear(m_U, m),"Assembly failed");
        return stabilityVec(m);
    }

protected:

    gsVector<T> _solveLinear();
    gsVector<T> _solveNonlinear();

    using Base::_computeStability;

    void _init();
    void _start();
    void _factorizeMatrix(const gsSparseMatrix<T> & jacMat) const;
    gsVector<T> _solveSystem(const gsVector<T> & F);

    gsVector<T> _computeResidual(const gsVector<T> & U);

    gsSparseMatrix<T> _computeJacobian(const gsVector<T> & U, const gsVector<T> & deltaU);

protected:

    const gsSparseMatrix<T> & m_linear;
    const gsVector<T> & m_force;
    const Jacobian_t      m_nonlinear;
          dJacobian_t     m_dnonlinear;
          Residual_t      m_residualFun;
    const ALResidual_t    m_ALresidualFun;

    using Base::m_R;

    using Base::m_U;
    using Base::m_DeltaU;
    using Base::m_deltaU;

    using Base::m_L;

    using Base::m_stabilityVec;

    using Base::m_verbose;

    bool m_NL;

    using Base::m_options;

    using Base::m_tolF;
    using Base::m_tolU;

    mutable T m_relax;

    using Base::m_maxIterations;
    using Base::m_numIterations;

    // Residual norms
    using Base::m_residual;
    using Base::m_residualIni;
    using Base::m_residualOld;

    using Base::m_start;
    using Base::m_headstart;

    using Base::m_solver;   // Cholesky by default
    
    using Base::m_stabilityMethod;

    using Base::m_indicator;

    using Base::m_dofs;

    // Solver status
    using Base::m_status;
};


} // namespace gismo


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