function triangulation_refine ( input_node_filename, ...
  input_triangulation_filename )

%% MAIN is the main program for TRIANGULATION_REFINE.
%
%  Discussion:
%
%    TRIANGULATION_REFINE refines a triangulation by doubling.
%
%  Modified:
%
%    10 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Usage:
%
%    triangulation_refine ( 'node_file', 'triangle_file' )
%
  debug = 1;
  
  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_REFINE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Read a "linear" or "quadratic" triangulation\n' );
  fprintf ( 1, '  and write out a refined triangulation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read a dataset of NODE_NUM1 points in 2 dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read an associated triangulation dataset of TRIANGLE_NUM \n' );
  fprintf ( 1, '  triangles which uses 3 or 6 nodes per triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Subdivide each triangle into 4 triangles,\n' );
  fprintf ( 1, '  generate new nodes as midpoints of current nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Write out the new node and triangulation data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  If the input triangulation was Delaunay, then\n' );
  fprintf ( 1, '  the output triangulation will be Delaunay.\n' );
%
%  Get the number of command line arguments.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_REFINE:\n' );
    input_node_filename = input ( '  Please enter the name of the node file.' );

  end
%
%  If at least two command line arguments, the second is the triangulation file.
%
  if ( nargin < 2 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_REFINE:\n' );
    input_triangulation_filename = input ( ...
      '  Please enter the name of the triangulation file.' );

  end
%
%  Read the node data.
%
  [ dim_num, node_num1 ] = dtable_header_read (  input_node_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', input_node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of points NODE_NUM1  = %d\n', node_num1 );

  node_xy1 = dtable_data_read ( input_node_filename, ...
    dim_num, node_num1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', input_node_filename );

  r8mat_transpose_print_some ( dim_num, node_num1, node_xy1, 1, 1, ...
    dim_num, 5, '  First 5 nodes:' );
%
%  Read the triangulation data.
%
  [ triangle_order, triangle_num1 ] = ...
    itable_header_read ( input_triangulation_filename );

  if ( triangle_order ~= 3 & triangle_order ~= 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_REFINE - Fatal error!\n' );
    fprintf ( 1, '  Data is not for a 3 or 6-node triangulation.\n' );
    error ( 'TRIANGULATION_REFINE - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of ""%s".\n', input_triangulation_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Triangle order = %d\n', triangle_order );
  fprintf ( 1, '  Number of triangles TRIANGLE_NUM1  = %d\n', triangle_num1 );

  triangle_node1 = itable_data_read ( ...
    input_triangulation_filename, triangle_order, triangle_num1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in ""%s".\n', input_triangulation_filename );

  i4mat_transpose_print_some ( triangle_order, triangle_num1, triangle_node1, ...
    1, 1, triangle_order, 5, '  First 5 triangles:' );
%
%  Create the output file names from the input file names.
%
  output_node_filename = file_name_ext_swap ( ...
    input_node_filename, 'ref.txt' );

  output_triangulation_filename = file_name_ext_swap ( ...
    input_triangulation_filename, 'ref.txt' );
%
%  Determine the size of the refined mesh.
%
  if ( triangle_order == 3 )

    [ node_num2, triangle_num2, edge_data ] = ...
      triangulation_order3_refine_size ( node_num1, triangle_num1, ...
      triangle_node1 );

  elseif ( triangle_order == 6 )

    [ node_num2, triangle_num2, edge_data ] = ...
      triangulation_order6_refine_size ( node_num1, triangle_num1, ...
      triangle_node1 );

  end 

  fprintf ( 1, '\n' );
  fprintf ( 1,'  Number of nodes in refined mesh =      %d\n', node_num2 );
  fprintf ( 1,'  Number of triangles in refined mesh =  %d\n', triangle_num2 );
%
%  Compute the refined mesh.
%
  if ( triangle_order == 3 )

    [ node_xy2, triangle_node2 ] = triangulation_order3_refine_compute ( ...
      node_num1, triangle_num1, node_xy1, triangle_node1, node_num2, ...
      triangle_num2, edge_data );

  else

    [ node_xy2, triangle_node2 ] = triangulation_order6_refine_compute ( ...
      node_num1, triangle_num1, node_xy1, triangle_node1, node_num2, ...
      triangle_num2, edge_data );

  end

  if ( debug )

    r8mat_transpose_print_some ( dim_num, node_num2, node_xy2, ...
      1, 1, dim_num, 5, '  First 5 output nodes:' );

    i4mat_transpose_print_some ( triangle_order, triangle_num2, ...
      triangle_node2, 1, 1, triangle_order, 5, '  First 5 output triangles:' );

  end
%
%  Write out the node and triangle data.
%
  header = 0;

  dtable_write ( output_node_filename, dim_num, node_num2, node_xy2, header );

  itable_write ( output_triangulation_filename, triangle_order, ...
    triangle_num2, triangle_node2, header );
%
%  Free up memory.
%
  edge_data = [];
  node_xy1 = [];
  node_xy2 = [];
  triangle_node1 = [];
  triangle_node2 = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_REFINE\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp;