Brezzi-Pitkaranta

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-\alpha\sum_{\kappa}|\kappa|^2(\nabla p,\nabla q)_\kappa &=& 0 &\mbox{on } \Omega \end{array}$$

Here $\kappa$ denotes any of the triangles in the partition of $\bar{\Omega}$ and $|\kappa|^2$ is its area. A value for the parameter $\kappa$ is given by the formula (taken from Gilles Fourestey, These in 2002)

$\displaystyle{ \alpha=\frac{|\kappa|^2}{5(c_1^2+c_2^2+c_3^2)}}$

where $c1,\,c_2,\,c_3$ are the lengths of the sides of the triangle $\kappa$.

// remark: the sign of p is correct
real s0=clock();
border base(t=0,1){ x=t; y=0; label=1;};
border right(t=0,1){ x=1; y=t; label=2;};
border top(t=0,1){ x=1-t; y=1; label=3;};
border left(t=0,1){ x=0; y=1-t; label=4;};
int m=60;
mesh Th= buildmesh(base(m)+right(m)+top(m)+left(m));

fespace Vh(Th,P1);
fespace Ph(Th,P1);
fespace Mh(Th,P0);

Mh h = hTriangle;// For each triangle tau h is the diameter of tau
// atau=h^2/2 approx area of tau
// Estabilization coefficient is
//    atau^2/(5*(c1^2+c2^2+c3^2))
// where ci are the length of each side of tau
//
// Take ci approx h and we find stabilization coefficient to be
//    h^2/60
//plot(h,value=1,fill=1,wait=1);
Vh u1,u2,v1,v2;
Ph p,q;

func g=(x)*(1-x)*4;

real alpha=1.0/60.0;//0.0166666
alpha=alpha;
real  nu=1;
solve Stokes (u1,u2,p,v1,v2,q) =
int2d(Th)( ( dx(u1)*dx(v1) + dy(u1)*dy(v1) + dx(u2)*dx(v2) + dy(u2)*dy(v2) )
-p*q*1e-13- p*dx(v1) - p*dy(v2) - dx(u1)*q - dy(u2)*q
-alpha*square(h)*(dx(p)*dx(q)+dy(p)*dy(q))