//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/driven_cavity_navier_stokes/driven_cavity_navier_stokes.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 MIT license.
//
//  Modified:
//
//    22 May 2020
//
//  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.
//
cout << "\n";
cout << "driven_cavity_stokes:\n";
cout << "  FreeFem++ version.\n";
cout << "  Solve for steady Navier Stokes flow in a cavity.\n";
//
//  eps: the penalty coefficient for the continuity equation.
//
real eps = 1.0E-10;
//
//  n: the number of nodes in each spatial direction.
//
int n = 30;
//
//  Th: the triangulation of the square.
//
mesh Th = square ( n, n );
//
//  Vh: the finite element space for 2D velocity vector (piecewise linear "bubble") 
//
fespace Vh ( Th, P1b );
// 
//  Qh: the finite element space for the scalar pressure (piecewise linear).
//
fespace Qh ( Th, P1 );
//
//  u1, u1o, u2, u2o: velocity trial functions.
//
Vh u1, u1o, u2, u2o;
//
//  p: pressure trial functions.
//
Qh p;
//
//  v1, v2: velocity test functions.
//
Vh v1, v2;
//
//  q: pressure test functions.
//
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.0;
  u2 = 0.0;
//
//  Seek a solution [u1,u2,p] of the nonlinear equations via a fixed point iteration.
//
  for ( int i = 0; i < 10; i++ )
  {
    u1o = u1;
    u2o = u2;
    navierstokes;
  }
//
//  Plot the velocity vectors, using 1 color for the vectors.
//
plot ( [ u1, u2 ], nbarrow = 1, ps = "driven_cavity_navier_stokes_velocity.ps" );
//
//  Plot the pressure.
//
plot ( p, ps = "driven_cavity_navier_stokes_pressure.ps" );
//
//  Terminate.
//
cout << "\n";
cout << "driven_cavity_navier_stokes:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );