// file heat-CN-2D
//
// The 2D heat equation in a plate:
//
// Diffusion coefficient: k = k0 + k1 * (y < ymax)
// Initial datum: u0 = u00 (constant)
// Right hand side: v = v00 * (x1c < x) * (x < x2c) * (y1c < y) * (y < y2c)
//
// P1 or P2 spatial finite elements
// Crank-Nicholson for time discretization

real u00 = 1., m0 = 1.;
real k0 = 0.2, k1 = 1.8, ymax = 0.5;
real T = 5., dt = 0.05;
real x1c = 4., x2c = 5.5, y1c = 0.1, y2c = 0.4, v00 = -5.;
real errmesh = .1;
int n = 0, nmesh = 1, jiter = 0;

func u0 = u00;
func m = m0;
func k = k0 + k1 * (y < ymax);
func v0 = 2. * v00 * (x1c < x) * (x < x2c) * (y1c < y) * (y < y2c);

real acn=2./dt;

mesh Th=square(60,10,[6*x,y]);
plot(Th, wait = true);

// fespace Vh(Th,P1);
fespace Vh(Th,P2);
Vh u = u0, w, uold, f0, fx, fy;

// for the flat plate:
problem thermic(u,w,init=jiter) = int2d(Th)(acn*u*w + k*(dx(u)*dx(w) + dy(u)*dy(w)))
- int2d(Th)(f0*w - fx*dx(w) - fy*dy(w))
- int2d(Th)(v0 * w)
+ on(1,2,3,4,u=0.);

// ofstream ff("thermic2D.dat");
for(real t=0;t<T;t+=dt){
	n = n+1;
	uold = u;
	f0 = acn*uold;
	fx = k*dx(uold); fy = k*dy(uold);
	thermic;
	jiter = 1;
//    ff << t << " "<< u(3,0.5) << endl;
    // plot(u, wait=true, fill=true);
	plot(u, value=false, wait=true, fill=true);
};

    plot(u, value=true, wait=true, fill=true);