/usr/share/doc/freefem++/examples/examples++-eigen/BeamEigenValue.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 | // Computation of the eigen value and eigen vector of the
// Dirichlet problem on square $]0,\pi[^2$
// ----------------------------------------
// 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 $$
verbosity=1;
int bottombeam = 2;
border aaa(t=2,0) { x=0; y=t ;label=1;}; // left beam
border bbb(t=0,10) { x=t; y=0 ;label=bottombeam;}; // bottom of beam
border ccc(t=0,2) { x=10; y=t ;label=1;}; // rigth beam
border ddd(t=0,10) { x=10-t; y=2; label=3;}; // top beam
real E = 21.5;
real sigma = 0.29;
real mu = E/(2*(1+sigma));
real lambda = E*sigma/((1+sigma)*(1-2*sigma));
real gravity = -0.05;
mesh Th = buildmesh( bbb(20)+ccc(5)+ddd(20)+aaa(5));
fespace Vh(Th,[P1,P1]);
Vh [uu,vv], [w,s];
cout << "lambda,mu,gravity ="<<lambda<< " " << mu << " " << gravity << endl;
// deformation of a beam under its own weight
real shift = 0; // value of the shift
varf a([uu,vv],[w,s])=
int2d(Th)(
2*mu*(dx(uu)*dx(w)+dy(vv)*dy(s)+ ((dx(vv)+dy(uu))*(dx(s)+dy(w)))/2 )
+ lambda*(dx(uu)+dy(vv))*(dx(w)+dy(s))
- shift* (uu*w + vv*s)
)
+ on(1,uu=0,vv=0)
;
varf b([uu,vv],[w,s])=
int2d(Th)(uu*w + vv*s) ;
matrix A= a(Vh,Vh,solver=Crout,factorize=1);
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=20; // number of computed eigen valeu close to sigma
real[int] ev(nev); // to store nev eigein value
Vh[int] [eV,eW](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;
for (int i=0;i<k;i++)
{
[uu,vv]=[eV[i],eW[i]];
plot([uu,vv],cmm="Eigen Vector "+i+" valeur =" + ev[i] ,wait=1,value=1,ps="eigen"+i+".eps");
}
|