//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/chrispell_freefem++/scripting.edp
//
//  Discussion:
//
//    We solve the Poisson equation in a square.
//
//      - uxx - uyy = u
//
//    and then:
//
//    * sample the solution, writing the data to a file.
//    * use EXEC() to run GNUPLOT and create a plot of the data.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    27 January 2015
//
//  Reference:
//
//    John Chrispell, Jason Howell,
//    Finite Element Approximation of Partial Differential Equations
//    Using FreeFem++, or, How I Learned to Stop Worrying and Love
//    Numerical Analysis.
//
//    Frederic Hecht,
//    New development in FreeFem++,
//    Journal of Numerical Mathematics,
//    Volume 20, Number 3-4, 2012, pages 251-265.
//

//
//  Define Th, the triangulation of the square.
//
int n = 10;

mesh Th = square ( n, n );
//
//  Define Vh, the finite element space for 
//    2D velocity vector (piecewise linear "bubble") and 
//    pressure (piecewise linear).
//
fespace Vh ( Th, P1 );
//
//  Define UH, the trial function.
//
Vh uh;
//
//  Define VH, test function.
//
Vh vh;
//
//  Define F(X,Y), the right hand side function.
//
func f = 
  +      sin(5*pi*x*(1-x)) *      sin(4*pi*y*(1-y)) * pow(5*pi*(1-x)-5*pi*x,2)
  + 10 * cos(5*pi*x*(1-x)) * pi * sin(4*pi*y*(1-y))
  +      sin(5*pi*x*(1-x)) *      sin(4*pi*y*(1-y)) * pow(4*pi*(1-y)-4*pi*y,2) 
  + 8  * sin(5*pi*x*(1-x)) * pi * cos(4*pi*y*(1-y));
//
//  Define G(X,Y), the boundary condition function.
//
func g = 0.0;
//
//  Define POISSON, the problem to be solved.
//
problem poisson ( uh, vh ) 
  = int2d ( Th ) 
  ( 
    dx(uh) * dx(vh) 
  + dy(uh) * dy(vh)
  )
  - int2d ( Th )
  (
    f * vh 
  )
  + on ( 1, 2, 3, 4, uh = g );
//
//  Solve the POISSON problem.
//  This puts the solution into UH.
//
poisson;
//
//  Define DATAFILE, the name to use for the file storing the data.
//
string datafile = "poisson" + n + ".txt";
//
//  Define PLOTUNIT, an output device associated with "poissonN.sol"
//
  ofstream plotunit ( datafile );
//
//  Write a list of X, Y, UH(X,Y) values to the datafile.
//
{
  for ( int i = 0; i <= n ; i++ )
  {
    real xx = real ( i ) / real ( n );
    for ( int j = 0; j <= n ; j++ )
    {
      real yy = real ( j ) / real ( n );
      plotunit << xx << " " << yy << " " << uh( xx, yy ) << endl;
    }
    plotunit << endl;
  }
}
//
//  Define PLOTFILE, the name of a PostScript file to be created.
//
string plotfile = "poisson" + n + ".eps";
//
//  Using EXEC ( "command" ) create a stream of input to gnuplot.
//
exec
(
  "echo"
  + " '"
  +   "set parametric\n" 
  +   "set term postscript eps enhanced color solid\n"
  +   "set hidden\n"
  +   "set contour base\n"
  +   "set style data lines\n"
  +   "set output \""
  +   plotfile
  +   "\"\n"
  +   "splot \""
  +   datafile
  +   "\"\n"
  + " ' | gnuplot"
);