// file ttemp-IG-2D
//
// The 2D heat equation in a square domain:
//
// \Omega=(-10.,10.) X (-10.,10.), T=5.
//
// Diffusion matrix: K, K_{11}=1+x^2, K_{22}=1+y^2, K_{12}=K_{21}=xy 
//
// Right hand side: v=v00 * (x1c < x) * (x < x2c) * (y1c < y) * (y < y2c)
//
// Initial temperature: u0=u00*exp(-(x^2+y^2))
// Homogeneous Dirichlet, Neumann or Fourier conditions
//
// P1 or P2 spatial finite elements
//
// Implicit Euler + implicit Gear for time discretization


//
// THE DATA:
//

real L=10., T=5.;
real v00=1., x1c=4., x2c=5., y1c=0.5, y2c=1.5;
real u00=1.;
real p=2.; // coefficient for Fourier condition
int iBC=2; // iBC=0 => Dirichlet, iBC=1 => Neumann, iBC=2 => Fourier
int NT=50;

func K11=1.+x*x;
func K22=1.+y*y;
func K12=x*y;
func v0=v00*(x1c < x)*(x < x2c)*(y1c < y)*(y < y2c);
func u0=u00*exp(-(x*x+y*y));

// AUXILIAR VARIABLES:

real L2=2.*L, dt=T/NT, a0=1./dt, a=2.*a0, aE=a, aG=aE, bG=.5*a0;
int ind=0;

//
// THE MESH:
//

mesh Th=square(50,50,[-L+L2*x,-L+L2*y]);
plot(Th, wait=true);

//
// THE FINITE ELEMENT SPACE:
//

// fespace Vh(Th,P1);
fespace Vh(Th,P2);
Vh u=u0, w, f0=v0, uold, uvold, aold;

//
// CASE iBC=0, DIRICHLET:
//

problem thermic0(u,w,solver=Cholesky,init=ind)
    =int2d(Th)(a*u*w + K11*dx(u)*dx(w) + K22*dy(u)*dy(w)
                     + K12*(dx(u)*dy(w) + dy(u)*dx(w)))
    -int2d(Th)(aold*w) 
    +on(1,2,3,4,u=0.);

//
// CASE iBC=1, NEUMANN:
//

problem thermic1(u,w,solver=Cholesky,init=ind)
    =int2d(Th)(a*u*w + K11*dx(u)*dx(w) + K22*dy(u)*dy(w)
                     + K12*(dx(u)*dy(w) + dy(u)*dx(w)))
    -int2d(Th)(aold*w);

//
// CASE iBC=2, FOURIER:
//

problem thermic2(u,w,solver=Cholesky,init=ind)
    =int2d(Th)(a*u*w + K11*dx(u)*dx(w) + K22*dy(u)*dy(w)
                     + K12*(dx(u)*dy(w) + dy(u)*dx(w)))
    +int1d(Th,1,2,3,4)(p*u*w)
    -int2d(Th)(aold*w);

//
// THE COMPUTATIONS	:
//

if (iBC==0){
    uold=u; aold=aE*uold+f0;
	thermic0; // FIRST EULER
	plot(u, value=true, wait=true, fill=true);
    ind=1;
    uold=u; aold=aE*uold+f0;
	thermic0; // SECOND EULER
	plot(u, value=true, wait=true, fill=true);
    ind=0;
    a=1.5*a0;

	for (real t=dt;t<T;t+=dt){ // (real t=0;t<T;t=t+dt){
		uvold=uold; uold=u;
		aold=aG*uold-bG*uvold+f0;
	    thermic0;
	    plot(u, value=true, wait=true, fill=true);
		ind=1;
	}
}
else if (iBC==1){
    uold=u; aold=aE*uold+f0;
	thermic1; // FIRST EULER
	plot(u, value=true, wait=true, fill=true);
    ind=1;
    uold=u; aold=aE*uold+f0;
	thermic1; // SECOND EULER
	plot(u, value=true, wait=true, fill=true);
    ind=0;
    a=1.5*a0;

	for (real t=0;t<T;t+=dt){
		uvold=uold; uold=u;
		aold=aG*uold-bG*uvold+f0;
	    thermic1;
	    plot(u, value=true, wait=true, fill=true);
		ind=1;
	}
}
else if (iBC==2){
    uold=u; aold=aE*uold+f0;
	thermic1; // FIRST EULER
	plot(u, value=true, wait=true, fill=true);
    ind=1;
    uold=u; aold=aE*uold+f0;
	thermic1; // SECOND EULER
	plot(u, value=true, wait=true, fill=true);
    ind=0;
    a=1.5*a0;

	for (real t=0;t<T;t+=dt){
		uvold=uold; uold=u;
		aold=aG*uold-bG*uvold+f0;
	    thermic2;
	    plot(u, value=true, wait=true, fill=true);
		ind=1;
	}
}