//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/chrispell_freefem++/ns_cavity.edp
//
//  Discussion:
//
//    Steady Navier-Stokes problem in a cavity.
//
//      - uxx - uyy + u ux + v uy + px = uf
//      - vxx - vyy + u vx + v vy + py = vf
//        ux  + vy                     = pf
//
//      u = 1 along top wall.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    26 January 2015
//
//  Reference:
//
//    John Chrispell, Jason Howell,
//    Finite Element Approximation of Partial Differential Equations
//    Using FreeFem++, or, How I Learned to Stop Worrying and Love
//    Numerical Analysis.
//
//    Frederic Hecht,
//    New development in FreeFem++,
//    Journal of Numerical Mathematics,
//    Volume 20, Number 3-4, 2012, pages 251-265.
//
real eps = 1.0E-10;
//
//  Define Th, the triangulation of the square.
//
int n = 30;

mesh Th = square ( n, n );
//
//  Define Vh, the finite element space for 
//    2D velocity vector (piecewise linear "bubble") and 
//    pressure (piecewise linear).
//
fespace Vh ( Th, P1b );
fespace Qh ( Th, P1 );
//
//  Define u1, u1o, u2, u2o, p, trial functions.
//
Vh u1, u1o, u2, u2o;
Qh p;
//
//  Define v1, v2, q, test functions.
//
Vh v1, v2;
Qh q;
//
//  Define the problem.
//
problem navierstokes ( [u1,u2,p], [v1,v2,q] ) 
  = int2d ( Th, qforder = 3 ) 
  ( 
    dx(u1)      * dx(v1) 
  + dy(u1)      * dy(v1)
  + dx(p)       *    v1
  + u1*dx(u1o)  *    v1
  + u2*dy(u1o)  *    v1
  + u1o*dx(u1)  *    v1
  + u2o*dy(u1)  *    v1
  + dx(u2)      * dx(v2) 
  + dy(u2)      * dy(v2)
  + dy(p)       *    v2
  + u1*dx(u2o)  *    v2
  + u2*dy(u2o)  *    v2
  + u1o*dx(u2)  *    v2
  + u2o*dy(u2)  *    v2
  + dx(u1)      *    q
  + dy(u2)      *    q
  - eps * p     *    q
  ) 
  - int2d ( Th, qforder = 3 )
  (
    u1o*dx(u1o) *    v1
  + u2o*dy(u1o) *    v1
  + u1o*dx(u2o) *    v2
  + u2o*dy(u2o) *    v2
  )
  + on ( 1, 2,    4, u1 = 0.0 )
  + on (       3,    u1 = 1.0 )
  + on ( 1, 2, 3, 4, u2 = 0.0 );
//
//  Initialize u1 and u2.
//
  u1 = 0;
  u2 = 0;
//
//  Seek a solution of the nonlinear equations via a fixed point iteration.
//
  for ( int i = 0; i < 10; i++ )
  {
    u1o = u1;
    u2o = u2;
    navierstokes;
  }
//
//  Draw the velocity vectors, using 1 color for the vectors.
//
plot ( [ u1, u2 ], nbarrow = 1, ps = "ns_cavity.ps" );