//
// resolution de l'équation de la chaleur dans un barreau fin
// PYL   dec 17
//http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/C2cond.ENSTA.pdf page 2
//http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/HTMLPC/PC1/FF++/ailette.edp
//
// On veut resoudre par FreeFem++ :                                
//                     d^2  T/dx^2 + d^2 T/dy^2 = 0
// dans un rectangle fin ,    T impose=1 sur la face de gauche (notee 4),  
//   sur les autres faces (notees 1,2,3) on met un coefficient d'échange, T=0 à l'extérieur
//    -k dT/dn = h T     
// geometrie  du rectangle longueur L0, epaisseur a, a/L0<<1
//
real a=.1;
real L0=1;
// definition des cotes Maillage
border b(t=0,1) { x= t;          y = 0  ;        };
border d(t=0,1) { x= L0;         y =a * t;      } ;
border h(t=0,1) { x= L0*(1 - t); y = a  ;  };
border g(t=0,1) { x= 0;          y =a * (1-t) ;   };
int n=5;
// Oh le beau maillage, faire CR pour continuer (wait=1)
mesh Th = buildmesh(b(15*n)+d(n)+h(15*n)+g(n)); 
plot(Th,wait=1);
// Espace et champs 
fespace Vh2(Th,P2); 
Vh2 T,Tt; 
//Temperature en entree
real T0 =1; 
//Nombre de Biot issu de la théorie: Bi = a h/k 
// c'est lui que l'on se donne car on veut jouer avec lui
// normalement c'est h qui est donne.                                      
real Bi=100;
// H coefficient d'echange reduit H=hL0/k:    donc H = h L0/k = Bi L0/a             
// la longueur L/L0 est     (a/L0)*Bi^(-1/2) 
real H= Bi*L0/a;
// resolution : ecriture sous forme variationelle du problème
// remarquer la condition mixte   sous la forme: int1d(Th, ) (H*T*Tt)
problem ailette(T,Tt)=
   int2d(Th)(dx(T)*dx(Tt)+ dy(T)*dy(Tt)) 
   + int1d(Th,1,2,3) (H*T*Tt)
   + on(4,T=T0)  ;
 
// boucle en Biot, divisé par 10  
for(int i=0;i<7;i++)
{ cout << "coef d echange H=" << H << "   Bi=" << Bi 
              << " L/L0:=" << (a/L0)*Bi^(-1/2)  << endl ;
// résolution              
   ailette;
// affichage de l'ailette et de sa température, faire CR pour continuer (wait=1)   
   plot(T,fill=1,cmm="H="+H + ", min=" + T[].min + ", max=" + T[].max,wait=1);
// sauvegarde en postscript de l'image   
   plot(T,fill=1,cmm="H="+H + ", min=" + T[].min + ", max=" + T[].max,wait=0,ps="TBi"+Bi+".eps");
// sauvegarde gnuplot  de la température au centre (x,y=a/2) et de la solution analytique (uniquement exp, modifier pour mettre le cosh)
  { ofstream gnu("plot.gp"); 
    int nx=50;
    for (int i=0;i<=nx;i++) 
      {  real x=i*1./nx;
         gnu << x << " " << T(x,a/2) << " "  << exp(-x/a*sqrt(2*Bi)) << endl; 
  }   } 
// update du nouveau Bi et du nouvel H
  Bi=Bi/10;
  H= Bi*L0/a; 
}


//
// run:
// freefem++ ailette.pde 
// 
// entre chaque division par 10 de Bi, le programme attend (wait=1) un retour charriot (CR)
// on peut en profiter pour plotter avec gnuplot
// gnuplot>  p[][0:1]'plot.gp' t 'Calcul' ,''u 1:3 t'Anal.'
//
//