// example 2
// Stokes equations : regularized cavity flow problem 
// saddle point solver by a preconditioned conjugate gradient method
// COE-tutorail 2007, Atsushi Suzuki 2007/03/08
//
mesh Th = square(20, 20);
fespace Xh(Th, P2), Mh(Th, P1);
Xh u1, u2;
Mh p, q, ppp, qqq, pprec;
real qnorms;
qqq = 1.0;
qnorms = qqq[] '* qqq[];

varf bx(u1,q) = int2d(Th)( dx(u1)*q ); 
varf by(u1,q) = int2d(Th)( dy(u1)*q );
varf a(u1,u2) = int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) )
                    +  on(1,2,4, u1=0)  +  on(3, u1=1) ;
varf l2innerp(p,q) = int2d(Th)(p*q);

matrix A= a(Xh, Xh, solver = CG); 
matrix Bx= bx(Xh, Mh);
matrix By= by(Xh, Mh);
matrix MassP = l2innerp(Mh, Mh, solver = CG);

Xh bc1; bc1[] = a(0,Xh);
Xh b;
Xh bcx, bcy;
bcx=(1.0-x)*x*4.0;
bcy=0;

func real[int] divup(real[int] & pp)
{
   real ptmp = 0.0;
   verbosity=0;
   b[]  = Bx'*pp; b[] += bc1[] .*bcx[];
   u1[] = A^-1*b[];
   b[]  = By'*pp; b[] += bc1[] .*bcy[];
   u2[] = A^-1*b[];
   ppp[] =   Bx*u1[];
   ppp[] +=  By*u2[];
// projection
   ptmp = ppp[] '* qqq[];
   ppp[] -= (ptmp / qnorms) * qqq[];
   return ppp[] ;
};

func real[int] PressurePrecond(real[int] & xx)
{
   real qtmp = 0.0;
   pprec[] = MassP^-1*xx;
// projection
   qtmp = pprec[] '* qqq[];
   pprec[] -= (qtmp / qnorms) * qqq[];
   return pprec[];
}
        
p = 0;
LinearCG(divup,p[],eps=1.e-6,nbiter=50,precon=PressurePrecond);

divup(p[]);
real area = int2d(Th)(1.0);
real meanp = int2d(Th)(p) / area;
cout << "mean pressure = " << meanp << endl;
p = p - meanp;

plot([u1,u2], p, wait=1, value=true, coef=0.1);