include "VTKout.edp"
// Primjer upotrebe programa Paraview za vizualizaciju. Ispisujemo 
// ulazne podatke za Paraview u klasičnom VTK formatu.

//=================================================================
// == Generiranje triangulacije ===================================


////////////////////////////////////////////////////////////////////
// Generalna napomena o labeliranju granice: UVIJEK TREBA LABELU  //
// PRIDRUŽITI SVAKOM DIJELU GRANICE! Ukoliko se to ne učini, onda //
// će FreeFem++ preostale granice labelirati sam. On će to        //
// učiniti dajući suksesivne brojeve po svom algoritmu (pokušaj   //
// ga odrediti eksperimentalno) tako da je moguće da nekom dijelu //
// nelabelirane granice pridruži neku labelu koju smo već pridru- //
// žili nekom drugom dijelu granice. Time se najčešće mijenjaju   //
// rubni uvjeti, što vodi do grubih grešaka.                      //
////////////////////////////////////////////////////////////////////

int DNO   = 1;
int DESNO1 = 21;
int DESNO2 = 22;
int KROV  = 3;
int LIJEVO2  = 42;
int LIJEVO1  = 41;
int KRUG  = 5;

border a1(t=0,1) { x=t;   y=0; label=DNO;}
border a21(t=0,1) { x=1;  y=t; label=DESNO1;}
border a22(t=1,4) { x=1;  y=t; label=DESNO2;}
border a3(t=0,1) { x=1-t; y=4; label=KROV;}
border a42(t=0,3) { x=0;   y=4-t; label=LIJEVO2;}
border a41(t=3,4) { x=0;   y=4-t; label=LIJEVO1;}
border b(t=0,2*pi) { x=0.5+0.1*cos(t);   y=3+0.1*sin(t); label=KRUG;}

mesh Th = buildmesh( a1(10) + a21(10) + a22(30)+ a3(10) + a42(30) + a41(10) + b(-40) );
 plot(Th,wait=1);
//=================================================================

//== Postavka problema ============================================
//  Jedinice u kojima su izrazene varijable
//  ========================================
//   duljina   ... m
//   akceleracija sile teze ... m/s^2
//   gustoca mase           ... kg/m^3
//   tlak                   ... bar
//   propusnost             ... mD (milidarcy)
//   viskoznost             ... cP (centipoise)
//   brzina                 ... m/dan 
//
//   Faktori konverzije 
//    1 mD = 9,86923E-16   m^2
//    1 cP = 1E-3          Pa s

////////////////////////////////////////////////////////////////////////////////
// Provjera kako je granica labelirana
// int nbtriangles=Th.nt;
// for (int i=0;i<nbtriangles;i++)
//   for (int j=0; j <3; j++){
//      if(Th[i][j].label != 0)
//          cout << i << " " << j << " Th[i][j] = "
//               << Th[i][j] << " x = "<< Th[i][j].x << " , y= "<< Th[i][j].y
//               << ", label=" << Th[i][j].label << endl;
//   }
////////////////////////////////////////////////////////////////////////////////


real mi  = 1;    // viskoznost (cP)
real kx  = 80;   // horizontalna propusnost (mD)
real ky  = 80;   // vertikalna propusnost  (mD)
real rho = 1000; // gustoća mase fluida (kg/m^3)
real g   = 9.81; // ubrzanje sile teze (m/s^2)

real Kx = 10;
real Ky = 1;


real qin  = 1;  // ulazna brzina (m/dan);  q * n = -qin

// Darcyjev zakon u danim jedinicama:
// q = - FACT2 * (1/mi) * K *( grad p - FACT1 * rho * g) 
// (K je tenzor propusnosti)
real FACT1  = 1.0E-5;          // faktor za rho . g
real dan    = 60*60*24;
real DARCY  = 9.86923E-8;      // 10^{-8} se skratilo!
real FACT2  = DARCY * dan;
real FACT3 = mi*qin/FACT2;

cout << "FACT2 = " << FACT2 << endl;
cout << "FACT3 = " << FACT3 << endl;

real pout1 = 1;
real pout2 = 4;

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


problem  tlak(u,v, solver=Cholesky) = int2d(Th)( Kx*dx(u)*dx(v) + Ky*dy(u)*dy(v) )
                   + int2d(Th)( FACT1*rho*g*Ky*dy(v))
                   - int1d(Th,KRUG)(FACT3 * v)
                   - int1d(Th,KROV)(0.1*FACT3 * v)
                  // + on(KROV,u=pout1)
                   + on(LIJEVO1,DESNO1,u=pout2);

// Rješavanje
tlak;
// plot(u,wait=1,fill=1,value=1);
//===================================================================


//== Crtanje rješenja ================================================
cout<< "Plotting the solution "<<endl;
  string s="rez2/pressure";
  string soltlak="tlak"; // Bez blankova
//  string bar ="bar";
//  saveP1UDPscal(Th, u[], s, soltlak, bar);
  saveP1VTKscal(Th, u[], s, soltlak);
//===================================================================

// Po dijelovima konstantne funkcije na profinjenoj triangulaciji
fespace Wh(Th,P0);

// Darcyjeva brzina [qx,qy]
Wh qx,qy;
qx = - FACT2 * (1/mi) * Kx * dx(u);
qy = - FACT2 * (1/mi) * Ky *(dy(u)  + FACT1 * rho * g);


string s1="rez2/darcyVel";
string darcyVel="DarcyjevaBrzina";  // Bez blankova
//string unitDar = "m/dan";
//saveP1UDPvect(Th, s1, qx[], qy[],darcyVel,unitDar);
saveP1VTKvect(Th, s1, qx[], qy[],darcyVel);