1//  Computation of the eigen value and eigen vector of the
2// Dirichlet problem  on square $]0,\pi[^2$
3// Functionnal interface
4// ----------------------------------------
5// we use the inverse shift mode
6// the shift is given with sigma real
7// -------------------------------------
8//  find $\lamda$ such that:
9// $$  \int_{\omega}  \nabla u_ \nabla v = \lamba \int_{\omega} u \nabla v  $$
10verbosity=1;
11mesh Th=square(20,20,[pi*y,pi*x]);
12fespace Vh(Th,P2);
13Vh u1,u2;
14int n=Vh.ndof;
15
16real sigma = 00;  // value of the shift
17
18varf  a(u1,u2)= int2d(Th)(  dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma* u1*u2 )
19                    + int1d(Th)(u1*u2) ;  // Boundary condition
20
21varf b([u1],[u2]) = int2d(Th)(  u1*u2 ) ; // no  Boundary condition
22
23matrix A= a(Vh,Vh,solver=sparsesolver);
24matrix B= b(Vh,Vh,solver=sparsesolver);
25
26func real[int] FA1(real[int] & u) { real[int] Au=A^-1*u;return Au;}
27func real[int] FB(real[int] & u) { real[int] Au=B*u;return Au;}
28func real[int] FA(real[int] & u) { real[int] Au=A*u;return Au;}
29func real[int] FB1(real[int] & u) { real[int] Au=B^-1*u;return Au;}
30// important remark:
31// the boundary condition is make with exact penalisation:
32//     we put 1e30=tgv  on the diagonal term of the lock degre of freedom.
33//  So take dirichlet boundary condition just on $a$ variationnal form
34// and not on  $b$ variationnanl form.
35// because we solve
36//  $$ w=A^-1*B*v $$
37
38int nev=20;  // number of computed eigen valeu close to sigma
39
40real[int] ev(nev); // to store nev eigein value
41Vh[int] eV(nev);   // to store nev eigen vector
42
43
44int k=EigenValue(n,A1=FA1,B=FB,A=FA,B1=FB1,sym=true,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=1000,ncv=200
45	,mode=3,which="LM");
46//   tol= the tolerace
47//   maxit= the maximal iteration see arpack doc.
48//   ncv   see arpack doc.
49//  the return value is number of converged eigen value.
50k=min(k,nev); //  some time the number of converged eigen value
51              // can be greater than nev;
52int nerr=0;
53for (int i=0;i<k;i++)
54{
55  u1=eV[i];
56  real gg = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1)) + int1d(Th)(u1*u1);
57  real mm= int2d(Th)(u1*u1) ;
58  real err = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) - (ev[i])*u1*u1) +  int1d(Th)(u1*u1) ;
59  if(abs(err) > 1e-6) nerr++;
60  cout << " ---- " <<  i<< " " << ev[i]  << " err= " << err << " --- "<<endl;
61
62  // FFCS: add 3D view capabilities
63  plot(eV[i],cmm="Eigen  Vector "+i+" valeur =" + ev[i]  ,wait=0,value=1,dim=3,fill=1,CutPlane=0,ShowAxes=0);
64}
65
66