/
LapEigen1DBeltrami.edp
39 lines (27 loc) · 1.09 KB
/
LapEigen1DBeltrami.edp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
// test to validate the addition of surfacic finite elements in FreeFEM
load "msh3"
load "medit"
/////variational form
/////////////////////////////////
// laplacian 2D
real R = 3, r=1;
real h = 0.1; //
int nx = R*2*pi/h;
func torex= (R+r*cos(y*pi*2))*cos(x*pi*2);
func torey= (R+r*cos(y*pi*2))*sin(x*pi*2);
func torez= r*sin(y*pi*2);
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++)
plot(eV[i],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);