//  Computation of the eigen value and eigen vector of the 
// 
// ----------------------------------------
// we use the inverse shift mode 
// the shift is given with sigma real
// -------------------------------------
// Linear elastisitic problem in axisym. model. 
// thank to 
// http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch02.d/AFEM.Ch02.pdf
//  for the model.
//  the benchmark  case is
// thank form comsol and ansys 
// the comsol sol testcaes:  Eigenfrequency Analysis of a Free Cylinder
// The model is taken from NAFEMS Free Vibration Benchmarks
// 1. F. Abassian, D.J. Dawswell, and N.C. Knowles, Free Vibration Benchmarks, 
// Volume 3, NAFEMS, Glasgow, 1987
//  the comsol data: 
// http://www.comsol.com/model/download/40173/free_cylinder.pdf 
// F. Hecht may 2014 ??
// warning: Axi / 0x   
//  r == y,  
real[int] freqref=[0,243.53,377.41,394.11,397.72,405.28];   
verbosity=0;
real ccc= 1; // coef scale of matrix 
real E =  2.e11*ccc ;// 72e9;
real nu = 0.3;
real rho = 8000*ccc;//2700;//
real freq0 = 100; 
real mu = E/(2.*(1.+nu));
real lambda = E*nu/((1.+nu)*(1.-2.*nu));
real R1 = 1.8, R2 = R1+ 0.4, L = 10;
int nn = 40, nz = 40;
border cc1(t=0,2*pi){ x = R1*cos(t); y=R1*sin(t); }
border cc2(t=0,2*pi){ x = R2*cos(t); y=R2*sin(t); }
int NN = 10; 
mesh Th=square(L*NN,(R2-R1)*NN,[x*L, (1-y)*R1+y*R2]); 
plot(Th, wait=1); 


real shift = 1; // make a shift because 0 mode ... 
real C11=  2*mu  + lambda  ; 
real C12 =  lambda ;
real C44=  2*mu;
func Amat = [  [C11, C12, C12, 0.  ],
               [C12, C11, C12,  0. ],
               [C12, C12, C11,  0. ],
               [0. , 0. , 0.,   C44]
             ];

fespace Vh(Th,[P2,P2]);
Vh [u1,u2], [v1,v2];
func dr=dy; 
func r=y; 

macro Strainaxi(ux,ur)  [dx(ux),dr(ur),ur/r, dr(ux)+dx(ur) ] // EOM
  
varf vLame([u1,u2],[v1,v2])=
  int2d(Th)(  (Strainaxi(v1,v2)'*(Amat*Strainaxi(u1,u2))  
  - ([u1,u2]'*[v1,v2]*rho*shift) )*y )
  ;
varf vMass([u1,u2],[v1,v2])=
  int2d(Th)( ( [v1,v2]'*[u1,u2]*rho )*y );
 
matrix A= vLame(Vh,Vh,solver=sparsesolver); 
matrix B= vMass(Vh,Vh,solver=CG,eps=1e-20); 

int nev=6;  // number of computed eigen valeu close to sigma

real[int] ev(nev); // to store nev eigein value
Vh[int] [eu1,eu2](nev);   // to store nev eigen vector


int k=EigenValue(A,B,sigma=shift,sym=true,value=ev,vector=eu1,tol=1e-20);

k=min(k,nev); //  some time the number of converged eigen value 
              // can be greater than nev;
int nerr=0;
for (int i=0;i<k;i++)
{
    real e = ev[i] ;
    real freq = sqrt(abs(e))/(pi*2) ;
    if( abs(freq-freqref[i]) > 0.3) nerr++;
  cout << i << " ev " <<  e << " freq = " << freq  << "Hz, err=" << abs(freq-freqref[i]) <<endl;   
  real coef= 0.5/eu1[i][].linfty;
  mesh Thm=movemesh(Th,[x+eu1[i]*coef,y+eu2[i]*coef]);
  //plot(Th,Thm, wait=1,cmm=i+"/"+k+"coef  amplification = "+coef );
 // medit("Thm-"+ev[i],Thm,wait=1);
 plot(Thm,Th,wait=1,cmm=i+"/"+k+"coef  amplification = "+coef +" freq " + freq + "Hz"  , fill=1);
}
assert(nerr==0);
