liptonbabuskabasis.hh

 
#ifndef DUNE_GENEO_LIPTONBABUSKA_HH
#define DUNE_GENEO_LIPTONBABUSKA_HH
 
#include "subdomainbasis.hh"
#include "arpackpp_fork.hh"
 
namespace Dune{
    namespace Geneo{
 
template<class GFS, class M, class X, class Y, int dim>
class LBBasis : public SubdomainBasis<X>
{
  typedef Dune::PDELab::Backend::Native<M> ISTLM;
  typedef Dune::PDELab::Backend::Native<X> ISTLX;
 
public:
  LBBasis(const GFS& gfs, const M& AF_neumann, const M& AF_ovlp, const double eigenvalue_threshold, X& part_unity, int& nev, int nev_arpack, double shift) {
    using Dune::PDELab::Backend::native;
 
    auto AF_interior = AF_neumann;
    native(AF_interior) -= native(AF_ovlp);
 
    ArpackFork::ArPackPlusPlus_Algorithms<ISTLM, ISTLX> arpack(native(AF_neumann));
    double eps = .0001;
 
    //if (nev_arpack < nev)
    //  DUNE_THROW(Dune::ISTLError,"eigenvectors_compute is less then eigenvectors or not specified!");
 
    std::vector<double> eigenvalues;
    std::vector<ISTLX> eigenvectors;
    eigenvectors.resize(nev_arpack);
    for (int i = 0; i < nev_arpack; i++) {
      eigenvectors[i] = native(X(gfs,0.0));
    }
 
    arpack.computeGenNonSymMinMagnitude(native(AF_interior), eps, nev_arpack, eigenvectors, eigenvalues, shift);
 
 
    /*for(int i = 0; i < nev_arpack; i++){
        auto check = eigenvectors[i];
        native(AF_neumann).mv(eigenvectors[i],check);
        auto check2 = eigenvectors[i];
        native(ovlp_mat).mv(eigenvectors[i],check2);
        check2 *= eigenvalues[i];
        check -= check2;
        if(native(check).infinity_norm() > 1e-6) std::cout << "Rank " << gfs.gridView().comm().rank() << " Error in EV calculation " << native(check).infinity_norm() << std::endl;
    }*/
 
 
    // Count eigenvectors below threshold
    int cnt = -1;
    if (configuration.get<bool>("threshold_eigenvectors",false)) {
      for (int i = 0; i < nev; i++) {
        if (eigenvalues[i] > eigenvalue_threshold) {
          cnt = i;
          break;
        }
      }
      if (configuration.get<bool>("verbose",false))
        std::cout << "Process " << gfs.gridView().comm().rank() << " picked " << cnt << " eigenvectors" << std::endl;
      if (cnt == -1)
        DUNE_THROW(Dune::ISTLError,"No eigenvalue above threshold - not enough eigenvalue solves!");
    } else {
      cnt = nev;
    }
 
    this->local_basis.resize(cnt);
 
    for (int base_id = 0; base_id < cnt; base_id++) {
      this->local_basis[base_id] = std::make_shared<X>(part_unity);
      for (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 (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();
    }
  }
 
};
 
}
}
 
#endif //DUNE_GENEO_GENEOBASIS_HH