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)
where \(c1,\,c_2,\,c_3\) are the lengths of the sides of the triangle \(\kappa\).
download example: stabilized_Stokes.edp or return to 2D examples
// 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));
macro Grad(u) [dx(u),dy(u)]// EOM
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))
//-alpha*square(h)*(Grad(p)'*Grad(q))
) + on(3,u1=g,u2=0)
+ on(1,2,4,u1=0,u2=0);
plot(cmm=" U with Brezzi-Pitkaranta alpha "+alpha,u1,value=1,fill=1,wait=1);
plot(cmm=" V with Brezzi-Pitkaranta alpha "+alpha,u2,value=1,fill=1,wait=1);
plot(cmm=" P with Brezzi-Pitkaranta alpha "+alpha,p,nbiso=50,value=1,fill=1,wait=1);