// PC du 30/01/17 
// Equation de la chaleur dans un rectangle
//  PYL + ENSTA à jour 12/17
// http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC1/FF++/pc170130.edp
// http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC1/FF++/index.html
//
// run:
// freefem++ pc170130.edp
//
// geometrie
real h0=1;   //hauteur domaine
real L0=1.;   //longueur
int n=32;    //nbre de points

// definition des cotes Maillage 
border b(t=0,1) { x= t*L0;    y = 0      ; };
border d(t=0,1) { x= L0;      y = h0 * t ; };
border h(t=1,0) { x= L0*(t);  y = h0     ; };
border g(t=1,0) { x= 0;       y = h0 * t ; };

// maillage "Th"
mesh Th= buildmesh(b(n)+d(n)+h(n)+g(n));

// affichage du beau maillage
plot(Th,wait=1);


//definition de l'espace Elements Finis en P2 sur maillage "Th"   
fespace Vh(Th,P2);
// definition des champs de temperature, de test et de temperature passee
Vh T,TT,Told;

// pas de temps et temps 
real dt = 0.0025;
real t=0;  
// initialement la temperature vaut T=1
T=1;
Told=1; 
 
// le probleme a resoudre 
// T temperature, TT champ test
// on reconnait la formulation variationelle
// la temperature au pas précédent passe en terme source
// sur le bord droit "d", on impose T=0
// sur les autres bords dT/dn = 0
// ce qui fait une invariance par translation en haut et en bas, 
// et une condition de symétrie ou d'adiabaticité à gauche (g) 
problem Ture (T,TT) =
    int2d(Th)( T*TT/dt +
             ( dx(T)*dx(TT) + dy(T)*dy(TT))  )
   + int2d(Th) ( -(1./dt)*(Told)*TT  
            )  
 + on(d,T=0);
//
// ouverture du fichier pour sauver toutes les températures
 ofstream gnu("NplotT.gp");
// boucle "temps que"  
  while((t<0.5))  
 { 
// on resout 
 Ture; 
// le temps avance 
 t=t+dt;
// l'ancien devient le nouveau 
 Told=T; 
// on sauve au milieu en h0/2 sur une ligne de 0 a L0
for (int i=0;i<=n;i++) 
{  y = h0/2;  x=(i*L0)/n; 
// un mode et deux modes
  real T1m = 1.2734*cos(1.5708*x)*exp(-((1.5708)^2)*t );
  real T2m = T1m - 0.424413* cos(4.71239*x)*exp(-((4.71239)^2)*t );
// on sauve
gnu<< x <<" "<<  T(x,y)    << " "  <<  t << " " << T2m << endl; }
gnu << endl; 
// on trace
 plot(T,cmm="T  t="+t,T,fill=1,wait=0);
 }
//
//
//
// pour lancer dans un terminal:
// FreeFem++ pc170130.edp ou freefem++ pc170130.edp
// dans une autre fenetre on lance gnuplot
// on trace tous les temps
// p[0:1][0:1]'NplotT.gp' u 1:2 w l
// on peut tracer en 3D 
//sp[0:1][0:1]'NplotT.gp' u 1:3:2 w l
// on trace le temps t=0.15 (dans la colonne 3 on a sauvé le temps)
//p[0:1][0:1]'NplotT.gp'   u 1:($3==0.15? $2:NaN) w l 
// on trace tous les temps en rouge, puis le temps t=0.1 en vert, et deux modes de la  solution exacte sauvés en colonne 4  en bleu
// p[0:1][0:1]'NplotT.gp' u 1:2 w l,'' u 1:($3==0.15? $2:NaN) w l,'' u 1:($3==0.15? $4:NaN) w p 
// on trace la solution autosemblable erf(eta/2), on regarde les temps > 0.01 car le maillage est un peu trop grossier
// et il faut un petit temps d'adaptation.
//on arrete de regarder pour les temps > 0.1 car la temperature a franchement diminue au centre
// p[0:6]'NplotT.gp' u ((1-$1)/sqrt($3)):(($3 <.1)&&($3>0.01) ? $2: NaN) w l,erf(x/2)
//
// plot de la température au centre et sauvegarde dans un fichier "ps"
// gnuplot> p[0:][0:]'NplotT.gp' u 3:($1==0.? $2:NaN) t'T(x=0,t)', 4/pi*exp(-pi*pi*x/4) t '4/pi*exp(-(pi/2)^2 t)'
// gnuplot> set term post
// gnuplot> set outpout 'Tcentre.ps'
// gnuplot> replot
// gnuplot> set term X11
//  ensuite faire "ps2pdf Tcentre.ps"
//
//