//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/s2de_freefem++/stokes1.edp
//
//  Discussion:
//
//    Steady Stokes problem in a square.
//
//      - uxx - uyy + px = uf
//      - vxx - vyy + py = vf
//        ux  + vy       = pf
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    24 January 2015
//
//  Reference:
//
//    Frederic Hecht,
//    New development in FreeFem++,
//    Journal of Numerical Mathematics,
//    Volume 20, Number 3-4, 2012, pages 251-265.
//
//    Junping Wang, Yanqiu Wang, Xiu Ye,
//    A robust numerical method for Stokes equations based on divergence-free
//    H(div) finite element methods,
//    SIAM Journal on Scientific Computing,
//    Volume 31, Number 4, 2009, pages 2784-2802.
//
real eps = 1.0E-10;
//
//  Define Th, the triangulation of the square using 11 nodes on a side,.
//
int n = 11;

mesh Th = square ( n, n );
//
//  Define Vh, the finite element space for 
//    2D velocity vector (piecewise quadratics) and 
//    pressure (piecewise linear).
//
fespace Vh ( Th, [P2,P2,P1] );
//
//  Define [u1,u2,p], a trial element of the space.
//
Vh [u1,u2,p];
//
//  Define [v1,v2,q], a test element of the space.
//
Vh [v1,v2,q];
//
//  Define the exact functions:
//
func u1exact = - 2.0 
  * x * x * ( x - 1.0 ) * ( x - 1.0 )
  * y * ( y - 1.0 ) * ( 2.0 * y - 1.0 );

func u2exact = + 2.0 
  * x * ( x - 1.0 ) * ( 2.0 * x - 1.0 )
  * y * y * ( y - 1.0 ) * ( y - 1.0 );

func pexact = 0.0;
//
//  Define the right hand side functions.
//
func u1f = 
  2.0 * ( 12.0 * x * x - 12.0 * x + 2 ) 
  * ( 2.0 * y * y * y - 3.0 * y * y + y )
  + 2.0 * ( x * x * x * x - 2.0 * x * x * x + x * x )
  * ( 12.0 * y - 6.0 );

func u2f =
  - 2.0 * ( 12.0 * x - 6.0 )
  * ( y * y * y * y - 2.0 * y * y * y + y * y )
  - 2.0 * ( 2.0 * x * x * x - 3.0 * x * x + x )
  * ( 12.0 * y * y - 12.0 * y + 2 );

func pf = 0.0;
//
//  Define and solve the problem.
//
//  The singularity in the pressure system is treated by
//  adding the term ( eps * q * p ) to the integral.
//  
//
solve Stokes ( [u1,u2,p], [v1,v2,q] ) 
  = int2d ( Th, qforder = 3 ) 
  (
    dx(u1) * dx(v1) 
  + dy(u1) * dy(v1)
  + dx(p)  *    v1
  + dx(u2) * dx(v2) 
  + dy(u2) * dy(v2)
  + dy(p)  *    v2
  + dx(u1) *    q
  + dy(u2) *    q
  - eps * p   * q
  )
  - int2d ( Th, qforder = 3 )
  (
    u1f * v1
  + u2f * v2
  )
  + on ( 1, 2, 3, 4, u1 = u1exact )
  + on ( 1, 2, 3, 4, u2 = u2exact );
//
//  Plot velocity vectors, using only one color for arrows,
//  and multiply default arrow sie by 4.
//
plot ( [u1, u2], nbarrow = 1, coef = 4, ps = "stokes1.ps" );