FreeFem++ users | Main >> StokesProblemEnriched

Stokes Problem Enriched

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) is

$\displaystyle{ a((u,\tau),(v,t))=(\nabla u,\nabla v)_\Omega+\tau\,t }$

and is the multiplier is

$\displaystyle{ b(q,(u,\tau))=-(q,div(u))_\Omega-\tau \int_\Omega q }$

Also, $\tau\in \mathbb{R}$ plays the role of the average of p, i.e.

$\displaystyle{\tau=\int_\Omega p. }$

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);

return to 2D examples