//////////////////////////////////////////////////////////
//                                                      //
//  Solve heat equation with Omega = circle of radius 1 //
//                                                      //
//  | u_t = u_xx + u_yy    in Omega x (0,T)             //
// <  u(0) = u0            in Omega                     //
//  | u = g                in dOmega x (0,T)            //
//                                                      //
//////////////////////////////////////////////////////////

verbosity=0; // disable freefem output

// Parametrization of the circle
border dOmega(t=0,2*pi){x=cos(t); y=sin(t);}

// Mesh generation
int nodes = 100;
mesh Omega = buildmesh(dOmega(nodes));

// Finite element space
fespace X(Omega,P1);

// Finite element variables
X u,uold,v;

// Functions
func u0=x^2+y^2;  // initial datum
func g=1;         // dirichlet

// Parameters
real dt=0.01;     // time step 
real T=1;         // final time
int  steps=T/dt;  // number of steps in time

// Isolevels for graphics
int nlevels=40;
real dlevels=1./(nlevels-1);
real[int] levels(nlevels);
for(int i=0; i< nlevels; i++) levels[i]=i*dlevels;

// Initialization (projection of u0 on X) and plot
uold=u0;
plot(uold,fill=true,value=false,wait=true, viso=levels); // plot initial datum

// Backward euler weak formulation of the problem
problem heat(u,v) =  int2d(Omega)( u*v + dt*(dx(u)*dx(v) + dy(u)*dy(v)) ) //  bilinear form
                   - int2d(Omega)( uold*v )                               //  linear form
                   + on( dOmega, u=g ) ;                                  //  dirichlet boundary condition 

// Solve iteratively
for(int i=0; i<steps; i++)
{
 heat; // solve the problem at current time 
 plot(u,fill=true,value=false, viso=levels); // plot solution
 uold=u; // swap u(n) <---> u(n+1) 
 cout<<endl<<"Iteration "<<i+1<<"  -  Time = "<<(i+1)*dt<<endl;
}