//  channel_stokes_steady.edp
//
//  Discussion:
//
//    Solve the steady Stokes equations.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/channel_stokes_steady/channel_stokes_steady.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    18 May 2015
//
//  Author:
//
//    Florian De Vuyst
//
//  Reference:
//
//    Florian De Vuyst,
//    Numerical modeling of transport problems using freefem++ software -
//    with examples in biology, CFD, traffic flow and energy transfer,
//    HAL id: cel-00842234
//    https://cel.archives-ouvertes.fr/cel-00842234
//
cout << "\n";
cout << "channel_stokes_steady:\n";
cout << "  FreeFem++ version.\n";
cout << "  Solve Stokes equations for steady flow in a channel.\n";

real Re = 600.0;
real nu = 1.0 / Re;
real Lx = 12.0;
real Ly = 5.0;
real dt = 0.5;
//
//  Set the boundary.
//
border c1 ( t=0,1 )    { x=t*Lx; y=0; }
border c2 ( t=0,1 )    { x=Lx; y=t*Ly; }
border c3 ( t=1,0 )    { x=t*Lx; y=Ly; }
border c4 ( t=1,0 )    { x=0; y=t*Ly; }
border c5 ( t=2*pi,0 ) { x=4+0.2*cos(t); y=Ly/2+0.2*sin(t); }
border c6 ( t=0,1 )    { x=5+Lx/2*t; y=Ly/2+0.4; }
border c7 ( t=0,1 )    { x=5+Lx/2*t; y=Ly/2-0.4; }
//
//  Make the mesh.
//
mesh Th = buildmesh 
(
  c1(80)+c2(40)+c3(80)+c4(20)+c5(60)+c6(100)+c7(100)
);

plot ( Th, wait = true, ps = "channel_stokes_steady_mesh.ps" );
//
//  Define discrete velocity space Uh and discrete pressure space Vh.
//
fespace Uh ( Th, P2 );
fespace Vh ( Th, P1 );
Uh u;
Uh v;
Uh uh;
Uh vh;
Uh uold;
Uh vold;
Vh p;
Vh ph;
Vh uplot;
Vh vplot;
Vh vort;
//
//  Define the weak form of the steady state Stokes equations.
//
real eps = 1.0E-10;

problem stokes ( [u,v,p], [uh,vh,ph] ) =
    int2d(Th) ( nu*dx(u)*dx(uh) + nu*dy(u)*dy(uh) )
  + int2d(Th) (nu*dx(v)*dx(vh) + nu*dy(v)*dy(vh) )
  - int2d(Th) (p*dx(uh))
  - int2d(Th) (p*dy(vh))
  - int1d(Th, c2) (nu*dx(u)*N.x*uh+nu*dy(u)*N.y*uh)
  - int1d(Th, c2) (nu*dx(v)*N.x*vh+nu*dy(v)*N.y*vh)
  + int2d(Th) (dx(u)*ph + dy(v)*ph)
  + int2d(Th) (eps*p*ph)
  + on ( c1, c3, c5, u=0, v=0)
  + on ( c4, u=4.0 * y/Ly * (1-y/Ly), v=0 );
//
//  Solve the equations.
//
stokes;
//
//  Plot the velocity vector field.
//
uplot = u;
vplot = v;
plot ( Th, [uplot,vplot], wait = true, fill = true, 
  ps = "channel_stokes_steady_velocity.ps", cmm = "channel_stokes_steady velocity" );
//
//  Plot the pressure field.
//
plot ( Th, p, wait = true, fill = true, ps = "channel_stokes_steady_pressure.ps", 
  cmm = "channel_stokes_steady pressure" );
//
//  Terminate.
//
cout << "\n";
cout << "channel_stokes_steady:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );