DUNE PDELab (git)

liptonbabuskabasis.hh
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
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
19 template<class GFS, class M, class X, class Y, int dim>
20 class LiptonBabuskaBasis : public SubdomainBasis<X>
21 {
22 typedef Dune::PDELab::Backend::Native<M> ISTLM;
23 typedef Dune::PDELab::Backend::Native<X> ISTLX;
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) {
28 using Dune::PDELab::Backend::native;
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
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 13, 23:29, 2024)