gsEigenProblemBase.h

Go to the documentation of this file.

#include <typeinfo>

#include <gsCore/gsLinearAlgebra.h>
#ifdef gsSpectra_ENABLED
#include <gsSpectra/gsSpectra.h>
#endif
#include <gsIO/gsOptionList.h>
#include <gsStructuralAnalysis/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h>

#pragma once


namespace gismo
{

template <class T>
class gsEigenProblemBase
{

public:

    gsEigenProblemBase()
    {
        m_options = defaultOptions();
        m_status  = gsStatus::NotStarted;
    }

    ~gsEigenProblemBase() {};

public:

    virtual gsOptionList defaultOptions()
    {
        gsOptionList options;
        options.addSwitch("verbose","Verbose output",false);
        options.addInt("solver","Spectra solver to be used (see Spectra documentation):"
                                    "0: Cholesky (default)"
                                    "1: RegularInverse"
                                    "2: ShiftInvert"
                                    "3: Buckling"
                                    "4: Cayley",0);

        options.addInt("selectionRule","Selection rule to be used (see Spectra documentation):"
                                        "0: LargestMagn"
                                        "1: LargestReal"
                                        "2: LargestImag"
                                        "3: LargestAlge"
                                        "4: SmallestMagn"
                                        "5: SmallestReal"
                                        "6: SmallestImag"
                                        "7: SmallestAlge"
                                        "8: BothEnds",4);

        options.addInt("sortRule","Sort rule to be used (see Spectra documentation):"
                                        "0: LargestMagn"
                                        "1: LargestReal"
                                        "2: LargestImag"
                                        "3: LargestAlge"
                                        "4: SmallestMagn"
                                        "5: SmallestReal"
                                        "6: SmallestImag"
                                        "7: SmallestAlge"
                                        "8: BothEnds",4);

        options.addInt("ncvFac","Factor for Spectra's ncv number. Ncv = ncvFac * numEigenvalues",3);
            options.addReal("tolerance","Tolerance for spectra and the power method",1e-10);
        options.addReal("shift","Shift for the eigenvalue solver",0.0);
        return options;
    }

    gsOptionList & options() {return m_options; };

    virtual void setOptions(gsOptionList & options) {m_options.update(options,gsOptionList::addIfUnknown); }

    virtual gsStatus compute();

    virtual gsStatus computeSparse(const index_t number = 10);

    virtual gsStatus computePower();

    virtual const gsMatrix<T> & values() const { return m_values; };
    virtual T value(int k) const { return m_values.at(k); };

    virtual const gsMatrix<T> & vectors() const { return m_vectors; };
    virtual gsMatrix<T> vector(int k) const { return m_vectors.col(k); };

    virtual std::vector<std::pair<T,gsMatrix<T>> > mode(int k) const {return makeMode(k); }

protected:

    virtual std::vector<std::pair<T,gsMatrix<T>> > makeMode(int k) const;

private:
    #ifdef gsSpectra_ENABLED
    template<Spectra::GEigsMode _GEigsMode>
    typename std::enable_if<_GEigsMode==Spectra::GEigsMode::Cholesky ||
                            _GEigsMode==Spectra::GEigsMode::RegularInverse
                            ,
                            gsStatus>::type computeSparse_impl(index_t number);
    #endif

    #ifdef gsSpectra_ENABLED
    template<Spectra::GEigsMode _GEigsMode>
    typename std::enable_if<_GEigsMode==Spectra::GEigsMode::ShiftInvert ||
                            _GEigsMode==Spectra::GEigsMode::Buckling ||
                            _GEigsMode==Spectra::GEigsMode::Cayley
                            ,
                            gsStatus>::type computeSparse_impl(index_t number);
    #endif

protected:

    gsSparseMatrix<T> m_A;
    gsSparseMatrix<T> m_B;

    gsOptionList m_options;

    gsEigen::GeneralizedSelfAdjointEigenSolver< typename gsMatrix<T>::Base >  m_eigSolver;

    gsMatrix<T> m_values,m_vectors;

    index_t m_num;

    gsStatus m_status;
};


} // namespace gismo


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