dune-pdelab  2.7-git
liptonbabuskabasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_LIPTONBABUSKABASIS_HH
2 #define DUNE_PDELAB_BACKEND_ISTL_GENEO_LIPTONBABUSKABASIS_HH
3 
4 #include <algorithm>
5 #include <functional>
6 
9 
10 #if HAVE_ARPACKPP
11 
12 namespace Dune {
13  namespace PDELab {
14 
19  template<class GFS, class M, class X, class Y, int dim>
20  class LiptonBabuskaBasis : public SubdomainBasis<X>
21  {
24 
25  public:
26  LiptonBabuskaBasis(const GFS& gfs, const M& AF_exterior, const M& AF_ovlp, const double eigenvalue_threshold, X& part_unity,
27  int& nev, int nev_arpack, double shift = 0.001, bool add_part_unity = false, int verbose = 0) {
29 
30  if (nev_arpack == -1)
31  nev_arpack = std::max(nev, 2);
32  if (nev_arpack < nev)
33  DUNE_THROW(Dune::Exception,"eigenvectors_compute is less then eigenvectors or not specified!");
34 
35  auto AF_interior = AF_exterior;
36  native(AF_interior) -= native(AF_ovlp);
37 
38  ArpackGeneo::ArPackPlusPlus_Algorithms<ISTLM, X> arpack(native(AF_exterior));
39  double eps = .0001;
40 
41  std::vector<double> eigenvalues(nev_arpack,0.0);
42  std::vector<X> eigenvectors(nev_arpack,X(gfs,0.0));
43 
44  arpack.computeGenNonSymMinMagnitude(native(AF_interior), eps, eigenvectors, eigenvalues, shift);
45 
46  // Count eigenvectors below threshold
47  int cnt = -1;
48  if (eigenvalue_threshold >= 0) {
49  for (int i = 0; i < nev; i++) {
50  if (eigenvalues[i] > eigenvalue_threshold) {
51  cnt = i;
52  break;
53  }
54  }
55  if (verbose > 0)
56  std::cout << "Process " << gfs.gridView().comm().rank() << " picked " << cnt << " eigenvectors" << std::endl;
57  if (cnt == -1)
58  DUNE_THROW(Dune::Exception,"No eigenvalue above threshold - not enough eigenvalues computed!");
59  } else {
60  cnt = nev;
61  }
62 
63  this->local_basis.resize(cnt);
64 
65  for (int base_id = 0; base_id < cnt; base_id++) {
66  this->local_basis[base_id] = std::make_shared<X>(part_unity);
67  // scale partition of unity with eigenvector
68  std::transform(
69  this->local_basis[base_id]->begin(),this->local_basis[base_id]->end(),
70  eigenvectors[base_id].begin(),
71  this->local_basis[base_id]->begin(),
72  std::multiplies<>()
73  );
74  }
75 
76  // Normalize basis vectors
77  for (auto& v : this->local_basis) {
78  *v *= 1.0 / v->two_norm2();
79  }
80 
81  if (add_part_unity && eigenvalues[0] > 1E-10) {
82  this->local_basis.insert (this->local_basis.begin(), std::make_shared<X>(part_unity));
83  this->local_basis.pop_back();
84  }
85  }
86 
87  };
88 
89  }
90 }
91 
92 #endif
93 
94 #endif //DUNE_PDELAB_BACKEND_ISTL_GENEO_LIPTONBABUSKABASIS_HH
arpackpp_geneo.hh
Dune::PDELab::Backend::Native
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper.
Definition: backend/interface.hh:176
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::Backend::native
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
subdomainbasis.hh