// test to validate the addition of surfacis finite elements in FreeFEM

load "medit"


/////variational form

/////////////////////////////////
// laplacian 1D
real R = 3, r=1; 
real h = 0.01; // 
int nx = R*2*pi/h;
func torex= (R+r*cos(x*pi*10))*cos(x*pi*2);
func torey= (R+r*cos(x*pi*10))*sin(x*pi*2);
func torez= r*sin(x*pi*4);


meshL Th=segment(nx,[torex,torey,torez],removeduplicate=true) ;
fespace Vh(Th,P1);
macro Grad3(uSVar) [dx(uSVar),dy(uSVar),dz(uSVar)]  // EOM

//with variational form
real sigma = 1; 
varf a(u,v) = int1d(Th)(Grad3(u)'*Grad3(v)- sigma* u*v);
varf m(u,v) = int1d(Th)( u*v);

matrix A =a(Vh,Vh);
matrix B =m(Vh,Vh,solver=CG);
int nev=100;  // number of computed eigen valeu close to sigma

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

int k=EigenValue(A,B,sym=true,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=0,ncv=0,which="LM"); // which="LM" default 

for (int i=0;i<k;i++) 
{
	meshL Thm = movemesh(Th,[x,y,z+3*eV[i]]);
	fespace Wh(Thm,P1);
	Wh ew; ew[]=eV[i][];
	plot(Th,ew,Thm,cmm=" (new plot in test) Eigen  Vector "+i+" valeur =" + ev[i]  ,wait=1,value=1,dim=3,fill=1,CutPlane=0,ShowAxes=1, nbiso=40);  
}
