///////////////////////////////////////////////////////////////////
// Mladen Jurak
// Praktikum primijenjene matematike 2, ak. god. 2004/05
// FreeFEM++
//
// aposteriori1.edp
//
// Primjer adaptacija mreže. 
//
// Poissonova jednažba s homogenim Neumannovim rubnim uvjetom.
// Profinjenje mreže automatski radi rutina adaptmesh().
// Mi računamo indikator greške kako bismo vidjeli gornju među
// za grešku i njenu razdiobu. Kod dobrog profinjenja očekujemo da je 
// greška približno uniformna.
//////////////////////////////////////////////////////////////

//  Domena je u obliku slova L 
border ba(t=0,1)   {x=t;   y=0;   label=1;}
border bb(t=0,0.5) {x=1;   y=t;   label=2;}
border bc(t=0,0.5) {x=1-t; y=0.5; label=3;}
border bd(t=0.5,1) {x=0.5; y=t;   label=4;}
border be(t=0.5,1) {x=1-t; y=1;   label=5;}
border bf(t=0,1)   {x=0;   y=1-t; label=6;}

mesh Th = buildmesh(ba(6) +  bb(4) + bc(4) + bd(4) + be(4) +bf(6));
savemesh(Th, "th.msh");
//plot(Th, wait=true);

// Profinjavanje mreže eksperimentiramo s P1 i P2 elementima.
//fespace Vh(Th,P1);  
fespace Vh(Th,P1);

// Indikator greške ćemo pretvoriti u P0 funkciju radi crtanja.
fespace Nh(Th,P0);
Vh u,v;
Nh rho;

real error=0.01;
func f=(x-y);

// Poissonova jednadžba s homogenim Neumannovim rubnim uvjetom.
// Dodali smo član +1.0e-10 * u da regulariziramo problem (koercitivnost)
// Red integracijske formule zadali smo eksplicitno (qforder=5).
problem Poisson(u,v,solver=CG,eps=1.0e-6) =
      int2d(Th,qforder=5) (dx(u)*dx(v) + dy(u)*dy(v) + u*v*1.0e-10)
     +int2d(Th,qforder=5) (-f*v);

// Indikator greške računamo pomoću varijacijske formulacije. 
// dummy je formalna varijabla koja nema nikakvu ulogu (nužna je radi sintakse)
// (Varijable u varf su posve formalne, što nije slučaj u problem i solve.)
// Ovdje je: 
//          intalledges = integral po svim rubovima svih trokuta. 
//          lenEdge     = globalna varijabla koja daje duljinu stranice
//          hTriangle   = globalna varijabla koja daje dijametar elementa
//          N           = globalna varijabla koja daje vanjsku normalu u danoj točki
//                        na krivulju definiranu s border i sranice trokuta.
//          jump() daje vrijednost skoka (lijevo -desno)
varf indikator(dummy, chiK) =
   intalledges(Th) (chiK*lenEdge*square(jump(N.x*dx(u)+N.y*dy(u))))
   +int2d(Th)( chiK*square(hTriangle*(f+dxx(u)+dyy(u))));

// Radimo 4 profinjenja mreže.
//
for(int i=0; i<4; ++i)
{
    // Riješimo problem
    Poisson;
    cout<< " Minimum rjesenja = "<< u[].min << ", maksimum = "<<u[].max<<endl;
//    plot(u,wait=true);
//  Izračunamo indikator greške. dummy varijabla je stavljena na nulu. 
//  Rezultat je vektor kojiminicijaliziramo P0 funkciju rho.
    rho[]=indikator(0,Nh);
    rho=sqrt(rho);
    cout<<" Minimum indikatora = "<<rho[].min<<", maksimum = "<< rho[].max<<endl;
    //    cmm daje labelu na crtežu 
    //    value = true crta vrijednosti izolinija
    plot(rho, fill=true,wait=true, cmm="Indikator greske", value=true);
//    plot(Th,wait=true); 
//  Adaptiramo mrežu prema gradijenu (P2 elementi)
//  err = P1 interpolacijska greška (0.01 default)
    Th=adaptmesh(Th, [dx(u),dy(u)], err=error);  // stara mreža je uništena
//    Th=adaptmesh(Th, u, err=error);               // stara mreža je uništena
    u=u;                                         // projekcija na novu mrežu 
    rho=rho;                                     // isto
    error=error/2;
}

//  Zadatak. 
//  Ponoviti gornju konstrukciju za zadaću     
//  -div(lambda(x)grad u) + b . grad u + a u = f 
//     + u=0  na granici
//     lambda(x)= 1+x^2
//     b = (1,1)
//     a = 1.0e-6