// resolution d equation chaleur  PC 1
// PYL  fev 09 dec 17
//
//http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC1/FF++/index.html
//http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC1/FF++/pc1_fourierB1.edp
// 
//
// On veut resoudre par FreeFem++ :                                
//                     d^2  T/dx^2 + d^2 T/dy^2 = 0
//
// avec avec Neumann ou Mixte ou Dirichlet 
//
// run:
// freefem++ pc1_fourierB1.edp
//
//
exec("echo chaleur ");
verbosity=-1;

real s0=clock();
real h0=1;    //hauteur domaine
real L0=1.;   //longueur
int n=25;     //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
mesh Th= buildmesh(b(n)+d(n)+h(n)+g(n));
 
//espace EF    
fespace Vh(Th,P2);
Vh T,TT,Told;
 
real dt = 0.01;// 0.0005;
real t=0;
real errl2=1;  
T=1;
Told=1;
int it=0;

// dirichlet a droite et gauche
problem TureT0T1 (T,TT) =
    int2d(Th)( T*TT/dt +
             ( dx(T)*dx(TT) + dy(T)*dy(TT))  )
   + int2d(Th) ( -(1./dt)*(Told)*TT  
            ) 
 + on(g,T=1)            
 + on(d,T=0) ;
 
 
 // dirichlet a droite et neumann 0 gauche
 problem TureT0P0 (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) 
 ;

// mixte a droite et neumann 0 gauche
real Bi=1;
 problem TureBi1P0 (T,TT) =
    int2d(Th)( T*TT/dt +
             ( dx(T)*dx(TT) + dy(T)*dy(TT))  )
   + int2d(Th) ( -(1./dt)*(Told)*TT  
            ) 
 + int1d(Th,d) (Bi*T*TT) 
 ; 
 
 

 // flux a droite  q=-k dT/dx (attention au signe) et neumann 0 gauche
 real q=.1;
 problem Tureq1P0 (T,TT) =
    int2d(Th)( T*TT/dt +
             ( dx(T)*dx(TT) + dy(T)*dy(TT))  )
   + int2d(Th) ( -(1./dt)*(Told)*TT  
            ) 
 + int1d(Th,d) (q*TT) 
 ; 
 
 
      
 while((t<2))  
 { 
// decommenter pour calculer un cas qui finit lineaire
//    TureT0T1;
// decommenter pour calculer un  qui finit a zero    
//    TureT0P0 ;
// decommenter pour calculer un cas a Bi 
      TureBi1P0;
// decommenter pour calculer un cas a flux constant
//      Tureq1P0;
 
 it++;
 t=t+dt;
 errl2=sqrt(int2d(Th)( (T -Told)^2) );	
 Told=T;
 cout << "t = "<< t <<" ---------------- "<< it <<" +++++++++++++++ errl2 " << errl2<< endl;
 
  if ( !(it % 10))  {
  	  plot(T,cmm="T  t="+t,T,fill=1 ,ps="FILM/T"+it+".eps");

// on sauve dans un fichier qui est ecrase a chaque pas de temps
{ ofstream gnu("NplotT.gp"); 
real x,y,TBi1; 
gnu << "# vals x  T" << endl; 
for (int i=0;i<n;i++) 
   {x =(i*L0)/n; 
    y=h0/2; 
  TBi1=  1.1192*exp(-(0.8603)^2 * t)*cos(0.8603*x)
        -0.1517*exp(-(3.4256)^2 * t)*cos(3.4256*x)  
        +0.0466*exp(-(6.4376)^2 * t)*cos(6.4376*x) 
        -0.0217*exp(-(9.5293)^2 * t)*cos(9.5293*x);
   gnu<< x   << " " << T(x,y) << " " << TBi1 << endl; }
   }
//   
// on peut lancer gnuplot a la volee dans un terminal 
//
exec("echo \" p[0:1][0:1]'NplotT.gp' u 1:2 t'calcul' w l,''u 1:3 t'anal.' w l ;pause .2;quit \"  |  gnuplot"); 
     	}

  } 
 
cout << "CPU " << clock()-s0 << "s " << endl;     

 

// consulter
//http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC1/FF++/index.html