//  traffic.edp
//
//  Discussion:
//
//    Pseudo two−dimensional vehicle trafic flow with road branching.
//
//    d/dt rho + del ( rho V(rho) n ) - nu * del^2 rho = 0
//    V(rho) = q(rho) / rho 
//    q(rho) = rho ∗ Vfree ∗ (1−rho/rhoc)
//
//    Here "n" is the unit normal vector.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/traffic/traffic.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    27 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
//
//  Parameters:
//
//    real DT, the time step.
//
//    real NU, the viscosity.
//
//    real T, current time.
//
cout << "\n";
cout << "traffic\n";
cout << "  FreeFem++ version\n";
cout << "  Simulate traffic flow.\n";

real [int] A(2), B(2), C(2), D(2), E(2);
real [int] F(2), G(2), H(2), I(2), J(2);
real [int] K(2), L(2), M(2), Ng(2);
A = [0,0]; B = [9, 1.5]; C = [12, 1.5]; D = [21,0];
E = [6, 2]; F = [15, 2]; G = [0, 4]; H = [9, 2.5];
I = [12, 2.5]; J = [21, 4]; K = [0, 1]; L = [0, 3];
M = [21, 1]; Ng = [21, 3];

border c1(t=0,1) {x=(1-t)*A[0]+t*B[0]; y=(1-t)*A[1]+t*B[1];}
border c2(t=0,1) {x=(1-t)*B[0]+t*C[0]; y=(1-t)*B[1]+t*C[1];}
border c3(t=0,1) {x=(1-t)*C[0]+t*D[0]; y=(1-t)*C[1]+t*D[1];}
border c4(t=0,1) {x=(1-t)*G[0]+t*L[0]; y=(1-t)*G[1]+t*L[1];}
border c5(t=0,1) {x=(1-t)*L[0]+t*E[0]; y=(1-t)*L[1]+t*E[1];}
border c6(t=0,1) {x=(1-t)*E[0]+t*K[0]; y=(1-t)*E[1]+t*K[1];}
border c7(t=0,1) {x=(1-t)*K[0]+t*A[0]; y=(1-t)*K[1]+t*A[1];}
border c8(t=0,1) {x=(1-t)*D[0]+t*M[0]; y=(1-t)*D[1]+t*M[1];}
border c9(t=0,1) {x=(1-t)*M[0]+t*F[0]; y=(1-t)*M[1]+t*F[1];}
border c10(t=0,1) {x=(1-t)*F[0]+t*Ng[0]; y=(1-t)*F[1]+t*Ng[1];}
border c11(t=0,1) {x=(1-t)*Ng[0]+t*J[0]; y=(1-t)*Ng[1]+t*J[1];}
border c12(t=0,1) {x=(1-t)*J[0]+t*I[0]; y=(1-t)*J[1]+t*I[1];}
border c13(t=0,1) {x=(1-t)*I[0]+t*H[0]; y=(1-t)*I[1]+t*H[1];}
border c14(t=0,1) {x=(1-t)*H[0]+t*G[0]; y=(1-t)*H[1]+t*G[1];}
//
//  Define the mesh.
//
mesh Th = buildmesh ( 
    c1(60) + c2(30) + c3(60) + c4(6) + c5(50) + c6(50) + c7(6) + c8(6) 
  + c9(50) + c10(50) + c11(6) + c12(60) + c13(30) + c14(60) );

plot ( Th, ps = "traffic_mesh.ps" );
//
//  Define the finite element spaces.
//
fespace Vh(Th, P1);
Vh rho;
Vh rhoold;
Vh sigma;
Vh sigmaold, rhoh, sigmah, p, ph;
Vh n1, n2, u1, u2, u1visu, u2visu, vrho;
fespace Wh(Th, P2);
Wh phi, phih;
//
real rho0 = 25; // [nb cars / km]
real rhoc = 300; // critical density

real vfree = 110; // [km / h]
real dt = 0.004; 
real nu = 10.0;
real t = 0.0;
//
// Initial field
//
rho = rho0;
//
// Step 0. Define a velocity unit vector by solving an independent Laplace problem. 
//
problem Laplace ( phi, phih ) =
  int2d (Th) ( dx(phi)*dx(phih) + dy(phi)*dy(phih) )
+ int1d (Th, c4) (phih) 
+ int1d (Th, c7) (phih)
+ on ( c8, c11, phi=0 );
//
//  Solve the problem for PHI.
//
Laplace;
//
//  Get the unit vector [N1,N2] of the gradient of the solution PHI.
//
n1 = dx(phi) / sqrt ( dx(phi)^2 + dy(phi)^2 );
n2 = dy(phi) / sqrt ( dx(phi)^2 + dy(phi)^2 );

rhoold = rho;
vrho = vfree * (1.0-rho / rhoc);
u1 = vrho * n1;
u2 = vrho * n2;

problem step1 ( rho, rhoh ) =
  int2d (Th) ( rho*rhoh/dt )
- int2d (Th) ( convect([u1,u2], -dt, rhoold)*rhoh/dt )
+ int2d (Th) ( nu*dx(rho)*dx(rhoh) + nu*dy(rho)*dy(rhoh) )
+ int2d (Th) ( rho*dx(u1)*rhoh + rho*dy(u2)*rhoh )
+ on (c4, rho=rho0 )
+ on (c7, rho=rhoc );

for ( int it = 0 ; it < 200 ; it++ )
{
  t = t + dt;
  step1;

  rhoold = rho;
  vrho = vfree * ( 1.0 - rho / rhoc);
  u1 = vrho * n1;
  u2 = vrho * n2;

  plot ( rho, nbiso = 50, fill = true, value = true, wait = false );
}

cout << "Final time = " << t << endl;

plot ( rho, nbiso = 40, fill = true, value = true, ps = "traffic_rho.ps", 
  cmm = "traffic density" );
u1visu = u1 / vfree;
u2visu = u2 / vfree;
plot ( [ u1visu, u2visu ], value = true, ps = "traffic_velocity.ps",
  cmm = "traffic velocity" );
//
//  Terminate.
//
cout << "\n";
cout << "traffic:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );