計算部分
real area= int2d(Th)(1.); //面積全体Ω
Vh uold = u, vold = v, pold=p; //流速u,v圧力p
Vh f=convect([u,v],-dt,uold), g=convect([u,v],-dt,vold); //移流項f,gとする。
solve pb4u(u,w,init=n,solver=LU)
=int2d(Th)((u-f)*w/dt + dx(p)*w +nu*(dx(u)*dx(w)+dy(u)*dy(w)))
+ on(1,u = 4*y*(1-y)) + on(2,4,u = 0)+ on(3,u=f); //境界条件
solve pb4v(v,w,init=n,solver=LU)
= int2d(Th)((v-g)*w/dt + dy(p)*w +nu*(dx(v)*dx(w)+dy(v)*dy(w)))
+on(1,2,3,4,v = 0); //境界条件
real meandiv = int2d(Th)(dx(u)+dy(v))/area; //uのdivの平均
solve pb4p(q,w,init=n,solver=LU)= int2d(Th)(dx(q)*dx(w)+dy(q)*dy(w))
- int2d(Th)((dx(u)+ dy(v)-meandiv)*w/dt)
+ on(3,q=0); //境界条件
real meanpq = int2d(Th)(pold - q)/area; //pold-qの平均
p = pold-q-meanpq;
u = u + dx(q)*dt;
v = v + dy(q)*dt;