//  chemotaxis.edp
//
//  Discussion:
//
//    Chemotaxis model for biological pattern generation.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/chemotaxis/chemotaxis.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 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 << "chemotaxis\n";
cout << "  FreeFem++ version\n";
cout << "  Model the motion of cells under the influence of an attractant.\n";
//
// Parameter definition
//
real r = 38.05;
real alpha = 285;
real rhoinf =1;
real D = 0.25;
real s = 1;

real dt = 0.05;
//
//  Define the mesh and plot it.
//
real Lx = 3.5;
real Ly = 4;

mesh Th = square ( 40, 40, [Lx*x, Ly*y] );

plot ( Th, wait = false, ps = "chemotaxis_mesh.ps" );
//
//  Define the finite element spaces.
//
fespace Vh(Th, P1);
Vh rho;
Vh rhoold;
Vh rhotest;
Vh u1;
Vh u2;
Vh output;
Vh sigma;
Vh sigmaold;
Vh sigmatest;

fespace Wh(Th, P2);
Wh c;
Wh cold;
Wh ctest;
//
//  Start from the unstable constant state.
//
rho = rhoinf;
rhoold = rho;
c = rhoinf/(1+rhoinf);
cold=c;
u1 = alpha * dx(c);
u2 = alpha * dy(c);
//
//  Define the weak form of the equation for the density.
//
problem step1(rho, rhotest) =
  int2d(Th)(rho*rhotest /dt)
- int2d(Th)( convect([u1,u2], -dt, rhoold)*rhotest /dt)
+ int2d(Th)(D*dx(rho)*dx(rhotest)+D*dy(rho)*dy(rhotest))
- int2d(Th)(r*s*rhoold*rhoinf*rhotest)
+ int2d(Th)(r*s*rhoold*rho*rhotest);
//
//  Define the weak form of the equation for sigma and c.
//
problem step2([sigma, c], [sigmatest, ctest]) =
  int2d(Th)( sigma*sigmatest /dt)
- int2d(Th)( sigmaold*sigmatest /dt)
- int2d(Th)(alpha*dx(c)*dx(sigmatest)+alpha*dy(c)*dy(sigmatest))
+ int2d(Th)(c*ctest /dt)
- int2d(Th)(cold*ctest /dt)
+ int2d(Th)(dx(c)*dx(ctest) + dy(c)*dy(ctest))
- int2d(Th)(s*rho/(1+rho)*ctest)
+ int2d(Th)(s*c*ctest);

for ( int it = 0; it < 200; it++ ) 
{
  cout << "it = " << it << endl;
//
// Step 1
//
  step1;
  rhoold = rho;
  sigmaold = log(rhoold);
//
//  Step 2
//
  step2;
  u1 = alpha * dx(c);
  u2 = alpha * dy(c);
  rho = exp ( sigma );
  rhoold = rho;
  cold = c;

  if ( ( it % 10 ) == 0 && 0 )
  {
    plot ( c, nbiso = 60, value = true, wait = true );
  }

  if ( ( it % 10 ) == 0 )
  {
    plot ( rho, nbiso = 40, grey = false, fill = true, value = true, 
      wait = false );
  }
}

cout << "cmin = " << c[].min << " cmax = " << c[].max << endl;
cout << "rhomin = " << rho[].min << " rhomax = " << rho[].max << endl;
//
plot( rho, nbiso = 60, grey = true, fill = true, value = false, wait = false,
  ps = "chemotaxis_rho.ps" );

plot ( c, nbiso = 60, grey = true, fill = true, value = false, wait = false,
  ps = "chemotaxis_c.ps");
//
//  Terminate.
//
cout << "\n";
cout << "chemotaxis:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );