//  convect.edp
//
//  Discussion:
//
//    Demonstrate a pure convection problem: the quantity U is convected by
//    a rotating velocity field.
//
//    Freefem++ has a "convect" command, of the form
//
//      convect ( ?, ?, ? )
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/convect/convect.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 July 2015
//
//  Author:
//
//    Olivier Pironneau
//
cout << "\n";
cout << "convect:\n";
cout << "  FreeFem++ version:\n";
cout << "  A simple demonstration of the convect() command.\n";

int n = 20;
int M = 4 * n;
real x0 = 0.25;
real y0 = 0.25;
real dt = 2.0 * pi / M;
//
//  Th: a square mesh on [-1,+1]x[-1,+1].
//
mesh Th = square ( n, n, [ 2 * x - 1, 2 * y - 1 ] );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "convect_mesh.ps" );
//
//  Vh: a finite element space on Vh, using piecewise quadratics.
//
fespace Vh ( Th, P2 );
//
//  u, uold: trial functions. 
//
Vh u;
Vh uold;
//
//  v: test function.
//
Vh v;
//
//  A: defines the problem in variational form.
//
problem A ( u, v ) = 
  int2d ( Th ) ( u * v )
- int2d ( Th ) ( convect ( [ y, -x ], -dt, uold ) * v )
+ on ( 1, 2, 3, 4, u = 0 );
//
//  uold: the initial condition as a Gaussian centered at (0.25,0.25).
//
uold = exp ( - 50.0 * ( ( x - x0 ) ^ 2 + ( y - y0 ) ^ 2 ) );
plot ( uold, dim = 32, fill = true, value = true, ps = "convect_uh.ps" );
//
//  Carry out convection M times, using a time step of 2 pi / M.
//
for ( int m = 0; m < M; m++ )
{
  A;
  if ( m < M - 1 )
  
  plot ( u, dim = 3d, fill = true, value = true );
  uold = u;
}
//
//  Terminate.
//
cout << "\n";
cout << "convect:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );