//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/s2de_freefem++/stokes2.edp
//
//  Discussion:
//
//    Stokes problem in a square.
//
//  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.
//

//
//  Parameters:
//
real eps = 1.0E-10;
real nu = 0.1;
real rho = 1.0;
real t = 0.0;
//
//  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 Grad(u), a macro for the gradient of a scalar field.
//
macro Grad(u) [dx(u),dy(u)] // EOM
//
//  Define div(u1,u2), a macro for the divergence of a vector field.
//
macro div(u1,u2) (dx(u1)+dy(u2))  // EOM
//
//  Define the exact functions:
//
func u1exact =   2.0 * sin ( 2.0 * pi * x ) * cos ( 2.0 * pi * y );

func u2exact = - 2.0 * cos ( 2.0 * pi * x ) * sin ( 2.0 * pi * y );

func pexact = x * x + y * y;
//
//  Define the right hand side functions:
//
func u1f = 2.0 * x
  - 16.0 * pi * pi * sin ( 2.0 * pi * x ) * cos ( 2.0 * pi * y );

func u2f = 2.0 * y
  - 16.0 * pi * pi * cos ( 2.0 * pi * x ) * sin ( 2.0 * pi * y );

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 )
  ( 
      Grad(u1)'      * Grad(v1) 
    + Grad(u2)'      * Grad(v2) 
    - div ( u1, u2 ) * q 
    - p              * div ( v1, v2 )
    + 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 ( [u1, u2], nbarrow = 1, coef = 0.125, ps = "stokes2.ps" );