// This FreeFEM script computes eigenvalues of the Laplace-Beltrami operator on the unit sphere
//
//---- A. DECLARATIONS ----
int npoints=30; 
int N=50; 
real[int] Evalues(N); 
real t; 
//---- END OF DECLARATIONS so far ----
//
//----  B. BOUNDARY DEFINITIONS ----
border dOmega1(t=0, 2*pi) { x = cos(t); y = 2+sin(t); label=1; };
border dOmega2(t=0, 2*pi) { x = cos(t); y = -2+sin(t); label=2; };
//---- END OF BOUNDARY DEFINITIONS ----
//
//---- C. MESH CREATION USING buildmesh COMMAND ----
mesh Th=buildmesh(dOmega1(npoints*2*pi)+dOmega2(npoints*2*pi));
//---- END OF MESH DEFINITIONS ----
//
//---- VISUALISE THE MESH, may be commented out
plot(Th,wait=1);
//
//---- D. DECLARE THE FEM SPACE AND FEM VARIABLES ----
fespace Vh(Th,P2);
Vh u,v;
//
//---- E. DECLARE THE QUADRATIC FORMS ----
varf q(u,v)=int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v) )+on(1,u=0);
varf b(u,v)=int2d(Th)(4/(1+x^2+(y-2*y/abs(y))^2)^2*u*v);
//
//----  F. CREATE THE MATRICES ----
matrix S=q(Vh,Vh);
matrix M=b(Vh,Vh);
//
//----  DECLARE THE ARRAY TO HOLD EIGENFUNCTIONS ----
Vh[int] Efunctions(N);
//
//---- G. SOLVE THE PROBLEM ----
int k=EigenValue(S,M,sym=true,value=Evalues,vector=Efunctions);
//
//---- H. PRINT THE EIGENVALUES ----
cout << Evalues;
//