// resolution d equation chaleur  PC 1 ENSTA
//
// PYL  fev 09 rev 12/10  12/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_fourierBinf.edp
// 
//
// On veut resoudre par FreeFem++ :                                
//                     d^2  T/dx^2 + d^2 T/dy^2 = 0
// avec Neumann en haut bas gauche, dirichelt 1 a droite
//
// run:
// freefem++ pc1_fourierBinf.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;
real t=0;  
T=1;
Told=1; 
 
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)
 ;
 
  
//  
      
 while((t<2))  
 { 
 Ture; 
 t=t+dt;
 Told=T;
 cout << "t = "<< t <<" ----------------   +++++++++++++++" << endl;
 
  	  plot(T,cmm="T  t="+t,T,fill=1,wait=1);
 

{ ofstream gnu("NplotT.gp"); 
real x,y,TBinf; 
gnu << "# vals x  T" << endl; 
for (int i=0;i<=n;i++) 
   {x =(i*L0)/n; 
    y=h0/2; 
  TBinf=    1.27324*cos(1.5708 *x)*exp(-(1.5708)^2* t) 
          -0.424413*cos(4.71239*x)*exp(-(4.71239^2)*t)
          +0.254648*cos(7.85398*x)*exp(-(7.85398^2)*t)
          -0.181891*cos(10.9956*x)*exp(-(10.9956^2)*t);
//  1 x 2 solution 3 quatre modes 5 mode principal          
   gnu<< x   << " " << T(x,y) << " " << TBinf << "  " <<  1.27324*cos(1.5708 *x)*exp(-(1.5708)^2* t)  << endl; }
   }
  } 
 
cout << "CPU " << clock()-s0 << "s " << endl;     

 

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