// MAP 562 (G. Allaire, T. Wick)
// Winter 2016/2017
//
// TD 2
// Exercise 1

// Task 1: Change the domain from square to rectangle
// Task 2: Change the boundary condition such that 
//   a part is a Neumann condition

// Number of boundary edges (as before)
int N=20;

// We declare a mesh
mesh Th;
Th=square(N,N);

// Numbering of the boundary parts
int[int] labels=[1,2,3,4];

// We associate a finite element space 
// with linear test functions
fespace Vh(Th,P1);

// Declaring trial and test functions
Vh u,v;

// Initializing of the eigenvalue 
real sigma = 0;

macro grad(u) [dx(u),dy(u)] //	gradient operator

// We realize the PDE as given in the problem statement
varf laplaceeigen(u,v)= int2d(Th)(grad(u)'*grad(v) - sigma* u*v) + on(labels,u=0) ;

// For later, namely calling `EigenValue' we need to 
// create a mass matrix of the right hand side
varf b(u,v) = int2d(Th)( u*v ) ;
 
// For solving the main problem, we simply 
// can use a CG method (the problem is symmetric!)
matrix M= laplaceeigen(Vh,Vh,solver=CG); 

// The second problem is again symmetric
matrix B= b(Vh,Vh,solver=CG);

// We initialize variables for 
// the number of eigenvalues, the eigenvalue,
// and the corresponding eigenfunction
int numberev=1;
real[int] eigv(numberev);
Vh[int] eigFunction(numberev); 

// This function is implemented in Arpack (external package)
EigenValue(M,B,sym=true,value=eigv,vector=eigFunction);

// Print the eigenvalue we have computed
cout<<"ev "<<eigv<<endl;

// Print (plot) the corresponding eigenfunction
plot(eigFunction);