//  gnuplot_surface.edp
//
//  Discussion:
//
//    Write data to a file for processing by gnuplot.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem++/gnuplot_surface/gnuplot_surface.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 August 2015
//
//  Author:
//
//    John Burkardt
//
  cout << "\n";
  cout << "GNUPLOT_SURFACE:\n";
  cout << "  FreeFem++ version:\n";
  cout << "  After solving a problem on a 2D region, write the data U(X,Y)\n";
  cout << "  to files that GNUPLOT can use to create a surface plot.\n";
//
//  Define the region.
//
  border domega ( t = 0.0, 2.0 * pi ) { x = cos(t); y = sin(t); }

  mesh Th = buildmesh ( domega ( 100 ) );

  plot ( Th, wait = true, ps = "gnuplot_surface_mesh.ps" );
//
//  Define the finite element space.
//
  fespace Vh ( Th, P1 );
  Vh uh;
  Vh vh;

  func f = 1.0;
  real epsilon = 1.0E-02;

  solve poisson ( uh, vh, solver = LU ) =
    int2d ( Th ) ( dx ( uh ) * dx ( vh ) + dy ( uh ) * dy ( vh ) )
  + int1d ( Th, domega ) ( epsilon * uh * vh )
  - int2d ( Th ) ( f * vh );
//
//  Make a contour plot of UH.
//
  plot ( uh, nbiso = 50, fill = false, value = 1.0, wait = true, 
  ps = "gnuplot_surface_uh.ps" );
//
//  Write a data file for gnuplot.
//
  {
    ofstream data ( "gnuplot_surface_data.txt" );

    for ( int i = 0; i < Th.nt; i++ )
    {
      for ( int j = 0; j <= 3; j++ )
      {
        int j2 = ( j % 3 );
        data  << "  " << Th[i][j2].x 
              << "  " << Th[i][j2].y
              << "  " << uh[][Vh(i,j2)] << "\n";
      }
      data << "\n";
      data << "\n";
    }
    cout << "\n";
    cout << "Wrote gnuplot data file.\n";
  }
//
//  Write a command file for gnuplot.
//
  {
    ofstream commands ( "gnuplot_surface_commands.txt" );

    commands << "# gnuplot_surface_commands.txt\n";
    commands << "#\n";
    commands << "#  Usage:\n";
    commands << "#    gnuplot gnuplot_surface_commands.txt\n";
    commands << "#\n";
    commands << "set term png\n";
    commands << "set output 'gnuplot_surface.png'\n";
    commands << "set xlabel '<--- X --->'\n";
    commands << "set ylabel '<--- Y --->'\n";
    commands << "set zlabel '<--- Uh(X,Y) -->'\n";
    commands << "set title 'Surface plot created from FreeFem++ data.'\n";
    commands << "unset key\n";
    commands << "set timestamp\n";
    commands << "splot 'gnuplot_surface_data.txt' with lines\n";
    commands << "quit\n";

    cout << "Wrote gnuplot command file.\n";
  }
//
//  Terminate.
//
  cout << "\n";
  cout << "gnuplot_surface:\n";
  cout << "  Normal end of execution.\n";

exit ( 0 );