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
|
assert(version>=1.24); // for big bug is non symetric matrix see HISTORY, and sign in int1d
include "beam.edp"
// Stokes on square b,e,f,g driven cavite on left side g
border e(t=0,10) { x=t; y=-10; label= 1; }; // bottom
border f(t=0,10) { x=10; y=-10+t ; label= 1; }; // right
border g(t=0,10) { x=0; y=-t ;label= 2;}; // left
border h(t=0,10) { x=t; y=vv(t,0)*( t>=0.001 )*(t <= 9.999); label=3;}; // top of cavity deforme
mesh sh = buildmesh(h(-20)+f(10)+e(10)+g(10));
plot(sh,wait=1);
fespace Xh(sh,P2),Mh(sh,P1);
fespace V1h(sh,P2);
Xh u1,u2,v1,v2;
Mh p,q,ppp;
varf bx(u1,q) = int2d(sh)( -(dx(u1)*q));
varf by(u1,q) = int2d(sh)( -(dy(u1)*q));
varf Lap(u1,u2)= int2d(sh)( dx(u1)*dx(u2) + dy(u1)*dy(u2) )
+ on(2,u1=1) + on(1,3,u1=0) ;
varf Mass(p,q)=int2d(sh)(p*q);
Xh bc1;
Xh brhs;
matrix A= Lap(Xh,Xh,solver=CG);
matrix M= Mass(Mh,Mh,solver=CG);
matrix Bx= bx(Xh,Mh);
matrix By= by(Xh,Mh);
Xh bcx,bcy;
func real[int] divup(real[int] & pp)
{ //
int verb=verbosity;
verbosity=1;
brhs[] = Bx'*pp; brhs[] += bc1[] .*bcx[];
u1[] = A^-1*brhs[];
brhs[] = By'*pp; brhs[] += bc1[] .*bcy[];
u2[] = A^-1*brhs[];
ppp[] = M^-1*pp;
ppp[] = 1.e-6*ppp[];
ppp[] = Bx*u1[];
ppp[] += By*u2[];
verbosity=verb;
return ppp[] ;
};
p=0;q=0;u1=0;v1=0;
int i;
for( i=0;i<2;i++)
{
bcx=0;
bcy= (-y)*(10.+y)/25.;
cout << " loop " << i << endl;
bc1[] = Lap(0,Xh);
q=0;
verbosity=50;
LinearCG(divup,p[],eps=1.e-3,nbiter=50);
divup(p[]);
verbosity=1;
plot([u1,u2],p,wait=1,value=true,coef=0.1,cmm="[u1,u2],p");
real coef=0.2;
V1h sigma11,sigma22,sigma12;
Vh [uu1,vv1];
[uu1,vv1]=[uu,vv];
sigma11([x+uu,y+vv]) = (2*dx(u1)-p);
sigma22([x+uu,y+vv]) = (2*dy(u2)-p);
sigma12([x+uu,y+vv]) = (dx(u2)+dy(u1));
solve bbst([uu,vv],[w,s],init=i) =
int2d(th)(
lambda*div(w,s)*div(uu,vv)
+2.*mu*( epsilon(w,s)'*epsilon(uu,vv) )
)
+ int2d(th) (-gravity*s)
+ int1d(th,bottombeam)( coef*(sigma11*N.x*w + sigma22*N.y*s + sigma12*(N.y*w+N.x*s) ) )
+ on(1,uu=0,vv=0)
;
cout << " max deplacement " << uu[].linfty << endl;
plot([uu,vv],wait=1);
real err = sqrt(int2d(th)( (uu-uu1)^2 + (vv-vv1)^2 ));
cout << " ----- Iteration = " << i << " Erreur L2 = " << err << "----------\n";
th1 = movemesh(th, [x+uu, y+vv]);
plot(th1,wait=1);
sh = buildmesh(h(-20)+f(10)+e(10)+g(10));
plot(sh);
p=p;q=q;u1=u1;u2=u2;brhs=brhs;ppp=ppp;
A= Lap(Xh,Xh,solver=CG);
M= Mass(Mh,Mh,solver=CG);
Bx= bx(Xh,Mh);
By= by(Xh,Mh);
bc1=0; // for resize
bc1[] = Lap(0,Xh);
}
cout << " max deplacement " << uu[].linfty << endl;
|
