///////////////////////////////////////////////////////////////////
// Mladen Jurak
// Praktikum primijenjene matematike 2, ak. god. 2004/05
// FreeFEM++
//
// Konstrukcija triangulacije 2.
////////////////////////////////////////////////////////////////////


//  NACA0012 airfoil

border upper(t=0,1){
           x=t;
           y=0.17735*sqrt(t)-0.075597*t-0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4);
       }

border lower(t=1,0){
          x=t;
          y=-(0.17735*sqrt(t)-0.075597*t-0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4));
       }

border vanjski(t=0,2*pi){
          x=0.5+0.8*cos(t);
          y=    0.8*sin(t);
       }

//   NAPOMENA:  FreeFem++ ima ugrađene standardne funkcije:
//              sqrt sin cos tan exp log (=ln) log10 (=log) asin acos atan sinh cosh asinh acosh
//              Postoji operator potenciranja (^).

mesh Th = buildmesh(vanjski(30) + upper(35) + lower(35));
plot(Th,wait=true);

// Triangulaciju možemo spremiti u datoteku pomoću metode savemesh

savemesh(Th,"mesh1.msh");

// Trinagulacija se može spremiti u različitim formatima.
// Format je određen s ekstenzijom koju dajemo datoteci.
// Ako je ekstenzija .msh onda će format zapisa biti FreeFem-ov format.
// Triangulacija se može pročitati iz datoteke pomoću readmesh funkcije.

mesh Th1=readmesh("mesh1.msh");
plot(Th1, wait=true);

// Drugi koristan format za spremanje triangulacije (između njih nekoliko) je
// medit-ov format. medit je program za vizualizaciju podataka na netrukturiranim
// mrežama. Datoteka koja sadrži triangulaciju u medif-formatu mora imati
// ekstenziju .mesh

savemesh(Th,"mesh2.mesh");

// Ova naredba generira dvije datoteke: mesh2.mesh i mesh2.mesh.gmsh
// medit pozivamo naredbom
//      medit mesh2
// s komandne linije.
// Pomoću funkcije  exec() možemo izvršiti bilo koju naredbu operacijskog sustava,
// tako da možemo medit pozvati na ovaj način:

//exec("medit mesh2");

//  Funkcija plot pored crtanja triangulacije u ekranu omogućava crtanje u
//  postscript datoteku.

plot(Th, ps="Th.ps");

//  Varijabla tima mesh je jedna "struktura" koja sadrži sve informacije o
//  triangulaciji i koje iz nje možemo pročitati. U sljedećem primjeru pokazujemo
//  sintaksu kojom to možemo postići. Najprije, konstruirajmo dovoljno mali
//  primjer triangulacije:

border C(t=0,2*pi){ x=cos(t); y=sin(t); label=17;}
mesh   Th2=buildmesh(C(10));
savemesh(Th2,"mesh2.mesh");
plot(Th2, wait=true);

//////////////////////////////////////////////////////////////
// TRIANGULACIJA:
//  Do podataka o triangulaciji dolazimo na ovaj način:
//////////////////////////////////////////////////////////////////
assert(version >= 2.19);

int  nt=Th2.nt;                                   // Broj trokuta
cout << "##### Broj trokuta = "<< nt << endl;

int  nv =Th2.nv;                                   // Broj vrhova
cout << "##### Broj vrhova  = "<< nv << endl;

cout << "#### Matrica incidencije #####\n";
cout << "     Za svaki trokut dajemo indekse njegova tri vrha. \n";
// Th2[i]    = trokut (s indeksom i)
// Th2[i][j] = indeks j-tog vrha i-tog trokuta, j=0,1,2.
for(int i = 0; i < nt; ++i)
{
     cout << "trokut "<<i<< ", vrhovi: (" << Th2[i][0]<<","<<Th2[i][1]<<","<<Th2[i][2]<<")\n";
}
cout << "#### Tablica vrhova ######\n";
cout << "     Koordinate svakog vrha i labela.\n";
// Th2(i) je i-ti vrh triangulacije. Dolje je pokazano kako dohvatiti x i y koordinatu
//        točke te njenu labelu. Unutarnje toče imaju labelu 0, a točke na granici dobivaju
//        ili automatski neku pozitivnu labelu, ili onu labelu koju smo koristili pri definiciji granice.
for(int i = 0; i < nv; ++i)
{
     cout << "vrh "<<i<< ", koordinate: (" << Th2(i).x<<","<<Th2(i).y <<"), labela = "<<Th2(i).label<<"\n";
}

cout << "#### Povrsine i labele trokuta ######\n";
// Th2[i].area = površina i-tog trokuta,
//  Th2[i].label = labela i-tog trokuta
for(int i = 0; i < nt; ++i)
{
     cout << "trokut "<<i<< ", povrsina = " << Th2[i].area << ", labela = " << Th2[i].label << "\n";
}

// Za točku (x,y) iz domene možemo naći trokut u kojem se nalazi na sljedeći način: uzmimo, npr. (x,y)=(0,0),
int n = Th2(0.0,0.0).nuTriangle;
cout << " Tocka (0,0) se nalazi u trokutu " << n << endl;


