// Computation of Steklov or Wentzell eigenvalues
// FreeFem implementation by Beniamin Bogosel
// 
// uses radial representation for the shape in terms of Fourier coefficients
// parametrization uses a vector of the form
// vec = [a0, as, bs] where as are coefficients of cos and bs coefficients of sin
//
// The code shows how to implement an eigenvalue problem
//     and how to handle tangential derivatives in weak form
//
// Please read and refer to the associated paper
// B. Bogosel, The method of fundamental solutions applied to boundary eigenvalue problems, 2016, Journal of Computational and Applied Mathematics
//

real beta = 0; // parameter for the Laplace-Beltrami coefficient 
                 // For Steklov eigenvalue choose 0    
real[int] vec;
int M = 300;
int eigCount = 11; // the number of eigenvalues to compute

vec = [1.0,0,0,0,0.3];  // vector of Fourier coefficients for the 
                        // boundary parametrization
vec = [1,0.1,0,0,0,0.1,0,0.1,0,0,-0.1];

vec = [1,0.1,0.2,0.3,-0.1,0.1,0.4];

// evaluation of a trigonometric polynomial
func real ptrig(real t, real[int] VV){
  real[int] vect = VV;
  int n = (vect.n-1)/2;
  real[int] as = vect(1:n) ;    // coeffs of cos
  real[int] bs = vect(n+1:2*n); // coeffs of sin
  real a0 =vect(0);
  real Sum = a0;   // initialize sum
  for(int i=0;i<n;i++){Sum+=(as[i]*cos((i+1)*t)+bs[i]*sin((i+1)*t));}
  return Sum;
}

// Construction of the mesh
border C(t=0,2*pi){x=cos(t)*ptrig(t,vec); 
                   y=sin(t)*ptrig(t,vec);
                   label=1;}
mesh Th = buildmesh (C(200));
 Th = adaptmesh(Th,0.05,IsMetric=1,nbvx=20000);
 Th = adaptmesh(Th,0.02,IsMetric=1,nbvx=30000);
plot(Th);

load "Element_P3";
fespace Vh(Th,P2);
Vh uh,vh;
// weak form of the Wentzell problem
varf va(uh, vh) = int2d(Th)( dx(uh)*dx(vh)+dy(uh)*dy(vh))+
                  int1d(Th,1)( beta*(dx(uh)*dx(vh)-
       dx(uh)*N.x*(N.x*dx(vh)+N.y*dy(vh))-
       dx(vh)*N.x*(N.x*dx(uh)+N.y*dy(uh))+
       N.x*(dx(vh)*N.x+dy(vh)*N.y)*N.x*(dx(uh)*N.x+dy(uh)*N.y)+
       dy(uh)*dy(vh)-
       dy(uh)*N.y*(dx(vh)*N.x+dy(vh)*N.y)-
       dy(vh)*N.y*(dx(uh)*N.x+dy(uh)*N.y)+
       (N.y)^2*(dx(vh)*N.x+dy(vh)*N.y)*(dx(uh)*N.x+dy(uh)*N.y)));

varf vb(uh, vh) = int1d(Th,1)( uh * vh); // linear form
matrix A = va(Vh, Vh ,solver = sparsesolver); // Matrix A on left side
matrix B = vb(Vh, Vh);                        // Matrix B on right side
real cpu=clock();  // get the clock in seconds

// Get first Eigenvalues
real[int] ev(eigCount); // Holds Eigenvalues
Vh[int] eV(eigCount);   // Holds Eigenfunctions

int numEigs = EigenValue(A,B,sym=true,sigma=0,value=ev,vector=eV);
for(int i=0;i<eigCount;i++) { // Plot the spectrum and show the eigenvalues
plot(eV[i],fill=true,value=true,cmm= ev[i]);
 cout << "Eigenvalue " << i << " :"  << ev[i] << endl;
}
cout << " Total CPU time = " << clock()-cpu << endl;