program main

!*****************************************************************************80
!
!! MAIN is the main program for FEM_SAMPLE_PRB.
!
!  Discussion:
!
!    FEM_SAMPLE_PRB tests routines from the FEM_SAMPLE library.
!
!  Modified:
!
!    12 October 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FEM_SAMPLE_PRB'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Test routines in the FEM_SAMPLE library.'

  call test01 ( 'nodes.txt', 'triangles.txt', 'coefs.txt' )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FEM_SAMPLE_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine test01 ( node_file, triangle_file, coef_file )

!*****************************************************************************80
!
!! 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.
!
  implicit none

  character ( len = * ) :: coef_file
  real ( kind = 8 ), allocatable, dimension ( : ) :: coef_node
  integer dim_num
  integer i
  integer j
  character ( len = * ) :: node_file
  integer node_num
  real ( kind = 8 ), allocatable, dimension ( :, : ) :: node_xy
  integer number
  integer triangle
  character ( len = * ) :: triangle_file
  integer, allocatable, dimension ( :, : ) :: triangle_neighbor
  integer, allocatable, dimension ( :, : ) :: triangle_node
  integer triangle_num
  integer triangle_order
  real ( kind = 8 ), allocatable, dimension ( :, : ) :: u_grid
  integer x_number
  real ( kind = 8 ), allocatable, dimension ( : ) :: x_vec
  integer y_number
  real ( kind = 8 ), allocatable, dimension ( : ) :: y_vec

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FEM_SAMPLE_TEST01:'
  write ( *, '(a)' ) '  Evaluate a finite element function on a uniform grid.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  This function expects to find three files to read:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  * node_file,      containing the node coordinates,'
  write ( *, '(a)' ) '  * triangle_file,  containing nodes that form triangles'
  write ( *, '(a)' ) '  * coef_file,      containing nodal coefficient values.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  It reads the files, constructs a grid of evenly spaced'
  write ( *, '(a)' ) '  sample points, and the evaluates the finite element'
  write ( *, '(a)' ) '  at the sample points.'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The node file is "' // trim ( node_file ) // '".'
  write ( *, '(a)' ) '  The triangle file is "' // trim ( triangle_file ) // '".'
  write ( *, '(a)' ) '  The coefficient file is "' // trim ( coef_file ) // '".'
!
!  Read the data.
!
  call dtable_header_read (  node_file, dim_num, node_num )

  allocate ( node_xy(1:dim_num,1:node_num) )

  call dtable_data_read ( node_file, dim_num, node_num, node_xy )

  write ( *, '(a,i8)' ) '  The number of nodes is ', node_num

  call itable_header_read ( triangle_file, triangle_order, triangle_num )

  allocate ( triangle_node(1:triangle_order,1:triangle_num) )
  allocate ( triangle_neighbor(1:3,1:triangle_num) )

  call itable_data_read ( triangle_file, triangle_order, triangle_num, &
    triangle_node )

  write ( *, '(a,i8)' ) '  The triangle order is ', triangle_order
  write ( *, '(a,i8)' ) '  The number of triangles is ', triangle_num

  allocate ( coef_node(1:node_num) )

  call dtable_data_read ( coef_file, 1, node_num, coef_node )

  write ( *, '(a)' ) '  The coefficient data has been read.'

  if ( triangle_order == 3 ) then

    call triangulation_order3_neighbor_triangles ( triangle_num, &
      triangle_node, triangle_neighbor )

  else if ( triangle_order == 6 ) then

    call triangulation_order6_neighbor_triangles ( triangle_num, &
      triangle_node, triangle_neighbor )

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FEM_SAMPLE_TEST01 - Fatal error!'
    write ( *, '(a)' ) '  The triangle order must be 3 or 6.'
    write ( *, '(a,i8)' ) '  But this data has triangle order = ', &
      triangle_order
    return

  end if

  write ( *, '(a)' ) '  The triangle neighbor array has been computed.'

  if ( .false. ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  TRIANGLE NEIGHBOR:'
    write ( *, '(a)' ) ' '
    do triangle = 1, triangle_num
      write ( *, '(2x,i8,2x,i8,2x,i8,2x,i8)' ) &
        triangle, triangle_neighbor(1:3,triangle)
    end do
  end if

  number = 4

  call grid_sample_size ( node_num, node_xy, number, x_number, y_number )

  allocate ( u_grid(x_number,y_number) )
  allocate ( x_vec(x_number) )
  allocate ( y_vec(y_number) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Requested   node density = ',   number
  write ( *, '(a,i8)' ) '  Suggested X node density = ', x_number
  write ( *, '(a,i8)' ) '  Suggested Y node density = ', y_number

  call 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 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  X sample coordinates:'
  write ( *, '(a)' ) ' '

  write ( *, '(5(2x,f10.4))' ) x_vec(1:x_number)
  
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Y sample coordinates:'
  write ( *, '(a)' ) ' '

  write ( *, '(5(2x,f10.4))' ) y_vec(1:y_number)
  
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  U evaluated at sample points:'
  write ( *, '(a)' ) ' '

  do j = y_number, 1, -1
    write ( *, '(5(2x,f10.4))' ) u_grid(1:x_number,j)
  end do

  deallocate ( coef_node )
  deallocate ( node_xy )
  deallocate ( triangle_neighbor )
  deallocate ( triangle_node )
  deallocate ( u_grid )
  deallocate ( x_vec )
  deallocate ( y_vec )

  return
end