Startpage >> Main >> 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\).

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
Vh u1,u2,v1,v2;
Ph p,q;

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

real alpha=1.0/60.0;//0.0166666
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
) + 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);

return to 2D examples

Page last modified on March 28, 2014, at 06:29 PM