gsBucklingSolver.h

Go to the documentation of this file.

 
#include <gsStructuralAnalysis/src/gsEigenSolvers/gsEigenProblemBase.h>
 
#include <gsIO/gsFileData.h>
#include <gsIO/gsReadFile.h>
 
#pragma once
 
namespace gismo
{
 
template <class T>
class gsBucklingSolver : public gsEigenProblemBase<T>
{
protected:
 
    typedef gsEigenProblemBase<T> Base;
 
    typedef std::function < bool ( gsVector<T> const &,                         gsSparseMatrix<T> & ) > Jacobian_t;
    typedef std::function < bool ( gsVector<T> const &, gsVector<T> const &,    gsSparseMatrix<T> & ) > dJacobian_t;
 
public:
 
    gsBucklingSolver(   gsSparseMatrix<T> &linear,
                        gsVector<T> &rhs,
                        Jacobian_t  &nonlinear,
                        T scaling = 1.0) :
    m_rhs(rhs),
    m_nonlinear(nonlinear),
    m_scaling(scaling)
    {
        m_A = linear;
        m_dnonlinear = [this](gsVector<T> const & x, gsVector<T> const & /*dx*/, gsSparseMatrix<T> & m) -> bool
        {
            return m_nonlinear(x,m);
        };
 
        m_solver = gsSparseSolver<T>::get( "SimplicialLDLT" );
 
        m_status = this->initializeMatrix();
    }
 
  gsBucklingSolver(     gsSparseMatrix<T> &linear,
                        gsVector<T> &rhs,
                        dJacobian_t &dnonlinear,
                        T scaling = 1.0) :
    m_rhs(rhs),
    m_dnonlinear(dnonlinear),
    m_scaling(scaling)
    {
        m_A = linear;
 
        m_solver = gsSparseSolver<T>::get( "SimplicialLDLT" );
 
        m_status = this->initializeMatrix();
    }
 
 
    gsBucklingSolver(     gsSparseMatrix<T> &linear,
                        gsSparseMatrix<T> &nonlinear )
    {
        m_A = linear;
        m_B = nonlinear-m_A;
    }
 
  // todo: add solver option
  // gsOptionList defaultOptions()
  // {
  //   gsOptionList options;
  //   options = Base::defaultOptions()
  //   options.addString("Solver","Specify the sparse solver","SimplicialLDLT");
  //   return options;
  // }
 
    gsStatus computeSparse(const index_t number = 10)
    {
        if (m_options.getInt("solver") != 3 )
            gsWarn<<"It is highly recommended to use the Buckling solver (option 4)\n";
 
        return Base::computeSparse(number);
    }
 
 
protected:
 
    gsStatus initializeMatrix()
    {
        bool verbose = m_options.getSwitch("verbose");
        if (verbose) { gsInfo<<"Computing matrices" ; }
        m_solver->compute(m_A);
        if (verbose) { gsInfo<<"." ; }
        m_solVec = m_solver->solve(m_scaling*m_rhs);
        if (verbose) { gsInfo<<"." ; }
        try
        {
            m_dnonlinear(m_solVec,gsVector<T>::Zero(m_solVec.rows()),m_B);
        }
        catch (...)
        {
            m_status = gsStatus::AssemblyError;
            return m_status;
        }
        m_B -= m_A;
        if (verbose) { gsInfo<<"." ; }
        if (verbose) { gsInfo<<"Finished\n" ; }
 
        return m_status;
    }
 
protected:
 
    using Base::m_A;
    const gsVector<T> m_rhs;
    const Jacobian_t m_nonlinear;
         dJacobian_t m_dnonlinear;
    T m_scaling;
    using Base::m_B;
 
 
    mutable typename gsSparseSolver<T>::uPtr m_solver;
    gsVector<> m_solVec;
 
    using Base::m_options;
 
    using Base::m_status;
};
 
} // namespace gismo