mesh Th=square(50,50,[x*pi,y*pi]);
macro grad(u) [dx(u),dy(u)]//
int[int] L=[1,2,3,4];
varf vlapWithBC(u,v)= int2d(Th)(grad(u)'*grad(v))+ on(L,u=0);
varf vlapWithoutBC(u,v)= int2d(Th)(grad(u)'*grad(v));
varf vMass(u,v)= int2d(Th)(u*v);
fespace Vh(Th,P1);
Vh u,v;
u=1; 
matrix Id=u[]; // vMass(Vh,Vh,solver=CG) ;
set(Id,solver=CG,verb=0);
matrix Abc= vlapWithBC(Vh,Vh,solver=sparsesolver);
matrix Ancl= vlapWithoutBC(Vh,Vh,solver=CG);
real[int] aevmin(2),aevmax(2);
int kmin=EigenValue(Abc,Id,sym=true,value=aevmin);
int kmax=EigenValue(Id,Ancl,sym=true,value=aevmax);
real evmin = aevmin(0);
real evmax = 1./aevmax(0);
cout << "VP min = " << evmin << "VP max = " << evmax << " cond number = " << evmax/evmin << endl; 
