// MAP 562 (G. Allaire, T. Wick)
// Winter 2016/2017
//
// TP 1
// Exercise 1

// We solve the problem:
// - \Delta u = f	in \Omega
//	u=0		on \Gamma_D
//	du/dn=g 	on \Gamma_N

// In general, we have the following structure of a finite
// element program solving a partial differential equation
//
// 0. Define values for right hand side, parameters etc.
// 1. Define the boundary and subsequently create a mesh
// 2. Define the polynomial spaces and finite elements
// 3. Write the specific bilinear form and right hand side functional
//    including realization of the boundary conditions
// 4. Solve the problem
// 5. Visualize the solution


// Right hand side
func f=1;	// volumetric flux

// Traction force (Neumann boundary)
func g=-1;	

// Define a string to provide a name for the solution
string save = "solution";

// Number of edges on the boundary
int N=10;	

// Define a `mesh' variable
mesh Sh;

// Define the boundaries. Here we deal with 
// Dirichlet and Neumann boundary parts
int GammaN,GammaD;
{
	GammaN=1;GammaD=2;
	// Recall the boundary colors: 1 bottom, 2 right, 3 top, 4, left
	// Therefore, we have here on the bottom and left Neumann boundaries
	//  and on the other two parts Dirichlet conditions
	int[int] labels=[GammaN,GammaD,GammaD,GammaN];

	// Our domain will be square
	Sh=square(N,N,label=labels);
}

// We define the polynomials space P1 (linear) for the test functions
macro vh() P1	// finite element space used

// We construct the finite element space in which the polynomials are
// associated with the mesh
fespace Vh(Sh,vh);	

// Next we declare trial and test functions
Vh u,v;			

// The next function is another macro of FreeFem, which much 
// simplifies the usage of derivatives
macro grad(u) [dx(u),dy(u)] //	gradient operator

// We now approach the heart of the problem by realizing the specific 
// problem. At this point 
// we need the weak (i.e., variational formulation) of 
// of our problem: 
//    a(u,v) = (\nabla u, \nabla v) 
//    l(v)   = (f,v) + (g,v)_{\Gamma_N}
//
// Then we bring everything on the left hand side to define a
// `root-finding'-like problem: 
// 
//   a(u,v) - l(v) = 0
problem Laplacian(u,v)=int2d(Sh)(grad(u)'*grad(v))-int2d(Sh)(f*v)-int1d(Sh,GammaN)(g*v)+on(GammaD,u=0);

// By calling the name (here `Laplacian') of our problem, we solve it.
Laplacian;


// Finally we want to get idea of what we actually solved. 
// Here we simply plot the solution
// But we also could have carried out some error analyse to compare our 
// discrete solution with a manufactured solution.
plot(u,cmm="potential",wait=1);

// This second solution, writes into *.esp file
plot(u,cmm="potential",wait=1,ps=save+".eps");