geneobasis.hh

#ifndef DUNE_GENEO_GENEOBASIS_HH
#define DUNE_GENEO_GENEOBASIS_HH
 
#include "subdomainbasis.hh"
#include "arpackpp_fork.hh"
 
namespace Dune{
    namespace Geneo{
 
        template<class GFS, class M, class X, class Y, class SI, int dim>
            class GenEOBasis : public SubdomainBasis<X>
        {
            typedef Dune::PDELab::Backend::Native<M> ISTLM;
            typedef Dune::PDELab::Backend::Native<X> ISTLX;
 
            public:
            GenEOBasis(const GFS& gfs, const M& AF_neumann, const M& AF_ovlp, const double eigenvalue_threshold, X& part_unity, SI& solver_info){//, SubdomainBasis<X>& localbasis2){
                using Dune::PDELab::Backend::native;
                M ovlp_mat(AF_ovlp);
                // X * A_0 * X
                for (auto row_iter = native(ovlp_mat).begin(); row_iter != native(ovlp_mat).end(); row_iter++) {
                    for (auto col_iter = row_iter->begin(); col_iter != row_iter->end(); col_iter++) {
                        *col_iter *= native(part_unity)[row_iter.index()] * native(part_unity)[col_iter.index()];
                    }
                }
 
                ArpackFork::ArPackPlusPlus_Algorithms<ISTLM, ISTLX> arpack(native(AF_neumann));
                double eps = 0.0001;
 
                if (solver_info.nev_arpack < solver_info.nev)
                  DUNE_THROW(Dune::ISTLError,"eigenvectors_compute is less then eigenvectors or not specified!");
 
                std::vector<double> eigenvalues;
                std::vector<ISTLX> eigenvectors;
                eigenvectors.resize(solver_info.nev_arpack);
                for (int i = 0; i < solver_info.nev_arpack; i++) {
                    eigenvectors[i] = native(X(gfs,0.0));
                }
 
                if (configuration.get<bool>("EigenvectorVTK",false)) {
                    for (int i = 0; i < solver_info.nev; i++) {
                        auto gv = gfs.gridView();
                        typedef typename GFS::Traits::GridView GV;
                        Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,0);
                        X eigenvector(gfs, eigenvectors[i]);
                        Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,eigenvector);
                        vtkwriter.write("eigenvector"+std::to_string(i),Dune::VTK::appendedraw);
                    }
                }
 
                arpack.computeGenNonSymMinMagnitude(native(ovlp_mat), eps, solver_info.nev_arpack, eigenvectors, eigenvalues, solver_info.eigen_shift);
 
                // Count eigenvectors below threshold
                int cnt = -1; //int cnt2 = -1;
                if (eigenvalue_threshold != 0) {
                    /*for (int i = 0; i < nev; i++) {
                        if (eigenvalues[i] > eigenvalue_threshold/configuration.get<double>("prec_threshhold",2)) {
                            cnt2 = i;
                            break;
                        }
                    }*/
                    for (int i = 0; i < solver_info.nev; i++) {
                        if (eigenvalues[i] > eigenvalue_threshold) {
                            cnt = i;
                            break;
                        }
                    }
                    std::cout << "Process " << gfs.gridView().comm().rank() << " picked " << cnt << " eigenvectors" << std::endl;
 
                    if (cnt == -1){
                        std:: cout << "No eigenvalue above threshold - not enough eigenvalue solves!" << std::endl;
                        cnt  = solver_info.nev;
                    }
                    /*if (cnt2 == -1){
                        cnt2 = nev/configuration.get<double>("prec_threshhold",2);
                    }*/
                }
                else {
                    cnt = solver_info.nev;
                    //cnt2 = nev;
                }
 
                /*X numEV(gfs,cnt);
                auto gv = gfs.gridView();
                typedef typename GFS::Traits::GridView GV;
                Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,0);
                Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,numEV);
                vtkwriter.write("numEV",Dune::VTK::appendedraw);*/
 
 
                this->local_basis.resize(cnt);
                //localbasis2.local_basis.resize(cnt2);
 
                for (int base_id = 0; base_id < cnt; base_id++) {
                    this->local_basis[base_id] = std::make_shared<X>(part_unity);
                    //if(base_id < cnt2) localbasis2.local_basis[base_id] = std::make_shared<X>(part_unity);
                    for (unsigned int it = 0; it < native(eigenvectors[base_id]).N(); it++) {
                        for(int j = 0; j < dim; j++){
                            native(*this->local_basis[base_id])[it][j] *= eigenvectors[base_id][it][j];
                            //if(base_id < cnt2) native(*localbasis2.local_basis[base_id])[it][j] *= eigenvectors[base_id][it][j];
                        }
                    }
                }
 
                if (configuration.get<bool>("AddPartUnityToEigenvectors",false) && eigenvalues[0] > 1E-10) {
                    auto part_unity_basis = std::make_shared<X>(gfs,0.0);
                    *part_unity_basis = part_unity;
                    this->local_basis.insert (this->local_basis.begin(), part_unity_basis);
                    this->local_basis.pop_back();
                    //localbasis2.local_basis.insert (localbasis2.local_basis.begin(), part_unity_basis);
                    //localbasis2.local_basis.pop_back();
                }
            }
 
        };
 
    }
}
 
#endif //DUNE_GENEO_GENEOBASIS_HH