// math formulation. // inlet label 4, outlet label 2 int inlet =2, outlet = 4; real R = 1,L=1./2.*pi/R, l=1 ; int n = 10; mesh Th = square(L*n*R,l*n,[x*L,y*l],flags=3); //plot(Th,wait=1); Th = movemesh(Th,[sin(x)*(R+(1-y)),-cos(x)*(R+(1-y))+R]);// y linear , x circle x=0 -> x =0 plot(Th,wait=1); real fin= 1 ;// normal Force inlet real fout = -1;// normal Force outlet fespace Vh(Th,P2); fespace Ph(Th,P1); macro grad(u) [dx(u),dy(u)] // macro Grad(u1,u2) [ grad(u1), grad(u2)]// macro epsilon(u1,u2) [dx(u1), dy(u2), (dx(u2)+dy(u1))*0.5, (dx(u2)+dy(u1))*0.5] // macro div(u1,u2) (dx(u1)+dy(u2)) // Vh u1=0,u2=0; Vh v1,v2; Ph p,q; real mu =1; verbosity=1; solve StokesMath([ u1,u2,p], [v1,v2,q] ) = int2d(Th) ( mu *( Grad(v1,v2):Grad(u1,u2) ) - p * div(v1,v2) - q *div(u1,u2) ) - int1d(Th,inlet)( [v1,v2]'*[N.x,N.y]*fin) - int1d(Th,outlet)( [v1,v2]'*[N.x,N.y]*fout) + on( 1,3,u1=0,u2=0) // no slip BC ; plot([u1,u2],cmm="u1,u2 math ",wait=1,coef=0.1/(u1[].linfty+u2[].linfty)); plot(p,cmm="p", wait=1);