Look for \(u\in H^1_0(\Omega)^2\) and \(p\in L^2_0(\Omega)\) such that
$$ \begin{array}{rcll} (\nabla u,\nabla v)_\Omega-(p,div(u))_\Omega& =& (f,v) & \mbox{on } \Omega\\ -(q,div(u))_\Omega&=& 0 &\mbox{on } \Omega \end{array} $$
We use Taylor-Hood and solve Stokes constructing the submatrices for each variable and ensambling to obtain the global matrix.
It seems that this is a faster way to compute the solution than solve and packing the variables in a vector. Moreover, this last procedure gives a scrambled, disordered, vector.
download this example: Stokes_testing.edp or return to 2D examples
//////////////////////////////////////////////////////////////
// We solve Stokes with Taylor-Hood in several simple ways
// and test computing times
////////////////////////////////////////////////////////
//
// 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 f=0;//-(x-0.4)^2*(y-1)^3;//sin(x*y+10)*cos(x*y-10);
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;
//
real cpu=clock();
varf at1(ut1,vt1)= int2d(th)(lap(ut1,vt1))+on(1,ut1=0)+on(2,ut1=1);// 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 mass(pt,qt)= int2d(th)(pt*qt);// penalization
varf ltu1(ut1,vt1)= int2d(th)(f*vt1)+on(1,ut1=0)+on(2,ut1=1);// right hand side for ut1
varf ltu2(ut1,vt1)= int2d(th)(f*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))
matrix M=mass(Ph,Ph); // Penalization matrix: it removes the zero mean restriccion on pressure
M=1e-13*M;
real[int] rhsut1=ltu1(0,Vh);// vector right hand side for ut1
real[int] rhsut2=ltu2(0,Vh);// vector right hand side for ut2
matrix E=[
[At1,0],
[0,At2]
]; // Matrix for elliptic part of the problem (\Nabla(ut1,ut2),\Nabla(vt1,vt2))
matrix TOTAL=[
[E,B],
[B',M]
];// 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; // size of the total unknown (ut1,ut2,pt)
real[int] rhs(ndof),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
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);
cout<<"Time splitting all variables= "<<clock()-cpu<<endl;
plot(cmm="Stokes solution vitesse ",[ut1,ut2],wait=1,fill=1,value=0,dim=2);
plot(cmm="Stokes solution pression ",pt,wait=1,fill=1,value=0,dim=2);
cpu=clock();
////////////////////////////
cout<<" USE OF SOLVE "<<endl;
////////////////////////////
Vh u,v,uu,vv;
Ph p,pp;
solve stokes(u,v,p,uu,vv,pp,solver=UMFPACK) = int2d(th)(dx(u)*dx(uu)+dy(u)*dy(uu) + dx(v)*dx(vv)+ dy(v)*dy(vv)
+ dx(p)*uu + dy(p)*vv + pp*(dx(u)+dy(v)) - 1e-10*p*pp)
+ on(1,u=0,v=0) + on(2,u=1,v=0);
cout<<"Time NOT splitting all variables= "<<clock()-cpu<<endl;
plot([u,v],p,wait=1);
///////////////////////////
cout<<" PACK UNKNOWS INTO A VECTOR "<<endl;
/////////////////////////////
cpu=clock();
fespace VVh(th,[P2,P2,P1]);
VVh [ug,vg,pg],[ush,vsh,psh];
varf Stokes ([ug,vg,pg],[ush,vsh,psh]) =
int2d(th)( dx(ug)*dx(ush) + dy(ug)*dy(ush) + dx(vg)*dx(vsh) + dy(vg)*dy(vsh) -
pg*psh*(1.e-10) - pg*(dx(ush)+dy(vsh)) - (dx(ug)+dy(vg))*psh) +
on(1,ug=0.,vg=0.)+on(2,ug=1,vg=0);
matrix AA=Stokes(VVh,VVh); //GLOBAL STOKES MATRIX
set(AA,solver=UMFPACK);
real[int] bb= Stokes(0,VVh); //GLOBAL RIGHT HAND SIDE
ug[]= AA^-1*bb;// CORRECT WAY TO STORE SOLUTION
cout<<"Time vectoring all variables= "<<clock()-cpu<<endl;
plot(cmm="Stokes solution vitesse ",[ug,vg],wait=1,fill=1,value=0,dim=2);
plot(cmm="Stokes solution pression ",pg,wait=1,fill=1,value=0,dim=2);
/////////////////////////////////////////////////////
// PACK UNKNOWS INTO A VECTOR: HERE WE SEE THAT THE
cout<<"SOLUTION OBTAINED IS SCRAMBLED (DISORDERED) AS A VECTOR "<<endl;
//////////////////////////////////////////////////////////
cpu=clock();
AA=Stokes(VVh,VVh); //GLOBAL STOKES MATRIX
set(AA,solver=UMFPACK);
bb= Stokes(0,VVh); //GLOBAL RIGHT HAND SIDE
sol= AA^-1*bb;// NOT CORRECT WAY TO STORE SOLUTION
ut1[]=sol(0:ndofu-1);
ut2[]=sol(ndofu:2*ndofu-1);
pt[]=sol(2*ndofu:ndof-1);
cout<<"Time vectoring all variables= "<<clock()-cpu<<endl;
plot(cmm="Stokes: vitesse extracted from scrambled vector", [ut1,ut2],wait=1,fill=1,value=0,dim=2);
plot(cmm="Stokes: pression extracted from scrambled vector",pt,wait=1,fill=1,value=0,dim=2);