//  migration.edp
//
//  Discussion:
//
//    Migration and proliferation of biological cells.
//
//    The region is a 3x2 rectangle, with periodic conditions
//    at the left and right ends.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/migration/migration.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 June 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 << "migration\n";
cout << "  FreeFem++ version\n";
cout << "  Model the migration and proliferation of biological cells.\n";

real dt = 0.05;
real uf = 1.0;
real rhoc = 100.0;
//
//  Define the mesh.
//
real Lx = 3.0;
real Ly = 2.0;
mesh Th = square ( 60, 40, [x*Lx,y*Ly] );
//
//  Define the finite element spaces using periodicity.
//
fespace Uh ( Th, P2, periodic=[[3,x],[1,x]] );
fespace Vh ( Th, P1, periodic=[[3,x],[1,x]] );
//
//  Declare variables.
//
Uh rho, lrho, rhoold, rhoh, c, ch, cold;
Vh u1, u2, u, n1, n2, v1, v2, rhop1, cp1;
//
//  Set the function rho(x,y).
//
rho = rhoc * exp ( -40.0*(x-Lx/4.0)^2      - 40.0*(y-Ly/2)^2 )
    + rhoc * exp ( -40.0*(x-3.0*Lx/4.0)^2  - 40.0*(y-Ly/2)^2 )
    + rhoc * exp ( -40.0*(x-0.55*Lx)^2     - 40.0*(y-Ly/2)^2 );
//
//  Adapt the mesh according to rho and periodicity.
//
Th = adaptmesh ( Th, rho, periodic=[[3,x],[1,x]] );
//
//  Plot the adapted mesh.
//
plot ( Th, ps = "migration_mesh_initial.ps" );
//
//  Transfer rho to the new mesh.
//
rho = rho;
c = rho / rhoc;
rhoold = rho;
cold = c;
//
//  Plot the initial density.
//
rhop1 = rho;
plot ( Th, rhop1, nbiso = 50, fill = false, value = true,
  ps = "migration_rho_initial.ps" );
//
//  Define the weak form of the migration equations.
//
problem migration ( [rho, c], [rhoh, ch]) =
  int2d ( Th ) ( rho*rhoh/dt )
- int2d ( Th ) ( convect([v1,v2], -dt,rhoold)*rhoh/dt )
+ int2d ( Th ) ( dx(v1)*rho*rhoh+dy(v2)*rho*rhoh )
+ int2d ( Th ) ( 0.01*dx(rho)*dx(rhoh)+0.01*dy(rho)*dy(rhoh) )
- int2d ( Th ) ( 0.01*rho*(rhoc-rhoold)*rhoh )
+ int2d ( Th ) ( c*ch/dt )
- int2d ( Th ) ( cold*ch/dt )
+ int2d ( Th ) ( dx(c)*dx(ch)+dy(c)*dy(ch) )
- int2d ( Th ) ( 10.0*(rho/rhoc-c)*ch );
//
//  Time loop.
//
for ( int it = 0; it < 20; it++ ) 
{
  for ( int substep = 0; substep < 2; substep++ )
  {
    u1 = - dx ( cold );
    u2 = - dy ( cold );
    v1 = 0.5 * u1 * rhoold / rhoc;
    v2 = 0.5 * u2 * rhoold / rhoc;

    migration;

    Th = adaptmesh ( Th, rho, periodic=[[3,x],[1,x]] );

    rho = rho;
    c = c;
    rhoold = rho;
    cold = c;
  }
//
//  Display the current state.
//
  rhop1 = rho;
  cp1 = c;
  plot ( Th, rhop1, nbiso = 50, fill = false, value = true );
}
//
//  Plot the final mesh.
//
  plot ( Th, ps = "migration_mesh_final.ps" );
//
//
//  Plot the final state.
//
  plot ( Th, rhop1, nbiso = 50, fill = false, value = true,
    ps = "migration_rho_final.ps" );
//
//  Terminate.
//
cout << "\n";
cout << "migration:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );