/usr/share/doc/freefem++/examples/examples++-eigen/Lap3dEigenValue.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 | // Computation of the eigen value and eigen vector of the
// Dirichlet problem on cube $]0,\pi[^3$
// ----------------------------------------
// we use the inverse shift mode
// the shift is given with sigma real
// -------------------------------------
// find $\lamda$ such that:
// $$ \int_{\omega} \nabla u_ \nabla v = \lamba \int_{\omega} u \nabla v $$
load "msh3"
int nn=15;
mesh Th2=square(nn,nn,[pi*x,pi*y]);
fespace Vh2(Th2,P1);
int[int] rup=[0,1], rdown=[0,1], rmid=[4,1,2,1, 1,1 ,3,1];
real zmin=0,zmax=pi;
mesh3 Th=buildlayers(Th2,nn,
zbound=[zmin,zmax],
// region=r1,
labelmid=rmid,
reffaceup = rup,
reffacelow = rdown);
cout << "Th : nv = " << Th.nv << " nt =" << Th.nt << endl;
fespace Vh(Th,P1);
Vh u1,u2;
real sigma = 00; // value of the shift
macro Grad(u) [dx(u),dy(u),dz(u)] // EOM
varf a(u1,u2)= int3d(Th)( Grad(u1)'*Grad(u2) - sigma* u1*u2 ) //')
+ on(1,u1=0.) ; // Boundary condition
varf b([u1],[u2]) = int3d(Th)( u1*u2 ) ; // no Boundary condition
matrix A= a(Vh,Vh,solver=UMFPACK);
cout << " fin A .. " << endl;
matrix B= b(Vh,Vh,solver=CG,eps=1e-20);
// important remark:
// the boundary condition is make with exact penalisation:
// we put 1e30=tgv on the diagonal term of the lock degre of freedom.
// So take dirichlet boundary condition just on $a$ variationnal form
// and not on $b$ variationnanl form.
// because we solve
// $$ w=A^-1*B*v $$
int nev=10; // number of computed eigen valeu close to sigma
real[int] ev(nev); // to store nev eigein value
Vh[int] eV(nev); // to store nev eigen vector
int k=EigenValue(A,B,sym=true,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=0,ncv=0);
// tol= the tolerace
// maxit= the maximal iteration see arpack doc.
// ncv see arpack doc.
// the return value is number of converged eigen value.
k=min(k,nev); // some time the number of converged eigen value
// can be greater than nev;
int nerr=0;
real[int] eev(6*6*6);
eev=1e100;
for(int i=1,k=0;i<6;++i)
for(int j=1;j<6;++j)
for(int l=1;l<6;++l)
eev[k++]=i*i+j*j+l*l;
eev.sort;
cout << eev << endl;
for (int i=0;i<k;i++)
{
u1=eV[i];
real gg = int3d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) + dz(u1)*dz(u1) );
real mm= int3d(Th)(u1*u1) ;
real err = int3d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) + dz(u1)*dz(u1) - (ev[i])*u1*u1) ;
if(abs(err) > 1e-6) nerr++;
if(abs(ev[i]-eev[i]) > eev[i]*1e-1) nerr++;
cout << " ---- " << i<< " " << ev[i] << " == " << eev[i] << " err= " << err << " --- "<<endl;
plot(eV[i],cmm="Eigen 3d Vector "+i+" valeur =" + ev[i]+ " == " + eev[i] ,wait=1,value=1,ps="eigen"+i+".eps");
}
assert(nerr==0);
|