/*  p'plot.dat'u 1:3,0.538*(1./x/invP)**(1./3.) */
/* OK 13 /03 /2002 trace de la solution semblable */
/* OK janvier 2005   New ++ 01/07 */ 

int nmx=10;
int nmy=10*2;  //multiplier par 4 les deux 

wait=1;
real L=15;
real L1=3;
/* effet d'entree dans un tube             */
border b1(t=0,1)  {x=-L1;       y=1-t;  };
border b2(t=0,L1) {x=-L1+t;     y=0;    };
border b3(t=0,L)  {x=t;         y=0;    };
border b4(t=0,1)  {x=L;         y=t;    };
border b5(t=0,L)  {x=L-t;       y=1;    };
border b6(t=0,L1) {x=-t;        y=1;    };
/* fin de la definition de la geometrie   */
mesh sh = 
 buildmesh(b1(nmy)+ b2(L1*nmx)+b3(L*nmx)+
           b4(nmy)+ b5(L*nmx)+b6(L1*nmx));
plot(sh);
/* construction du maillage               */

fespace Vh2(sh,P2);  

/* Temperature en entree                  */
real T0 =1; 
/* construction du maillage               */
/* vitesse de Poiseuille entre deux plans */
func u0=(1.-y)*y;

/* derivee vitesse de Poiseuille          */
func du02=(1-2*y)*(1-2*y);

/* inverse de Peclet                      */
real invP=1./1000;
/* nombre de Eckert et Peclet unis        */
real Ep=0.1*invP;

Vh2 T,Tp;

 
problem therm(T,Tp) =              
    int2d(sh)(dx(T)*dx(Tp)*invP + dy(T)*dy(Tp)*invP)   
  
  + int2d(sh)(dx(T)*Tp*u0)                          
  
  + int2d(sh)(Tp*Ep)
  
  
  + on(b2,T=0) 
  + on(b3,T=T0)
  + on(b5,T=T0)
  + on(b6,T=0)  
   ;  

therm;

plot(sh,cmm=" Pe=" + 1/invP ,T,fill=1);

//flux=dy(T);
/*plot(flux);*/
//save( 'T.c',T);


{ ofstream ff("T"+ 1/invP+".txt");
   real x,y;
   int i,np=100;
   for (i=0;i<np;i++)
     { x=-L1+ i*(L+L1)/np; 
       y=T(x,0.5);
    ff << x<< " "<< y <<endl;}
}    
// ex p[-0.1:.2]'T1000.txt' u ($1/1000):2,'T100.txt' u ($1/100):2,'T10.txt' u ($1/10):2
 
    
 
/*

p[0:3][0:1]'leveq1.txt','leveq10.txt' ,'leveq100.txt','leveq1000.txt'w l,'' u 3:4,''u 5:6,''u 7:8\
   ,'leveq10000.txt' w l,'' u 3:4 w l,''u 5:6,''u 7:8 w l,1-0.538*x
p[0:3][0:1] 'leveq1000.txt'w l,'' u 3:4,''u 5:6 \
   ,'leveq10000.txt' w l,'' u 3:4 w l,''u 5:6,''u 7:8 w l,1-0.538*x
p[0:3][0:1] 'leveq1000.txt'u ($1:2w l,'' u ($3*10):4,''u 5:6 \
   ,'leveq10000.txt' w l,'' u 3:4 w l,''u 5:6,''u 7:8 w l,1-0.538*x

*/

 { ofstream ff("ture"+ 1/invP+".txt");
   real x,y,eta;
   int i,np=100;
   for (i=0;i<np;i++)
      { x=L/2; 
       y=i*0.5/100 ;
       eta= y/pow(x,(1./3.));
    ff <<  y  << " " << T(L/32,y) << " " 
              << " " << T(L/16,y) << " " 
              << " " << T(L/8 ,y) << " "
              << " " << T(L/4 ,y) << " " 
              << " " << T(L/2 ,y) <<   endl;}
 }             
{ ofstream ff("leveq"+ 1/invP+".txt");
   real x,y,eta;
   int i,np=100;
   for (i=0;i<np;i++)
      { x=L/2; 
       y=i*0.5/10 ;
       eta= y/pow(x,(1./3.));
    ff <<  y/pow(L/16*invP,1./3) << " " << T(L/16,y) << " " 
       <<  y/pow(L/8 *invP,1./3) << " " << T(L/8 ,y) << " "
       <<  y/pow(L/4 *invP,1./3) << " " << T(L/4 ,y) << " " 
       <<  y/pow(L/2 *invP,1./3) << " " << T(L/2 ,y) <<   endl;}
}