///////////////////////////////////////////////////////////////////
// Mladen Jurak
// Praktikum primijenjene matematike 2, ak. god. 2004/05
// FreeFEM++
//
//   Primjer evolucijske jednadžbe. Animacija rješenja u medit-u.
////////////////////////////////////////////////////////////////////

include "saveMeditAnim.edp"
include "ispisGreskeT.edp"

mesh Th = square(32,200, [x,5*y]);
string ime="heat1";
savemesh(Th,ime+".mesh");

fespace Vh(Th,P1);

real mi = 0.01, t;


Vh u,v,uu;      //  u=u^{n+1}, uu=u^{n}
Vh f,g,exact;   // Funkcija f ovisi o vremenu. Kako t nije globalna varijabla
                // ne možemo ju zadati formulom.
real dt =0.1;   // vremenski korak 

problem Heat(u,v) =
        int2d(Th)( u*v +mi*dt*(dx(u)*dx(v)+dy(u)*dy(v)) )
       -int2d(Th)( dt*f*v + uu*v)
       +on(1,u=g)
       +on(3,u=0);

real T=3; // T = završno vrijeme simulacije
real param=0.8;  // Parametar u definiciji egzaktnog rješenja

real[int]  vrijeme(100);         // vremenski trenuci
real[int]  error(100);           // apsolutna L^2 greška
real[int]  errorRel(100);        // relativna L^2 greška
t=0;
uu=0;
assert(T/dt <= 100);
int nn=int(T/dt);
for(int m=0; m <= nn; ++m)
{
     t=t+dt;
     vrijeme[m]=t;
     f = param*cos(param*t)*(5-y)*(0.5/5. +y*(x-1)^2 * x^2 ) -mi * sin(param*t) *   // definiramo f 
         ( 2*y*(5-y)*(6*x^2 -6*x+1)-2*(x-1)^2 *x^2); // t je sada dobro definiramo
     g = 0.5*sin(param*t);
     exact = sin(param*t)*(5-y)*(0.5/5. + y*(x-1)^2 * x^2);
     Heat;
     uu=u;
     saveMeditAnim(u[],ime,m+1);
     error[m]= sqrt(int2d(Th)(u-exact)^2);
     errorRel[m]=error[m]/sqrt(int2d(Th)(exact)^2);
     cout << "####### t= "<<t<<"#### L^2 greska = "<< error[m] <<" #######"<<endl;
}
bool zaglavlje=true;
string imeDat="L2err.tex";

ispis(vrijeme, error, errorRel, nn, zaglavlje, imeDat);

// Zadatak: 1) Utvrditi eksperimentalno red konvergencije.
//          2) Staviti f=0 i rješavati problem na (0,1)x(0,1).
//          3) Ispisati evoluciju H^1 greške.