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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
load "msh3"
newconvect=1;
real nu=0.01,dt=0.3;
real alpha=1./dt,alpha2=sqrt(alpha);
int nn=5;
mesh Th2=square(nn,nn);
fespace Vh2(Th2,P2);
Vh2 ux,uz,p2;
int[int] rup=[0,2], rdown=[0,1], rmid=[1,1,2,1,3,1,4,1];
real zmin=0,zmax=1;
mesh3 Th=buildlayers(Th2,nn,
zbound=[zmin,zmax],
// region=r1,
labelmid=rmid,
reffaceup = rup,
reffacelow = rdown);
fespace VVh(Th,[P23d,P23d,P23d,P13d]);
fespace Vh(Th,P23d);
fespace Ph(Th,P13d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM
varf vStokes([u1,u2,u3,p],[v1,v2,v3,q]) =
int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) + Grad(u2)'*Grad(v2) + Grad(u3)'*Grad(v3)
- div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(2,u1=1.,u2=0,u3= 0)
+ on(1,u1=0,u2=0,u3=0)
;
cout << "b mat " << endl;
matrix A=vStokes(VVh,VVh);
cout << "e mat " << endl;
set(A,solver=UMFPACK);
cout << "e fac mat " << endl;
real[int] b= vStokes(0,VVh);
VVh [u1,u2,u3,p];
VVh [X1,X2,X3,Xp];
VVh [x1,x2,x3,xp]=[x,y,z,0];
u1[]= A^-1 * b;
ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5",wait=1);
macro XX1() (x-u1*dt)//
macro XX2() (y-u2*dt)//
macro XX3() (z-u3*dt)//
varf vNS([uu1,uu2,uu3,p],[v1,v2,v3,q]) =
int3d(Th)( alpha*(uu1*v1+uu2*v2+uu3*v3) + nu*(Grad(uu1)'*Grad(v1) + Grad(uu2)'*Grad(v2) + Grad(uu3)'*Grad(v3))
- div(uu1,uu2,uu3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(2,uu1=1,uu2=0,uu3=0)
+ on(1,uu1=0,uu2=0,uu3=0)
+ int3d(Th,optimize=1,qforder=4)( alpha*( convect([u1,u2,u3],-dt,u1)*v1 + convect([u1,u2,u3],-dt,u2)*v2 + convect([u1,u2,u3],-dt,u3)*v3 ) ) ;
// + int3d(Th,optimize=1)( alpha*( u1(X1,X2,X3)*v1 + u2(X1,X2,X3)*v2 + u3(X1,X2,X3)*v3 ) ) ;
// + int3d(Th,optimize=1)( alpha*( u1(XX1,XX2,XX3)*v1 + u2(XX1,XX2,XX3)*v2 + u3(XX1,XX2,XX3)*v3 ) ) ;
//+ int3d(Th,optimize=1)( alpha*( u1(x,y,z)*v1 + u2(x,y,z)*v2 + u3(x,y,z)*v3 ) ) ;
//+ int3d(Th,optimize=1)( alpha*( u1*v1 + u2*v2 + u3*v3 ) ) ;
cout << " build A" << endl;
A = vNS(VVh,VVh);
cout << " fac A" << endl;
set(A,solver=UMFPACK);
real t=0;
for(int i=0;i<10;++i)
{
t += dt;
cout << " iteration " << i << " t = " << t << endl;
X1[]=x1[]+u1[]*(-dt);
// verbosity=1000;
b=vNS(0,VVh);
verbosity=2;
u1[]= A^-1 * b;
ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=0);
if(i%5==6)
{
exec("mkdir dd");
string prefu="dd/u-"+(100+i);
string prefp="dd/p-"+(100+i);
savemesh(Th,prefu+".mesh");
savemesh(Th,prefp+".mesh");
ofstream file(prefu+".bb");
ofstream filep(prefp+".bb");
Ph up1=u1,up2=u2,up3=u3,pp=p;
file << "3 1 3 "<< up1[].n << " 2 \n";
filep << "3 1 1 "<< pp[].n << " 2 \n";
for (int j=0;j<up1[].n ; j++)
{
file << up1[][j] <<" " <<up2[][j] <<" "<< up3[][j] <<"\n";
filep << pp[][j] << endl;
}
}
}
plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=1);
|
