DUNE PDELab (git)

geneobasis.hh
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
7#include <dune/pdelab/backend/istl/geneo/subdomainbasis.hh>
8#include <dune/pdelab/backend/istl/geneo/arpackpp_geneo.hh>
9
10#if HAVE_ARPACKPP
11
12namespace Dune {
13 namespace PDELab {
14
21 template<class GFS, class M, class X, int dim>
22 class GenEOBasis : public SubdomainBasis<X>
23 {
24 typedef Dune::PDELab::Backend::Native<M> ISTLM;
25 typedef Dune::PDELab::Backend::Native<X> ISTLX;
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) {
44 using Dune::PDELab::Backend::native;
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
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)