///////////////////////////////////////////////////////////////////////////////////////////
//  REGIONS:
// Pri formiranju triangulacije možemo zadavati unutarnje granice.
// To ima 3 posljedice:
//     1) Triangulacija poštuje unutarnje granice;
//     2) Pomoću unutarnjih granica možemo lokalno profiniti mrežu;
//     3) Domenu možemo podijeliti na subdomene.
// Kao primjer uzmimo jedinični kvadrat s unutarnjim kvadratom

//  (0,1)x(0,1) -- domena
border A5(t=0,1){ x=t; y=0; label = 1;}
border B5(t=0,1){ x=1; y=t; label = 1;}
border C5(t=0,1){ x=1-t; y=1; label = 1;}
border D5(t=0,1){ x=0; y=1-t; label = 1;}

// Unutarnji kvadrat (0.3,0.7)x(0.3,0.7)
border E5(t=0.3,0.7){ x=t; y=0.3; label = 2;}
border F5(t=0.3,0.7){ x=0.7; y=t; label = 2;}
border G5(t=0.3,0.7){ x=1-t; y=0.7; label = 2;}
border H5(t=0.3,0.7){ x=0.3; y=1-t; label = 2;}

mesh Th5 = buildmesh( A5(20) + B5(20) +C5(20) + D5(20) + E5(10) + F5(10) + G5(10) + H5(10) );

plot(Th5, wait=true);

// Kada buildmesh prepozna da je domena podijeljena unutarnjim granicama on će svakoj
// poddomeni pridružiti različit "region" broj. Do njega dolazimo tako da odaberemo
// jednu točku iz poddomene i donjom sintaksom dohvatimo "region" broj.

int reg1 = Th5(0.5,0.5).region; // Uzmi region broj područja u kojem je točka (0.5,0.5)
                                // To će biti unutarnji kvadrat
int reg2 = Th5(0.1,0.1).region; // Uzmi region broj područja u kojem je točka (0.1,0.1)

cout << "region 1 = " << reg1 << ", region 2 = " << reg2 << endl;

// Region broj nam najčešće služi za konstrukciju po dijelovima konstantnih funkcija.
// Globalna varijabla region znade "region" broj u svakoj točki domene.
// Prostor funkcija koje su konstantne po trokutima (na svakom trokutu fj. je konstantna)
fespace Vh(Th5,P0);

real k0 = 1;
real k1 = 2;
// Konstrukcija po dijelovima konstantne funkcije.
// material = k0 u području čiji je "region"  broj  reg1  (unutarnji kvadrat)
// material = k1 u području čiji je "region"  broj  reg2  (vanjski dio)
Vh   material = k0 * (region == reg1) + k1 * (region == reg2);

plot(material, fill=true, value=true, wait=true);

// Zadatak. Napraviti primjer područja s 3 pododmene i step funkciju koja
//          prima različite vrijednosti u različitim poddomenama.

///////////////////////////////////////////////////////////////////////////////////////////
//   FUNKCIJE:  Jedan broj funkcija ugrađen je u FreeFem jezik kao što smo vidjeli.
//              Ostale funkcije dvije varijable možemo zadavati formulama. U tim formulama
//              x i y predstavljaju argumente funkcije. To su globalne varijable koje
//              ne treba deklarirati.
//   Na primjer:

func f = sin(x)*sin(y);
func g = 2*x^2 * y^2;

//   Uočite da se ne piše f(x,y) = ...
//   x i y su realne varijable. Ako želimo kompleksnu fukciju, onda možemo pisati

func z =x + y*1i;       // 1i je oznaka za imaginarnu jedinicu
func h = imag(sqrt(z)); // imag je funkcija koja daje imaginarni dio broja

//   FUNKCIJE:  Konačno, moguće je funkciju definirati kao u C-u, s tom razlikom da se
//              mora koristiti ključna riječ func.

func real myFunc(real t, real s) // ne koristiti x i y!!
{
      if(t < s) return t;
      else      return t*s;
}

//  NAPOMENA: Analogno se definiraju funkcije jedna varijable.


//////////////////////////////////////////////////////////////////////////////////////////
//  TRIANGULACIJA TRANSFORMIRANE DOMENE:  Ako imamo triangulaciju domene U i preslikavanje
//                                        F, onda možemo automatski dobiti triangulaciju
//                                        domene F(U) pomoću funkcije movemesh. Primjer

border G1(t=0,1){x=t; y=0; label=1;}
border G2(t=0,1){x=1; y=t; label=1;}
border G3(t=0,1){x=1-t;y=1;label=1;}
border G4(t=0,1){x=0;y=1-t;label=1;}

mesh Th3 = buildmesh(G1(10)+G2(10)+G3(10)+G4(10));

plot(Th3, wait=1);

real alpha=1;
func Fx = x;
func Fy = y+alpha*y*x*(1-x);

// Prije nego što pozovemo movemesh treba provjeriti hoćemo li dobiti regularnu
// triangulaciju. Funkcija checkmovemesh će dati minimalnu površinu trokuta u
// novoj triangulaciji. Ako je taj broj negativan, onda triangulacija nije korektna.
//
real minT = checkmovemesh(Th3,[Fx,Fy]);
if(minT > 0){
  mesh Th4 = movemesh(Th3,[Fx,Fy]);
  plot(Th4);
  cout << "minT = "<< minT << endl;
}

//  Zadatak. Deformirajte krug tako da dobijete "kružnu domenu" s valovitom
//           granicom. Kontrolirati amplitudu i valnu duljinu valova.
//           Sugestija: polarne koordinate.