// example26.edp, from file chapter3/potential.edp

real S=99;
border C(t=0,2*pi) {  x=3*cos(t);  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));

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

solve potential(psi, w) = int2d(Th)(dx(psi)*dx(w) + dy(psi)*dy(w)) +
  on(C, psi=y) + on(S, psi=0);

plot(psi, wait=true);
plot(Th, wait=true);

// solve for heat flow
border D(t=0,2){x=1+t; y=0;}
mesh Sh = buildmesh(C(25) + Splus(-90) + Sminus(-90) + D(200));

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=50;
bool factoredMatrix = false;
problem thermic(v, vv, init=factoredMatrix, 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(int i=0;i<nbT;i++){
    v=vold; 
    thermic;
    factoredMatrix = true;
}

plot(v,wait=1);
plot(Th,wait=1);
plot(Sh,wait=1);