dune-pdelab  2.7-git
geneobasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_PDELAB_BACKEND_ISTL_GENEO_GENEOBASIS_HH
2 #define DUNE_PDELAB_BACKEND_ISTL_GENEO_GENEOBASIS_HH
3 
4 #include <algorithm>
5 #include <functional>
6 
9 
10 #if HAVE_ARPACKPP
11 
12 namespace Dune {
13  namespace PDELab {
14 
21  template<class GFS, class M, class X, int dim>
22  class GenEOBasis : public SubdomainBasis<X>
23  {
26 
27  public:
28 
42  GenEOBasis(const GFS& gfs, const M& AF_exterior, const M& AF_ovlp, const double eigenvalue_threshold, X& part_unity,
43  int& nev, int nev_arpack = -1, double shift = 0.001, bool add_part_unity = false, int verbose = 0) {
45 
46  if (nev_arpack == -1)
47  nev_arpack = std::max(nev, 2);
48  if (nev_arpack < nev)
49  DUNE_THROW(Dune::Exception,"nev_arpack is less then nev!");
50 
51  const int block_size = X::block_type::dimension;
52 
53  // X * A_0 * X
54  M ovlp_mat(AF_ovlp);
55  for (auto row_iter = native(ovlp_mat).begin(); row_iter != native(ovlp_mat).end(); row_iter++) {
56  for (auto col_iter = row_iter->begin(); col_iter != row_iter->end(); col_iter++) {
57  for (int block_row = 0; block_row < block_size; block_row++) {
58  for (int block_col = 0; block_col < block_size; block_col++) {
59  (*col_iter)[block_row][block_col] *= native(part_unity)[row_iter.index()][block_row] * native(part_unity)[col_iter.index()][block_col];
60  }
61  }
62  }
63  }
64 
65  // Setup Arpack for solving generalized eigenproblem
66  ArpackGeneo::ArPackPlusPlus_Algorithms<ISTLM, X> arpack(native(AF_exterior));
67  double eps = 0.0;
68 
69  std::vector<double> eigenvalues(nev_arpack,0.0);
70  std::vector<X> eigenvectors(nev_arpack,X(gfs,0.0));
71 
72  arpack.computeGenNonSymMinMagnitude(native(ovlp_mat), eps, eigenvectors, eigenvalues, shift);
73 
74  // Count eigenvectors below threshold
75  int cnt = -1;
76  if (eigenvalue_threshold >= 0) {
77  for (int i = 0; i < nev; i++) {
78  if (eigenvalues[i] > eigenvalue_threshold) {
79  cnt = i;
80  break;
81  }
82  }
83  if (verbose > 0)
84  std::cout << "Process " << gfs.gridView().comm().rank() << " picked " << cnt << " eigenvectors" << std::endl;
85  if (cnt == -1)
86  DUNE_THROW(Dune::Exception,"No eigenvalue above threshold - not enough eigenvalues computed!");
87  } else {
88  cnt = nev;
89  }
90 
91  // Write results to basis
92  this->local_basis.resize(cnt);
93  for (int base_id = 0; base_id < cnt; base_id++) {
94  this->local_basis[base_id] = std::make_shared<X>(part_unity);
95  // scale partition of unity with eigenvector
96  std::transform(
97  this->local_basis[base_id]->begin(),this->local_basis[base_id]->end(),
98  eigenvectors[base_id].begin(),
99  this->local_basis[base_id]->begin(),
100  std::multiplies<>()
101  );
102  }
103 
104  // Normalize basis vectors
105  for (auto& v : this->local_basis) {
106  *v *= 1.0 / v->two_norm2();
107  }
108 
109  // Optionally add partition of unity to eigenvectors
110  // Only if there is no near-zero eigenvalue (that usually already corresponds to a partition of unity!)
111  if (add_part_unity && eigenvalues[0] > 1E-10) {
112  this->local_basis.insert (this->local_basis.begin(), std::make_shared<X>(part_unity));
113  this->local_basis.pop_back();
114  }
115  }
116  };
117 
118 
119  }
120 }
121 
122 #endif
123 
124 #endif //DUNE_PDELAB_BACKEND_ISTL_GENEO_GENEOBASIS_HH
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
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