//  cloud.edp
//
//  Discussion:
//
//    The file "cloud_points.txt" lists points that define the boundary.
//
//    Create the corresponding region and mesh it.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 August 2015
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Alessandra Agostini, Chiara Piazzola,
//    Mesh Generation with FreeFem++ and Triangle,
//    Master's degree presentation,
//    University of Verona, 2014.
//
cout << "\n";
cout << "cloud:\n";
cout << "  FreeFem++ version.\n";
cout << "  Read a file of point coordinates that outline a region.\n";
cout << "  Create the region and mesh it.\n";
//
//  Attach the file as input.
//
string filename = "cloud_points.txt";

ifstream file ( filename );
cout << "\n";
cout << "Reading data from '" << filename << "'\n";
//
//  nn: the number of points.
//
int nn = 0; 
file >> nn;
cout << "  Number of nodes = "<< nn << endl;
//
//  Gxx, Gyy: the X and Y coordinates of the boundary.
//
real[int] Gxx(nn);
real[int] Gyy(nn);
for ( int i = 0; i < nn; i++ ) 
{
  file >> Gxx[i] >> Gyy[i];
}
//
//  List the nodes.
//
cout << "\n";
cout << "  Boundary nodes:\n";
cout << "\n";
for ( int i = 0; i < nn; i++ ) 
{
  cout << "  " << i
       << "  " << Gxx[i]
       << "  " << Gyy[i] << endl;
}
//
//  Gsup: the border defined by the points.
//  All points get label 1.
//
border Gsup ( t = 0, nn - 1 ) 
{
  int i = min ( int ( t ), Gxx.n - 2 ); 
  real t1 = t - i;
  x = Gxx[i] * ( 1.0 - ( t - i ) ) + Gxx[i+1] * ( t - i ); 
  y = Gyy[i] * ( 1.0 - ( t - i ) ) + Gyy[i+1] * ( t - i ); 
  label = 1; 
}
plot ( Gsup ( nn * 1 ), wait = true, ps = "cloud_border.ps" );
//
//  Th: the mesh.
//
mesh Th = buildmesh ( Gsup ( nn * 1 ) );
//
//  Plot the mesh.

//
plot ( Th, wait = true, ps = "cloud_mesh.ps" );
//
//  Solve a problem on this mesh.
//
//  Vh: the finite element space on mesh Vh, with piecewise linear polynomials.
//
fespace Vh ( Th, P1 );
//
// u: trial function
//
Vh u;
//
//  v: test function.
//
Vh v;
//
//  Define and solve the problem.
//  The boundary condition is associated with the points labeled 1.
//
solve Problem1 ( u, v ) =
  int2d ( Th ) ( dx(u) * dx(v) + dy(u) * dy(v) ) 
+ int2d ( Th ) ( -v )
+ on ( 1, u = 0.0 );
//
//  Plot the solution.
//
plot ( u, wait = true, ps = "cloud_uh.ps" );
//
//  Terminate.
//
cout << "\n";
cout << "cloud:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );