To avoid the cero mean restriction on pressure we use the variational formulation
Look for \(u\in H_0^1(\Omega)^2\), \(\tau\in \mathbb{R}\) and \(p\in L^2(\Omega)\) such that
$$ \begin{array}{rcll} (\nabla u,\nabla v)_\Omega+\tau\,t-(p,div(u))_\Omega-t\int_\Omega p& =& (f,v) & \mbox{on } \Omega\\ -(q,div(u))_\Omega-\tau \int_\Omega q&=& 0 &\mbox{on } \Omega \end{array} $$
for all \(v\in H_0^1(\Omega)^2\), \(t\in \mathbb{R}\) and \(q\in L^2(\Omega)\).
Here the elliptic part (reason for the name) isAlso, \(\tau\in \mathbb{R}\) plays the role of the average of p, i.e.
download example: Stokes_enriched.edp or return to 2D examples
// // Definition of the mesh // border base(t=0,2){ x=t; y=0; label=1;}; border right(t=0,1){ x=2; y=t; label=1;}; border top(t=0,2){ x=2-t; y=1; label=2;}; border left(t=0,1){ x=0; y=1-t; label=1;}; int m=20; mesh th= buildmesh(base(2*m)+right(m)+top(2*m)+left(m)); func f1=-(x-0.4)^2*(y-1)^3;//sin(x*y+10)*cos(x*y-10); func f2=sin(x*y+10)*cos(x*y-10); func g=0;//boundary data on top lid for u1 fespace Vh(th,P2),Ph(th,P1);// Taylor-Hood Vh ut1,ut2,vt1,vt2;// vitesse Ph pt,qt;// pression macro lap(ut,vt) (dx(ut)*dx(vt)+dy(ut)*dy(vt)) // fin macro // cout<<" WITH MATRICES AND EXPLICITLY DISPLAY ALL UNKNOWS" <<endl; cout<<" USE OF MULTIPLIER FOR MEAN CONDITION "<<endl; // real cpu=clock(); varf at1(ut1,vt1)= int2d(th)(lap(ut1,vt1))+on(1,ut1=0)+on(2,ut1=g);// vitesse first component ut1 varf at2(ut2,vt2)= int2d(th)(lap(ut2,vt2))+on(1,ut2=0)+on(2,ut2=0);// vitesse second component ut2 varf bx(pt,vt1)= int2d(th)(-pt*dx(vt1));// incompressibility varf by(pt,vt1)= int2d(th)(-pt*dy(vt1));// incompressibility varf press(pt,qt)= int2d(th)(-qt);// to compute the integral of each basis function varf ltu1(ut1,vt1)= int2d(th)(f1*vt1)+on(1,ut1=0)+on(2,ut1=g);// right hand side for ut1 varf ltu2(ut1,vt1)= int2d(th)(f2*vt1)+on(1,ut1=0)+on(2,ut1=0);// right hand side for ut2 real tgv = 1e30; matrix At1=at1(Vh,Vh,tgv=tgv);// Matrix for ut1 matrix At2=at2(Vh,Vh,tgv=tgv);// Matrix for ut2 matrix Bx=bx(Ph,Vh); // Matrix for -(p,dx(vt1)) matrix By=by(Ph,Vh);// Matrix for -(p,dx(vt1)) matrix B=[[Bx], [By] ];// Matrix for -(p,div(vt1,vt2)) real[int] rhsut1=ltu1(0,Vh);// vector right hand side for ut1 real[int] rhsut2=ltu2(0,Vh);// vector right hand side for ut2 real[int] pressv=press(0,Ph);// row vector with mean of each basis function matrix E=[ [At1,0], [0,At2] ]; // Matrix for elliptic part of the problem (\Nabla(ut1,ut2),\Nabla(vt1,vt2)) matrix TOTAL=[ [E, B, 0 ], [B', 0, pressv], [0, pressv', 1 ] ];// Matrix for full Stokes problem int ndofu=ut1.n;// size of ut1 and ut2 int ndofp=pt.n;// size of pt int ndof=2*ndofu+ndofp+1; // size of the total unknown (ut1,ut2,pt,tau) real[int] rhs(ndof); real[int] sol(ndof);// vectors to keep solution and rhs rhs(0:ndofu-1)=rhsut1;// rhs part for ut1 rhs(ndofu:2*ndofu-1)=rhsut2;// rhs part for ut2 rhs(2*ndofu:ndof-1)=0; // rhs part for pt and tau cout<<"Talla de rhs "<<rhs.n<<endl; cout<<"Talla de sol "<<sol.n<<endl; set(TOTAL, solver=UMFPACK); sol=TOTAL^-1*rhs;// Solving the system ut1[]=sol(0:ndofu-1); ut2[]=sol(ndofu:2*ndofu-1); pt[]=sol(2*ndofu:ndof-1); real tau=sol(ndof-1); cout<<"Time splitting all variables and multiplier for cero mean= "<<clock()-cpu<<endl; cout<<"Value for tau= "<<tau<<endl; plot(cmm="Stokes solution vitesse: Use of multiplier ",[ut1,ut2],wait=1,fill=1,value=0,dim=2); plot(cmm="Stokes solution pression: Use of multiplier ",pt,wait=1,fill=1,value=0,dim=2);