///////////////////////////////////////////////////////////////////
// Mladen Jurak
// Praktikum primijenjene matematike 2, ak. god. 2004/05
// FreeFEM++
//
// primjer2.edp
//
// Primjer neregularnog rješenja i adaptacija mreže.
//
// Mješoviti problem za eliptičku jednadžbu može imati rješenje
// čiji je gradijent neograničen čak ako je domena konveksna. 
// U takvom slučaju rješenje dobiveno konačnim elementima neće biti
// dovoljno precizno ako ne koristimo triangulaciju koja je 
// profinjena u okolini singulariteta.
//////////////////////////////////////////////////////////////

// saveMedit = funkcija za ispis rješenja u medit datoteku
include "saveMedit.edp"

// Domena je paravokutnik (-1,1)x(0,1)
border Neumann(t=0,1)  { x=-1+t;  y=0;   label=1; };
border Dirichlet1(t=0,1){ x=t;    y=0;   label=2; };
border Dirichlet2(t=0,1){ x=1;    y=t;   label=2; };
border Dirichlet3(t=0,2){ x=1-t;  y=1;   label=2; };
border Dirichlet4(t=0,1){ x=-1;   y=1-t; label=2; };


mesh Th0 = buildmesh(Neumann(10)     +
                     Dirichlet1(10) +
                     Dirichlet2(10) +
                     Dirichlet3(20) +
                     Dirichlet4(10)
                    );

//plot(Th0, wait=true);

fespace Vh0(Th0,P1);
Vh0 u0, v0;

// Desna strana 
func f = -60*(x^2+y^2);


// Singularni dio točnog rješenja. phi = atan2(y,x) = arctg(y/x)
// (vidi man stranicu za atan2)
func us = sqrt(sqrt(x^2+y^2)) * sin(atan2(y,x)/2);

real K = 10.0;
// Točno rješenje i Dirichletov rubni uvjet 
func ue = K*us + 30*(x^2*y^2);

solve Poisson0(u0,v0) =
      int2d(Th0) ( dx(u0)*dx(v0) +  dy(u0)*dy(v0) )
    - int2d(Th0) (f*v0)
    + on(2, u0=ue);

//plot(u0, Th0, wait=true, value=true, nbiso=100);
string ime = "adapt-0";
saveMedit(Th0, u0[], ime);


//  Adaptiramo mrežu prema singularnom dijelu točnog rješenja:
// Prvo nađimo hmax stare triangulacije
fespace Mh(Th0,P0); Mh htr = hTriangle; real hMax=htr[].max;
// Adaptiramo mrežu ali ne dozvoljavamo da trokuti budu veći od hMax
mesh Th = adaptmesh(Th0,us,hmax=hMax);

// Ponovo riješimo problem, sada  na novoj mreži

fespace Vh(Th,P1);
Vh u, v;

solve Poisson(u,v) =
      int2d(Th) ( dx(u)*dx(v) +  dy(u)*dy(v) )
    - int2d(Th) (f*v)
    + on(2, u=ue);

//plot(u0, Th0, wait=true, value=true, nbiso=100);
string ime1 = "adapt-1";
saveMedit(Th, u[], ime1);

// Izračunajmo relativnu H^1 grešku neadaptiranog i adaptiranog rješenja.
// Prvo izračunamo H1 normu točnog rješenja na finijoj mreži
Vh uue=ue;
real H1e = sqrt( int2d(Th)( dx(uue)^2 +  dy(uue)^2 +uue^2 ) );

// Zatim računamo grešku rješenja na obje triangulacije (u0 i u) i to  
// projicirajući rješenje s neadaptirane triangulacije na adaptiranu.
Vh err0 = u0 - ue;
Vh err  = u  - ue;

real H1err0 = sqrt( int2d(Th)( dx(err0)^2 +  dy(err0)^2 + err0^2 ) );
real H1err  = sqrt( int2d(Th)( dx(err)^2 +  dy(err)^2 + err^2 ) );

cout << "Relativna H1 greska na neadaptiranoj mrezi (br. trokuta = "
         <<Th0.nt<<") = "<<H1err0/H1e<<endl;
cout << "Relativna H1 greska na adaptiranoj mrezi (br. trokuta = "
         <<Th.nt<<") = "<<H1err/H1e<<endl;

//////////////////////////////////////////////////////////////////////
//  Zadatak:  
//          U domeni u obliku slova L konstruirati singularno rješenje 
//          Dirichletove zadaće za Lapalaceovu jednadžbu i usporediti
//          rješenja na uniformnoj i adaptiranoj mreži. Adaptaciju 
//          raditi prema singularnom dijelu točnog rješenja.
///////////////////////////////////////////////////////////////////////