/////////////////////////////////////////////////////////////
//                                                         //
//  Solve poisson equation with Omega = square			   //
//                                                         //
//  | - Lap u = f  in Omega                                //
// <                                                       //
//  |       u = g  on dOmega                               //
//                                                         //
/////////////////////////////////////////////////////////////


// Square mesh generation
mesh griglia=square(160,160);
//The command ''square'' triangulates the unit square.
//Generates a 40x40 grid in the unit square [0,1]^2
//mesh Th=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
//Generates a 40x40 grid in the unit square [x0,x1]x[y0,y1]

// FEM space
fespace Vh(griglia,P1); //	The command ''fespace''defined the FEM space
Vh h = hTriangle;
cout << "size of mesh = " << h[].max << endl;
// FE variables
Vh u,v,ues,errore, sorgente;

// Functions
//func f=((x-0.5)^2+(y-0.5)^2<0.25^2) ? 1: 0;  // source (se espressione vera alla f=1 altrimenti 0)
func f=2*16*pi^2*(cos(4*pi*x)*sin(4*pi*y));  // source
func g=sin(4*pi*y);                              // dirichlet
func uesatta=cos(4*pi*x)*sin(4*pi*y);         
// Colors for plot (blue-red in HSV)
real[int] colorhsv=[ 4./6,1,1,  1,1,1 ];
// Colors for plot (green-yellow in HSV)
real[int] colorhsv2=[0.16,0.5,1, 0.5,1,1];


sorgente=f;

plot(sorgente,fill=1,value=0,wait=1,hsv=colorhsv,dim=3); // plot source term 2d
cout << "GRAFICO SORGENTE"<< endl;
//fill= do a color fill between isovalues.
//value=1 to plot the value of isolines 
//wait=1 to stop the plot at each plot
//hsv=(array of float) to defined color of 3*n value in HSV color model




// Weak formulation of the problem
problem poisson(u,v) =  int2d(griglia)( dx(u)*dx(v) + dy(u)*dy(v) ) //  bilinear form
                      - int2d(griglia)( f*v )                       //  linear form
                      + on(1,2,3,4, u=g ) ;                       //  dirichlet boundary condition 

//int2d implement a Gauss quafrature formula  of order 6, exatc on P_5 (int2d(griglia,qforder=2))

poisson;  // solve the problem 
//The default solver is sparsesolver UMFPACK

// Isolevels
real a,b;
a=u[].min;
b=u[].max;
//cout << " min u="<< a<< endl;
//cout << " max u="<< b<< endl;


int nlevels=100;
real dlevels=(b-a)/(nlevels-1);
real[int] viso(nlevels);
for(int i=0; i< nlevels; i++) viso[i]=a+i*dlevels;
//viso= sets the array of isovalues (an array real[int] of increasing values)

plot(u,fill=0,value=0,wait=1,dim=2);// plot solution 2d
plot(u,fill=1,value=1,wait=1, viso=viso(0:viso.n-1),hsv=colorhsv,dim=3); // plot solution 3d
cout << "GRAFICO SOLUZIONE APPROSSIMATA"<< endl;


ues=uesatta;
plot(ues,wait=1);//dim=2 by default
plot(ues,fill=1,value=0,wait=1, viso=viso(0:viso.n-1),hsv=colorhsv2,dim=3); // plot solution 3d    
cout << "GRAFICO SOLUZIONE ESATTA"<< endl;

errore=ues-u;
plot(errore,fill=1,value=0,wait=1, viso=viso(0:viso.n-1),hsv=colorhsv,dim=3); // plot solution 3d
cout << "GRAFICO errore"<< endl;

real eL2=sqrt(int2d(griglia)(errore^2)) ;
real eH1=sqrt(int2d(griglia)( dx(errore)^2 +dy(errore)^2 )+eL2^2);
cout<<"errore in norma L2="<<eL2<<endl;
cout<<"errore in norma H1="<<eH1<<endl;