//
//
//  A fast  way to solve : 
//  $  \partial_t u - \Delta u = f  $ on $ \Omega \times ]0,T[$
//  $ u = 0 $   for $t=0$ and $ u = g$ on $\partial \Omega$
//  the variationnal form :
//    u^0 = 0;
//   for n = 0 .... 
//   $ \int_\Omega (u^{n+1} - u^n ) v + \nabla u^{n+1}. \nabla v = \int_\Omega fv $
//   ----------------------------------------
func bool GetNoBC(matrix & A,real[int] & in)
{ // def a array in such what
  //  on unkwnon i
  //  in[i] = 1 if no boundary condition 
  //  in[i] = 0 if    boundary condition

  in = A.diag; // take the daig of the matrix
  real tgv = in.max;
  for(int i=0;i<in.n;i++)
    in[i]= in[i] < tgv; // 
  return true; 
}

real dt = 0.01;
mesh Th=square(100,100);
fespace Vh(Th,P1);     // P1 FE space
Vh uh;              // unkown and test function. 
func f=1;                 //  right hand side function 
func g=1;                 //  boundary condition function


varf  vlaplace(uh,vh) =                    //  definion of  the problem 
  int2d(Th)( uh*vh+ dt*(dx(uh)*dx(vh) + dy(uh)*dy(vh)) ) //  bilinear form
  + int2d(Th)( dt*vh*f)
  + on(1,2,3,4,uh=g) ;                      //  boundary condition form

varf vmasse(u,v) = int2d(Th)(u*v);

matrix A = vlaplace(Vh,Vh);		
set(A,solver=UMFPACK);
matrix M = vmasse(Vh,Vh);		
real [int] b(A.n);
real[int] bcl(A.n);
bcl = vlaplace(0,Vh);
real[int] in(A.n); //    1 si interne 0 si frontiere  
GetNoBC(A,in);


Vh p; p[]=in; plot(p,wait=1,cmm="in "); 
for(int i=0;i<10;i++)
  {
    b = M*uh[];
    b = b.* in;
    b += bcl; 
    uh[] = A^-1*b;
    plot(uh,value=true);
    cout << i << endl;
  }