// file idealNACA.edp
//
// Computation of the heat generation inside a NACA airfoil
// External potential fluid, parallel to the airfoil
// 
//

real S=99;
real alpha=0.2, uinf=cos(alpha), vinf=sin(alpha);

border C(t=0,2*pi) {  x=3*cos(t)+.5;  y=3*sin(t);}
border Splus(t=0,1){  x = t; y = 0.17735*sqrt(t)-0.075597*t
        - 0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4); label=S;}
border Sminus(t=1,0){  x =t; y= -(0.17735*sqrt(t)-0.075597*t
		-0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4)); label=S;}
mesh Th= buildmesh(C(50)+Splus(70)+Sminus(70));
plot(Th,cmm="FIRST MESH: FLOW COMPUTATION",wait=1);

fespace Vh(Th,P2); Vh psi,w;

solve potential(psi,w)=int2d(Th)(dx(psi)*dx(w)+dy(psi)*dy(w))+
  on(C,psi = uinf*y-vinf*x) + on(S,psi=0.);
plot(psi,cmm="THE STREAMLINES",wait=1);

border D(t=0,2){x=1+t;y=0;}
mesh Sh = buildmesh(C(25)+Splus(-90)+Sminus(-90)+D(200));
plot(Sh,cmm="SECOND MESH: HEAT COMPUTATION",wait=1);

fespace Wh(Sh,P1); Wh v,vv;

int steel=Sh(0.5,0).region, air=Sh(-1,0).region;

fespace W0(Sh,P0);
W0 k=0.01*(region==air)+0.1*(region==steel);
W0 u1=dy(psi)*(region==air), u2=-dx(psi)*(region==air);
Wh vold = 120.*(region==steel);

real dt=0.05, nbT=500;

int i;
problem thermic(v,vv,init=i,solver=LU) = 
		int2d(Sh)(v * vv/dt + k*(dx(v) * dx(vv) + dy(v) * dy(vv))
        + 10. * (u1 * dx(v) + u2 * dy(v)) * vv) - int2d(Sh)(vold * vv/dt);

for (i=0;i<nbT;i++){
	v=vold;
	thermic;
	plot(v,fill=true,cmm="THE TEMPERATURE AND HEAT GENERATION");
}

plot(v,wait=true,fill=true,value=true,cmm="FINAL TEMPERATURE DISTRIBUTION");