//Matlab jet colormap
real[int] colorJet=[
    0.6667,    1.0000,    0.5625,
    0.6667,    1.0000,    0.6250,
    0.6667,    1.0000,    0.6875,
    0.6667,    1.0000,    0.7500,
    0.6667,    1.0000,    0.8125,
    0.6667,    1.0000,    0.8750,
    0.6667,    1.0000,    0.9375,
    0.6667,    1.0000,    1.0000,
    0.6563,    1.0000,    1.0000,
    0.6458,    1.0000,    1.0000,
    0.6354,    1.0000,    1.0000,
    0.6250,    1.0000,    1.0000,
    0.6146,    1.0000,    1.0000,
    0.6042,    1.0000,    1.0000,
    0.5938,    1.0000,    1.0000,
    0.5833,    1.0000,    1.0000,
    0.5729,    1.0000,    1.0000,
    0.5625,    1.0000,    1.0000,
    0.5521,    1.0000,    1.0000,
    0.5417,    1.0000,    1.0000,
    0.5313,    1.0000,    1.0000,
    0.5208,    1.0000,    1.0000,
    0.5104,    1.0000,    1.0000,
    0.5000,    1.0000,    1.0000,
    0.4889,    0.9375,    1.0000,
    0.4762,    0.8750,    1.0000,
    0.4615,    0.8125,    1.0000,
    0.4444,    0.7500,    1.0000,
    0.4242,    0.6875,    1.0000,
    0.4000,    0.6250,    1.0000,
    0.3704,    0.5625,    1.0000,
    0.3333,    0.5000,    1.0000,
    0.2963,    0.5625,    1.0000,
    0.2667,    0.6250,    1.0000,
    0.2424,    0.6875,    1.0000,
    0.2222,    0.7500,    1.0000,
    0.2051,    0.8125,    1.0000,
    0.1905,    0.8750,    1.0000,
    0.1778,    0.9375,    1.0000,
    0.1667,    1.0000,    1.0000,
    0.1562,    1.0000,    1.0000,
    0.1458,    1.0000,    1.0000,
    0.1354,    1.0000,    1.0000,
    0.1250,    1.0000,    1.0000,
    0.1146,    1.0000,    1.0000,
    0.1042,    1.0000,    1.0000,
    0.0938,    1.0000,    1.0000,
    0.0833,    1.0000,    1.0000,
    0.0729,    1.0000,    1.0000,
    0.0625,    1.0000,    1.0000,
    0.0521,    1.0000,    1.0000,
    0.0417,    1.0000,    1.0000,
    0.0312,    1.0000,    1.0000,
    0.0208,    1.0000,    1.0000,
    0.0104,    1.0000,    1.0000,
         0,    1.0000,    1.0000,
         0,    1.0000,    0.9375,
         0,    1.0000,    0.8750,
         0,    1.0000,    0.8125,
         0,    1.0000,    0.7500,
         0,    1.0000,    0.6875,
         0,    1.0000,    0.6250,
         0,    1.0000,    0.5625,
         0,    1.0000,    0.5000
	];

real kf = 100; //Fluid conductivity
real ks = 1; //Solid conductivity

func Qf=0;
func Qs=0;
func h=1;
func T0=0;
func f=3*x+4*y;

mesh Th=readmesh("Th.mesh");
plot(Th,cmm="Th");

fespace Fh1(Th,P1);
Fh1 T,S;

macro grad(u) [dx(u),dy(u)]//

solve thermic(T,S)=
	int2d(Th,3)(kf*grad(T)'*grad(S))
	+int2d(Th,2)(ks*grad(T)'*grad(S))
	-int2d(Th,3)(Qf*S)
	-int2d(Th,2)(Qs*S)
    -int1d(Th,1)(h*S)
	+on(5,T=T0);	


plot(T,fill=1,cmm="T",value=true,hsv=colorJet);

varf therm(T,S)=int2d(Th,3)(kf*grad(T)'*grad(S))
                    +int2d(Th,2)(ks*grad(T)'*grad(S))
                    +on(5,T=T0);

varf rhs(T,S)=int2d(Th,3)(Qf*S)
             +int2d(Th,2)(Qs*S)
             +int1d(Th,1)(h*S)
             +on(5,T=T0);	

matrix A=therm(Fh1,Fh1,tgv=-1);
real[int] b=rhs(0,Fh1,tgv=-1);

T[] = A^-1*b;
plot(T,fill=1,cmm="T2",value=true,hsv=colorJet);