beameigenvalue.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

//  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.

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");
}