function triangulation_l2q ( input_node_filename, input_triangulation_filename )

%% MAIN is the main program for TRIANGULATION_L2Q.
%
%  Discussion:
%
%    TRIANGULATION_L2Q makes a quadratic triangulation from a linear one.
%
%  Modified:
%
%    01 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Usage:
%
%    triangulation_l2q ( 'node_file', 'tri_file' )
%
  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_L2Q\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Read a "linear" triangulation and\n' );
  fprintf ( 1, '  write out a "quadratic" triangulation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read a dataset of NODE_NUM1 points in 2 dimensions.\n' );
  fprintf ( 1, '  Read an associated triangulation dataset of TRIANGLE_NUM \n' );
  fprintf ( 1, '  triangles which uses 3 nodes per triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Create new nodes which are triangle midpoints,\n' );
  fprintf ( 1, '  generate new node and triangulation data for\n' );
  fprintf ( 1, '  quadratic 6-node triangles, and write them out.\n' );
%
%  Get the number of command line arguments.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_L2Q:\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_L2Q:\n' );
    input_triangulation_filename = input ( ...
      '  Please enter the name of the triangulation file.' );

  end
%
%  Read the data.
%
  [ m, 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 M = %d\n', m );
  fprintf ( 1, '  Number of points NODE_NUM1  = %d\n', node_num1 );

  node_xy1(1:m,1:node_num1) = dtable_data_read ( input_node_filename, ...
    m, node_num1 );

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

  r8mat_transpose_print_some ( m, node_num1, node_xy1, 1, 1, 5, 5, ...
    '  5 by 5 portion of data read from file:' );
%
%  Read the triangulation data.
%
  [ m3, triangle_num ] = itable_header_read ( input_triangulation_filename );

  if ( m3 ~= 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_L2Q - Fatal error!\n' );
    fprintf ( 1, '  Data is not for a 3-node triangulation.\n' );
    error ( 'TRIANGULATION_L2Q - 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', m3 );
  fprintf ( 1, '  Number of triangles TRIANGLE_NUM  = %d\n', triangle_num );

  triangle_node1(1:m3,1:triangle_num) = itable_data_read ( ...
    input_triangulation_filename, m3, triangle_num );

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

  i4mat_transpose_print_some ( m3, triangle_num, triangle_node1, ...
    1, 1, 3, 10, '  3 by 10 portion TRIANGLE_NODE1:' );
%
%  Create the output file names from the input file names.
%
  output_node_filename = file_name_ext_swap ( input_node_filename, 'l2q.txt' );

  output_triangulation_filename = file_name_ext_swap ( ...
    input_triangulation_filename, 'l2q.txt' );
%
%  Determine the number of midside nodes that will be added.
%
  boundary_num = triangulation_order3_boundary_edge_count ( triangle_num, ...
    triangle_node1 );

  interior_num = ( 3 * triangle_num - boundary_num ) / 2;
  edge_num = interior_num + boundary_num;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of midside nodes to add = %d\n', edge_num );

  node_num2 = node_num1 + edge_num;
%
%  Build the triangle neighbor array.
%
  triangle_neighbor = triangulation_order3_neighbor_triangles ( ...
    triangle_num, triangle_node1 );

  i4mat_transpose_print ( 3, triangle_num, triangle_neighbor, ...
    '  Triangle_neighbor' );
%
%  Create the midside nodes.
%
  triangle_node2(1:3,1:triangle_num) = triangle_node1(1:3,1:triangle_num);
  triangle_node2(4:6,1:triangle_num) = -1;
  node_xy2(1:2,1:node_num1) = node_xy1(1:2,1:node_num1);

  node_num2 = node_num1;

  for triangle = 1 : triangle_num

    for i = 1 : 3

      triangle2 = triangle_neighbor(i,triangle);

      if ( 0 < triangle2 & triangle2 < triangle )
        continue
      end

      ip1 = i4_wrap ( i + 1, 1, 3 );

      k1 = triangle_node2(i,triangle);
      k2 = triangle_node2(ip1,triangle);

      node_num2 = node_num2 + 1;
      node_xy2(1:m,node_num2) = 0.5 ...
        * ( node_xy1(1:m,k1) + node_xy1(1:m,k2) );

      fprintf ( 1, '%8d  %14f  %14f\n', node_num2, node_xy2(1:m,node_num2) );

      triangle_node2(3+i,triangle) = node_num2;

      if ( 0 < triangle2 )
        for ii = 1 : 3
          if ( triangle_neighbor(ii,triangle2) == triangle )
            triangle_node2(ii+3,triangle2) = node_num2;
          end
        end
      end

    end
  end

  i4mat_transpose_print ( 6, triangle_num, triangle_node2, ...
    '  TRIANGLE_NODE2' );
%
%  Write out the node and triangle data for the quadratic mesh
%
  r8mat_transpose_print ( m, node_num2, node_xy2, '  NODE_XY2:' );

  header = 1;
  dtable_write ( output_node_filename, m, node_num2, node_xy2, header );

  header = 1;
  itable_write ( output_triangulation_filename, 6, triangle_num, ...
    triangle_node2, header );

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

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