//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/driven_cavity_stokes/driven_cavity_stokes.edp
//
//  Discussion:
//
//    Steady Stokes problem in a cavity.
//
//      - uxx - uyy + px = uf
//      - vxx - vyy + 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";
//
//  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, u2: velocity trial functions.
//
Vh u1, u2;
//
//  p: pressure trial function.
//
Qh p;
//
//  v1, v2: velocity test functions.
//
Vh v1, v2;
//
//  q: pressure test function.
//
Qh q;
//
//  Define and solve the problem.
//
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
  ) 
  + on ( 1, 2,    4, u1 = 0 )
  + on (       3,    u1 = 1 )
  + on ( 1, 2, 3, 4, u2 = 0 );
//
//  Plot the velocity vectors, using 1 color for the vectors.
//
plot ( [ u1, u2 ], nbarrow = 1, ps = "driven_cavity_stokes_velocity.ps" );
//
//  Plot the pressure.
//
plot ( p, ps = "driven_cavity_stokes_pressure.ps" );
//
//  Terminate.
//
cout << "\n";
cout << "driven_cavity_stokes:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );