// This a the rotating hill problem with 2 turns.
// = convection equation 
border a(t=0, 2*pi){x = cos(t);y = sin(t);};// the unit circle
int M=40; 
mesh th = buildmesh(a(M));          // triangulates the disk
fespace Vh(th,P1);

//func r2=(x-0.5)^2 + y^2; 
func real hill(real r1, real r2){return exp(-5*r1)*(1-r2);}

real t1=0; 
func xrot= cos(-t1)*x - sin(-t1)*y;
func yrot= sin(-t1)*x + cos(-t1)*y; 
func rr1=(xrot-0.5)^2 + yrot^2;
func rr2=x^2 + y^2;
func sol= hill(rr1,rr2);

Vh v = hill(  (x-0.5)^2 + y^2  ,  x^2+ y^2);                  // initial condition
plot(th,wait=1);
plot(v,wait=1);


//cout << "Adaptation de maillage" << endl;
//th=adaptmesh(th,v);
//v=v;
//plot(th,v,cmm="nouveau v",wait=1);

real WAIT=0; // 0 : pas de pause apres graphique

real N=10;	// N time steps. 
real Npause=1;	// Pause tous les Npause iterations.
real T=pi; 	// (2*pi);
real dt = T/N;    // time step

Vh f1 = -y, f2 = x;   // rotation velocity

int n;
Vh vold, vex; // work  Finite element function 

real t=0;

for ( n=0; n< N ; n++) {
  t += dt;
  vold=v;

  //- methode avec "convect" de FreeFem++
  v=convect([f1,f2],-dt,vold);	

  //- methode "semi-lagrangienne"
  //t1=dt; v=v(xrot,yrot);    
  
  //- adpatation de maillage 
  //cout << "Adaptation de maillage" << endl;
  //th=adaptmesh(th,v);
  //v=v;
  //plot(th,v,cmm="nouveau v",wait=1);

  if (abs((n+1)% Npause) ==0){

    //cout << " n+1="<< n+1 << ";" ;// << endl;
    //cout << " t="<< t << ";" ;// << endl;

    //- Solution exacte:
    t1=t; vex=sol;
    //plot(vex,cmm="Solution convection: t="+t + 
    // ", min=" + vex[].min + ", max=" +  vex[].max,wait=WAIT);

    //- erreur L2:
    real normvL2=sqrt(int2d(th)((v-vex)*(v-vex)));
    cout << "M=" << M << "; t=" << t << "; ||erreur||_L2 = " << normvL2 << endl;

    //- Solution approchee + AFFICHAGES
    plot(v,cmm="convection: t="+t + ", n="+n+"; min=" + v[].min + ", max=" +  v[].max + ", errL2="+normvL2,wait=WAIT);

  };

};

//- Solution approchee + AFFICHAGES
//plot(v,cmm="convection: t="+t + ", n="+n,wait=1);
plot(v,wait=1,value=1);