# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>

using namespace std;

# include "fem_sample.H"

int main ( void );
void test01 ( char *node_file, char *triangle_file, char *coef_file );

//****************************************************************************80

int main ( void )

//****************************************************************************80
//
//  Purpose:
//
//    FEM_SAMPLE_PRB tests routines from the FEM_SAMPLE library.
//
//  Modified:
//
//    12 October 2006
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );

  cout << "\n";
  cout << "FEM_SAMPLE_PRB\n";
  cout << "  C++ version\n";
  cout << "  Test the FEM_SAMPLE library.\n";

  test01 ( "nodes.txt", "triangles.txt", "coefs.txt" );

  cout << "\n";
  cout << "FEM_SAMPLE_PRB\n";
  cout << "  Normal end of execution.\n";

  cout << "\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

void test01 ( char *node_file, char *triangle_file, char *coef_file )

//****************************************************************************80
//
//  Purpose:
//
//    TEST01 demonstrates the use of GRID_SAMPLE.
//
//  Discussion:
//
//    This program assumes that you have computed the value of some scalar
//    quantity (such as pressure or temperature) at a set of nodes, and that
//    these nodes have been triangulated, forming a network of triangles.
//
//    This program can read that data, and evalute the finite element
//    function at a uniformly spaced set of nodes.
//
//    The program reads three data files:
//
//    * node_file, containing X, Y coordinates of points;
//      
//    * triangle_file, containing the nodes that make up
//      a particular triangle;
//      
//    * coef_file, containing the value of a scalar function
//      U(X,Y) evaluated at each node.
//
//  Modified:
//
//    13 October 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string NODE_FILE, the name of the node file.
//
//    Input, string TRIANGLE_FILE, the name of the triangle file.
//
//    Input, string COEF_FILE, the name of the nodal coefficient file.
//
{
  double *coef_node;
  int dim_num;
  int i;
  int j;
  int node_num;
  double *node_xy;
  int number;
  int triangle;
  int *triangle_neighbor;
  int *triangle_node;
  int triangle_num;
  int triangle_order;
  double *u_grid;
  int x_number;
  double *x_vec;
  int y_number;
  double *y_vec;

  cout << "\n";
  cout << "FEM_SAMPLE_TEST01:\n";
  cout << "  C++ version\n";
  cout << "  Evaluate a finite element function on a uniform grid.\n";
  cout << "\n";
  cout << "  This function expects to find three files to read:\n";
  cout << "\n";
  cout << "  * node_file,      containing the node coordinates,\n";
  cout << "  * triangle_file,  containing nodes that form triangles\n";
  cout << "  * coef_file,      containing nodal coefficient values.\n";
  cout << "\n";
  cout << "  It reads the files, constructs a grid of evenly spaced\n";
  cout << "  sample points, and the evaluates the finite element\n";
  cout << "  at the sample points.\n";

  cout << "\n";
  cout << "  The node file is \"" << node_file << "\".\n";
  cout << "  The triangle file is \"" << triangle_file << "\".\n";
  cout << "  The coefficient file is \""<< coef_file << "\".\n";
//
//  Read the data.
//
  dtable_header_read (  node_file, &dim_num, &node_num );

  node_xy = dtable_data_read ( node_file, dim_num, node_num );

  cout << "  The number of nodes is " << node_num << "\n";

  itable_header_read ( triangle_file, &triangle_order, &triangle_num );

  triangle_node = itable_data_read ( triangle_file, triangle_order, 
    triangle_num );

  if ( false )
  {
    cout << "\n";
    cout << "  TRIANGLE NEIGHBOR:\n";
    cout << "\n";

    for ( triangle = 0; triangle < triangle_num; triangle++ )
    {
      cout << "  " << setw(8) << triangle
           << "  " << setw(8) << triangle_node[0+triangle*3]
           << "  " << setw(8) << triangle_node[1+triangle*3]
           << "  " << setw(8) << triangle_node[2+triangle*3] << "\n";
    }
  }

  cout << "  The triangle order is " << triangle_order << "\n";
  cout << "  The number of triangles is " << triangle_num <<"\n";

  coef_node = dtable_data_read ( coef_file, 1, node_num );

  cout << "  The coefficient data has been read.\n";

  triangle_neighbor = new int[3*triangle_num];

  if ( triangle_order == 3 )
  {
    triangulation_order3_neighbor_triangles ( triangle_num, 
      triangle_node, triangle_neighbor );
  }
  else if ( triangle_order == 6 )
  {
    triangulation_order6_neighbor_triangles ( triangle_num, 
      triangle_node, triangle_neighbor );
  }
  else
  {
    cout << "\n";
    cout << "FEM_SAMPLE_TEST01 - Fatal error!\n";
    cout << "  The triangle order must be 3 or 6.\n";
    cout << "  But this data has triangle order = " << triangle_order << "\n";
    return;
  }

  for ( triangle = 0; triangle < triangle_num; triangle++ )
  {
    cout << "  " << setw(8) << triangle
         << "  " << setw(8) << triangle_neighbor[0+triangle*3]
         << "  " << setw(8) << triangle_neighbor[1+triangle*3]
         << "  " << setw(8) << triangle_neighbor[2+triangle*3] << "\n";
  }
  cout << "  The triangle neighbor array has been computed.\n";

  number = 4;

  grid_sample_size ( node_num, node_xy, number, &x_number, &y_number );

  u_grid = new double[x_number*y_number];
  x_vec = new double[x_number];
  y_vec = new double[y_number];

  cout << "\n";
  cout << "  Requested   node density = " <<   number << "\n";
  cout << "  Suggested X node density = " << x_number << "\n";
  cout << "  Suggested Y node density = " << y_number << "\n";

  grid_sample ( node_num, node_xy, triangle_order, triangle_num, 
    triangle_node, triangle_neighbor, coef_node, x_number, y_number, 
    x_vec, y_vec, u_grid );

  cout << "\n";
  cout << "  X sample coordinates:\n";
  cout << "\n";

  for ( i = 0; i < x_number; i++ )
  {
    cout << "  " << setw(10) << x_vec[i];
  }
  cout << "\n";
  
  cout << "\n";
  cout << "  Y sample coordinates:\n";
  cout << "\n";

  for ( j = 0; j < y_number; j++ )
  {
    cout << "  " << setw(10) << y_vec[j];
  }
  cout << "\n";
  
  cout << "\n";
  cout << "  U evaluated at sample points:\n";
  cout << "\n";

  for ( j = y_number - 1; 0 <= j; j-- )
  {
    for ( i = 0; i < x_number; i++ )
    {
      cout << "  " << setw(10) << u_grid[i+j*x_number];
    }
    cout << "\n";
  }

  delete [] coef_node;
  delete [] node_xy;
  delete [] triangle_neighbor;
  delete [] triangle_node;
  delete [] u_grid;
  delete [] x_vec;
  delete [] y_vec;

  return;
}