/////////////////////////////////////////////////////////////
//                                                         //
//  Solve poisson equation with Omega=cirlce			   //
//                                                         //
//  | - Lap u = f  in Omega                                //
// <    du/dn = gN  on Gamma_N                             //
//  |       u = gD  on Gamma_D							   //
//														   //
//									                       //
//                                                         //
/////////////////////////////////////////////////////////////


// Parametrization of the domain
border dOmega1(t=0,pi){x=cos(t);y=sin(t);}
border dOmega2(t=pi,2*pi){x=cos(t);y=sin(t);}

int nodes=100;

mesh griglia=buildmesh(dOmega1(nodes)+dOmega2(nodes));

// buildmesh calls a triangulation subroutine based on the Delaunay test, which first triangulates with only the boundary points, 
// then adds internal points by subdividing the edges.

plot(dOmega1(nodes)+dOmega2(nodes),wait=1);

plot(griglia,wait=1);


// 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 gD=sin(x);                              // dirichlet
func gN=exp(x*y);     
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,wait=1);

// 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)
						+ on(dOmega1,u=gD)    //  dirichlet boundary condition 
					    -int1d(griglia,dOmega2)(v*gN); //  Neumann boundary condition 



poisson;  // solve the problem 

// 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=1,value=1,wait=1, viso=viso(0:viso.n-1),hsv=colorhsv,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;