// M. Jurak, 22.11.2006.
// Ispis funkcija na triangulaciji u leagacy VTK formatu.
// Za Visualization ToolKit (VTK) vidi njegovu web-stranicu  http://www.vtk.org/
// Za VTK file format vidi http://www.vtk.org/documents.php
// 

func int saveP1VTKscal(mesh & Th, real[int] & u, string & ime, string  & var)
{
// ========= Ispis rješenja u VTK formatu =====================
// Rjesenje je skalarna velicina za koju se pretpostavlja da je FE funkcija tipa P1.
//=================================================================
//  Th  = triangulacija;
//  u   = vektor rješenja. Za sada funkcionira samo za P1 elemente (TODO);
//  ime = bazično ime datoteke u koju spremamo rješenje;
//  var = string koji daje naziv varijable; NE SMIJE IMATI BLANKOVE
//  Napomena: Svi stringovi moraju biti preneseni kao reference, jer se inače javljaju greške.


   ofstream  file(ime+".vtk");
  
   // Trenutna verzija je 3.0; Moguće je dati 1.0, 2.0 i 3.0. Verzija 3.0 je
   // kompatibilna sa starijim verzijama.

   file << "# vtk DataFile Version 2.0"<<endl; 
   
   // Naziv podataka (max 255 znakova)
   
   file << "saveP1VTKscal:"+ime+":"+var<<endl;
   file << "ASCII"<<endl;

   // Poliginalni podaci

   file << "DATASET POLYDATA"<<endl;


   int nbt = Th.nt;            // broj trokuta 
   int nbv = Th.nv;            // broj vrhova 
   int nbvdata = nbv;          // broj podataka u vrhovima
   int nbtdata = 0;            // broj podataka u trokutima
   int nbmdata = 0;            // broj modelnih podataka ("model data")

   // Točke, broj točaka (=broj vrhoba trokuta) i dip koordinata
   file << "POINTS " << nbv << " double" << endl;

   fespace ftmp1(Th, P1);  // Trik da se dobiju koordinate vrhova trokuta
   ftmp1 Thx=x, Thy=y;

   //   x-koordinata y-koordinata z-koordinata(0.0)
   for(int i=0; i < nbv; i++) 
      file <<Thx[][i]<<" "<<Thy[][i]<<" "<< "0 \n"; 

   // Trokuti. Za trokute treba uzeti size=br_tokuta*(nr_bvrhova +1). U Deal.II isto
   // uzimaju za četverokute i tetraedre. 
   file << "POLYGONS " << nbt << " " << nbt*4 <<endl;

  // Ispisuje
   // broj vrhova    cell_vert_1     cell_vert_2     cell_vert_3 
   for(int i=0; i < nbt; i++) 
      file <<"3 "<<Th[i][0]<<" "<<Th[i][1]<<" "<<Th[i][2]<<endl;

   // Ovdje počinje skcija s podacima
   file << "POINT_DATA " << nbv << endl;

   // Skalarno polje, ime je u var, tip doouble
   file << "SCALARS " << var << " double"<< endl;

   // Ovo mora biti -- vidi dokumentaciju
   file << "LOOKUP_TABLE default"<<endl;

   // Podaci u vrhovima 
   for(int i=0; i < nbv; i++) 
      file <<u[i]<<endl;

   return 1;
}

/////////////////////////////////////////////////////////////////////////////////////
// Druga verzija funkcije koja ispisuje i vektorsko polje gradijentnog tipa.
// Na primjer, Darcyjevu brzinu.

func int saveP1VTKvect(mesh      & Th, string    & ime, 
                       real[int] & qx, real[int] & qy,
                       string    & var) 
{
 // ========= Ispis rjesenja u leagacy VTK formatu =====================
// Rješenje je skalarna veličina za koju se pretpostavlja da je FE funkcija tipa P1.
// Ispisuje se  "gradijent rješenja", odnosno jedna vektoraska veličina
// tipa P0, tj, konstanta po trokutima. To je vektor [qx,qy].
//=================================================================
//  Th  = triangulacija;
//  u   = vektor rješenja. Za sada funkcionira samp za P1 elemente;
//  ime = bazično ime datoteke u koju spremamo rješenje;
//  var = string koji daje naziv varijable;
//  qx, qy = "gradijent" rješenja - veličina konstantna po trokutima;
//  var = string koji daje naziv varijable [qx,qy]; NE SMIJE IMATI BLANKOVE
//  Napomena: Svi stringovi moraju biti preneseni kao reference, jer se inače javljaju greške.

   ofstream  file(ime+".vtk");
  
   // Trenutna verzija je 3.0; Moguće je dati 1.0, 2.0 i 3.0. Verzija 3.0 je
   // kompatibilna sa starijim verzijama.

   file << "# vtk DataFile Version 2.0"<<endl; 
   
   // Naziv podataka (max 255 znakova)
   
   file << "saveP1VTKscal:"+ime+":"+var<<endl;
   file << "ASCII"<<endl;

   // Poliginalni podaci

   file << "DATASET POLYDATA"<<endl;


   int nbt = Th.nt;            // broj trokuta 
   int nbv = Th.nv;            // broj vrhova 
   int nbvdata = nbv;          // broj podataka u vrhovima
   int nbtdata = 0;            // broj podataka u trokutima
   int nbmdata = 0;            // broj modelnih podataka ("model data")

   // Točke, broj točaka (=broj vrhoba trokuta) i dip koordinata
   file << "POINTS " << nbv << " double" << endl;

   fespace ftmp1(Th, P1);  // Trik da se dobiju koordinate vrhova trokuta
   ftmp1 Thx=x, Thy=y;

   //   x-koordinata y-koordinata z-koordinata(0.0)
   for(int i=0; i < nbv; i++) 
      file <<Thx[][i]<<" "<<Thy[][i]<<" "<< "0 \n"; 

   // Trokuti. Za trokute treba uzeti size=br_tokuta*(nr_bvrhova +1). U Deal.II isto
   // uzimaju za četverokute i tetraedre. 
   file << "POLYGONS " << nbt << " " << nbt*4 <<endl;

  // Ispisuje
   // broj vrhova    cell_vert_1     cell_vert_2     cell_vert_3 
   for(int i=0; i < nbt; i++) 
      file <<"3 "<<Th[i][0]<<" "<<Th[i][1]<<" "<<Th[i][2]<<endl;

   // Ovdje počinje skcija s podacima
   file << "CELL_DATA " << nbt << endl;

   // Skalarno polje, ime je u var, tip doouble
   file << "VECTORS " << var << " double"<< endl;

   // Podaci u elementima 
   for(int i=0; i < nbt; i++) 
      file <<qx[i]<<" "<<qy[i]<<" 0"<< endl;

   
   //file.close();  --> ne postoji

   return 1;
}