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
39
40
41
42
43
| real s0=clock();
int n=25;
mesh Th=square(n,n);
real[int,int] Bb((n+1)*(n+1)+2*(n-1)*(n-1)+8*n,(n+1)*(n+1)+2*(n-1)*(n-1)+8*n);
fespace Xh(Th,P2) ;
Xh u1,u2,v1,v2 ;
fespace Mh(Th,P1) ;
Mh p,q ;
fespace XhxXhxMh(Th,[P1,P1,P1]);
int i=0,iter=0;
solve stokes([u1,u2,p],[v1,v2,q]) = int2d(Th)((dx(u1)*dx(v1)+dy(u1)*dy(v1)
+ dx(u2)*dx(v2)+dy(u2)*dy(v2))
- p*q*(0.000001)
- p*(dx(v1)+dy(v2))
- q*(dx(u1)+dy(u2)))
+ on(1,2,4,u1=0,u2=0)
+ on(3,u1=1,u2=0) ;
plot([u1,u2],p,ps="stokes_exple.eps") ;
matrix A= vfstokes(XhxXhxMh,XhxXhxMh,solver=Crout,factorize=1);
Bb=0;
for(int i=0;i<(n+1)*(n+1)+2*(n-1)*(n-1)+8*n;++i)
Bb(i,i)=1;
matrix B;
B=Bb;
real sigma=0;
int nev=200; //nombre de valeurs propres autour de sigma
real[int] ev(nev); // stocker les 200 valeurs propres
XhxXhxMh [int] [eu1,eu2,ep](nev);// stocker les 200 vecteurs propres
int k=EigenValue(A,B,sym=true,sigma=0,value=ev,vector=eu1,tol=1e-10,maxit=0,ncv=0);
k=min(k,nev);
for(int i=0;i<k;i++)
{
cout<<"valeur propre:"<<i<<":"<<ev[i]<<endl;
u1=eu1[i];
u2=eu2[i];
p=ep[i];
plot([u1,u2],p,cmm="Eigen Vector "+i+"valeur="+ ev[i],wait=1,value=1,ps="stockes-eigen-p1-p1"+i+".eps");
}
cout <<"CPU"<<clock()-s0 <<"s"<<endl;
assert(abs(ev[0]-52.3471<0.1); |
