//  mesh_points.edp
//
//  Discussion:
//
//    The file "mesh_points_input.txt" lists points that define the boundary.
//
//    Create the corresponding region and mesh it.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mesh_points/mesh_points.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 August 2015
//
cout << "\n";
cout << "mesh_points:\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 = "mesh_points_input.txt";

ifstream file ( filename );
cout << "\n";
cout << "Reading data from '" << filename << "'\n";
//
//  Read the number of points.
//
int nn = 0; 
file >> nn;
cout << "  Number of nodes = "<< nn << endl;
//
//  Read the X and Y coordinates.
//
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;
}
//
//  Create the corresponding border.
//
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 = i + 1; 
}
plot ( Gsup ( nn * 20 ), wait = true, ps = "mesh_points_border.ps", 
  cmm = "mesh_points boundary" );
//
//  Create the mesh.
//
mesh Th = buildmesh ( Gsup ( nn * 20 ) );
plot ( Th, wait = true, ps = "mesh_points_mesh.ps", cmm = "mesh_points mesh" );
//
//  Solve a problem on this mesh.
//
//  Define the finite element space.
//
fespace Vh ( Th, P1 );
//
//  Define test and trial functions.
//
Vh u;
Vh v;
//
//  Define and solve the problem.
//
solve Problem1 ( u, v ) =
  int2d ( Th )( dx(u) * dx(v) + dy(u) * dy(v) ) 
+ int2d ( Th ) ( -v )
+ on ( 1, u = 0 )
+ on ( 2, u = 1 )
+ on ( 3, u = 0 )
+ on ( 4, u = 0 );

plot ( u, wait = true, fill = true, ps = "mesh_points_uh.ps", 
  cmm = "mesh points solution" );
//
//  Terminate.
//
cout << "\n";
cout << "mesh_points:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );