function free_fem_heat ( node_file_name, element_file_name )

%% MAIN is the main routine of FREE_FEM_HEAT.
%
%  Discussion:
%
%    This program solves the heat equation
%
%      dUdT - Laplacian U(X,Y,T) + K(X,Y,T) * U(X,Y,T) = F(X,Y,T)
%
%    in a triangulated region in the plane.  
%
%    Along the boundary of the region, Dirichlet conditions 
%    are imposed:
%
%      U(X,Y,T) = G(X,Y,T)
%
%    At the initial time T_INIT, the value of U is given
%    at all points in the region:
%
%      U(X,Y,T_INIT) = H(X,Y)
%
%    The code uses continuous piecewise linear basis functions on 
%    triangles.
%
%    The backward Euler approximation is used for the time derivatives.
%
%  Problem specification:
%
%    The user defines the geometry by supplying two data files
%    which list the node coordinates, and list the nodes that make up
%    each element.
%
%    The user specifies the coefficient function K(X,Y,T) 
%    by supplying a routine of the form
%
%      k = k_coef ( node_num, node_xy, time )
%
%    The user specifies the right hand side
%    by supplying a routine of the form
%
%     f = rhs ( node_num, node_xy, time )
%
%    The user specifies the right hand side of the Dirichlet boundary
%    conditions by supplying a function
%
%      u = dirichlet_condition ( node_num, node_xy, time )
%
%    The user specifies the initial condition by supplying a function
%
%      u = initial_condition ( node_num, node_xy, time )
%
%  Usage:
%
%    free_fem_heat ( node_file, element_file )
%
%    invokes the program:
%
%    * "node_file" contains the coordinates of the nodes;
%    * "element_file" contains the indices of nodes that make up each element.
%
%    Files created include:
%
%    * "nodes.eps", an image of the nodes;
%    * "elements.eps", an image of the elements;
%    * "u0000.txt", the value of the solution at the initial condition.
%    * "u0001.txt" through "UNNNN.txt", the value of the solution at
%      each time step;
%    * "time.txt", the value of time at each step, from the initial to 
%      final times.
%
%  Modified:
%
%    08 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Local parameters:
%
%    Local, real A(3*IB+1,NODE_NUM), the coefficient matrix.
%
%    Local, integer DIM_NUM, the spatial dimension, which is 2.
%
%    Local, string ELEMENT_FILE_NAME, the name of the
%    input file containing the element information.
%
%    Local, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); 
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Local, integer ELEMENT_NUM, the number of elements.
%
%    Local, integer ELEMENT_ORDER, the order of each element.
%
%    Local, real F(NODE_NUM), the right hand side.
%
%    Local, integer IB, the half-bandwidth of the matrix.
%
%    Local, integer NODE_NUM, the number of nodes.
%
%    Local, logical NODE_BOUNDARY(NODE_NUM), is TRUE if a given
%    node is on the boundary.
%
%    Local, integer NODE_CONDITION(NODE_NUM), indicates the type of 
%    boundary condition being applied to nodes on the boundary.
%
%    Local, string NODE_FILE_NAME, the name of the
%    input file containing the node coordinate information.
%
%    Local, real NODE_U(NODE_NUM), the finite element coefficients.
%
%    Local, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Local, integer QUAD_NUM, the number of quadrature points used for 
%    assembly.  This is currently set to 3, the lowest reasonable value.
%    Legal values are 1, 3, 4, 6, 7, 9, 13, and for some problems, a value 
%    of QUAD_NUM greater than 3 may be appropriate.
%
  debug = 0;
  quad_num = 3;

  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FREE_FEM_HEAT\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution of the time dependent heat equation\n' );
  fprintf ( 1, '  on an arbitrary triangulated region D in 2 dimensions.n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ut - Uxx - Uyy + K(x,y,t) * U = F(x,y,t) in D;\n' );
  fprintf ( 1, '                              U = G(x,y,t) on boundary;\n' );
  fprintf ( 1, '                              U = H(x,y,t) at initial time.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The finite element method is used, with\n' );
  fprintf ( 1, '  6 node quadratic triangular elements ("T6").\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The time derivative is approximated using the\n' );
  fprintf ( 1, '  backward Euler formula.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current status:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * Time step information currently set internally!\n' );
  fprintf ( 1, '  * Would be easy to do linear triangles as well.\n' );
  fprintf ( 1, '  * Do you want ability to compare to an exact solution?\n' );
  fprintf ( 1, '\n' );
%
%  Get the file names.
%
  [ node_file_name, element_file_name ] = file_name_specification ( ...
    node_file_name, element_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node file is "%s".\n', node_file_name );
  fprintf ( 1, '  Element file is "%s".\n', element_file_name );
%
%  Read the node coordinate file.
%
  [ dim_num, node_num ] = dtable_header_read ( node_file_name );

  fprintf ( 1, '  Number of nodes =          %d\n', node_num );

  node_xy = dtable_data_read ( node_file_name, dim_num, node_num );

  r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, 2, 10, ...
    '  First 10 nodes' );
%
%  Read the element description file.
%
  [ element_order, element_num ] = itable_header_read ( element_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Element order =            %d\n', element_order );
  fprintf ( 1, '  Number of elements =       %d\n', element_num );

  if ( element_order ~= 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FREE_FEM_HEAT - Fatal error!\n' );
    fprintf ( 1, '  The input triangulation has order %d.\n', element_order );
    fprintf ( 1, '  However, a triangulation of order 6 is required.' );
    error ( 'FREE_FEM_POISSON - Fatal error!' );
  end

  element_node = itable_data_read ( element_file_name, element_order, ...
    element_num );

  i4mat_transpose_print_some ( element_order, element_num, ...
    element_node, 1, 1, element_order, 10, '  First 10 elements' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature order =          %d\n', quad_num );
%
%  Determine which nodes are boundary nodes and which have a
%  finite element unknown.  Then set the boundary values.
%
  node_boundary = triangulation_order6_boundary_node ( node_num, ...
    element_num, element_node );

  if ( debug )
    lvec_print ( node_num, node_boundary, '    Node  Boundary?' );
  end
%
%  Determine the node conditions.
%  For now, we'll just assume all boundary nodes are Dirichlet.
%
  node_condition(1:node_num) = 1;

  for node = 1 : node_num
    if ( node_boundary(node) )
      node_condition(node) = 2;
    end
  end
%
%  Determine the bandwidth of the coefficient matrix.
%
  ib = bandwidth ( element_order, element_num, element_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix half bandwidth is %d\n', ib );
  fprintf ( 1, '  The matrix bandwidth is      %d\n', 2 * ib + 1 );
  fprintf ( 1, '  The storage bandwidth is     %d\n', 3 * ib + 1 );
%
%  Make a picture of the nodes.
%
  if ( node_num <= 1000 )

    node_file_name = 'nodes.eps';

    if ( node_num <= 100 )
      node_label = 1;
    else
      node_label = 0;
    end

    points_plot ( node_file_name, node_num, node_xy, node_label );

  end
%
%  Make a picture of the elements.
%
  if ( node_num <= 1000 )

    triangulation_eps_file_name = 'elements.eps';

    if ( node_num <= 100 )
      node_show = 2;
    elseif ( node_num <= 250 )
      node_show = 1;
    else
      node_show = 0;
    end

    if ( element_num <= 100 )
      element_show = 2;
    else
      element_show = 1;
    end

    triangulation_order6_plot ( triangulation_eps_file_name, node_num, ...
      node_xy, element_num, element_node, node_show, element_show );

  end
%
%  Set time stepping quantities.
%
  time_init = 0.0;
  time_final = 0.5;
  time_step_num = 10;
  time_step_size = ( time_final - time_init ) / time_step_num;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial time = %f\n', time_init );
  fprintf ( 1, '  Final time =   %f\n', time_final );
  fprintf ( 1, '  Step size =    %f\n', time_step_size );
  fprintf ( 1, '  Number of steps = %d\n', time_step_num );
%
%  Initialize the names of the time and solution file.
%
  time_file_name = 'time.txt';
  solution_file_name = 'u0000.txt';
%
%  Set the value of U at the initial time.
%
  time = time_init;
  u_exact(1:node_num) = initial_condition ( node_num, node_xy, time );
  u(1:node_num) = u_exact(1:node_num);

  time_unit = fopen ( time_file_name, 'wt');
  fprintf ( time_unit, '%14f\n', time );

  solution_write ( node_num, u, solution_file_name, time );
%
%  Time looping.
%
  for time_step = 1 : time_step_num

    time_old = time;
    u_old(1:node_num) = u(1:node_num);

    time = ( ( time_step_num - time_step ) * time_init    ...
           + (                 time_step ) * time_final ) ...
           / ( time_step_num             );
%
%  Assemble the finite element coefficient matrix A and the right-hand side F.
%
  [ a, f ] = assemble_heat ( node_num, node_xy, node_condition, ...
    element_order, element_num, element_node, quad_num, ib, time );

  if ( debug ) 
    dgb_print_some ( node_num, node_num, ib, ib, a, 10, 1, 12, 25, ...
      '  Initial block of matrix A:' );

    r8vec_print_some ( node_num, f, 1, 10, '  Part of right hand side F:' );
  end
%
%  Adjust the linear system for the dU/dT term, which we are treating
%  using the backward Euler formula.
%
    [ a, f ] = assemble_backward_euler ( node_num, node_xy, element_order, ...
      element_num, element_node, quad_num, ib, time, time_step_size, ...
      u_old, a, f );

    if ( debug )

      dgb_print_some ( node_num, node_num, ib, ib, a, 10, 1, 12, 25, ...
        '  A after DT adjustment:' );

      r8vec_print_some ( node_num, f, 1, 10, ...
        '  F after DT adjustment:' );

    end
%
%  Adjust the linear system to account for boundary conditions.
%
    [ a, f ] = assemble_boundary ( node_num, node_xy, node_condition, ib, ...
      time, a, f );

    if ( debug )
 
      dgb_print_some ( node_num, node_num, ib, ib, a, 10, 1, 12, 25, ...
        '  A after BC adjustment:' );

      r8vec_print_some ( node_num, f, 1, 10, ...
        '  F after BC adjustment:' );

    end
%
%  Solve the linear system using a banded solver.
%
    [ a, pivot, ierr] = dgb_fa ( node_num, ib, ib, a );

    if ( ierr ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'FREE_FEM_HEAT - Fatal error!\n' );
      fprintf ( 1, '  DGB_FA returned an error condition.\n' );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  The linear system was not factored, and the\n' );
      fprintf ( 1, '  algorithm cannot proceed.\n' );
      error ( 'FREE_FEM_HEAT - Fatal error!' );
    end

    job = 0;
    node_u(1:node_num) = f(1:node_num);

    node_u = dgb_sl ( node_num, ib, ib, a, pivot, node_u, job );

    if ( debug ) 
      r8vec_print_some ( node_num, node_u, 1, 10, ...
        '  Part of the solution vector:' );
    end
%
%  Can compute errors versus exact solution here.
%

%
%  Increment the file name, and write the new solution.
%
    fprintf ( time_unit, '%14f\n', time );

    solution_file_name = file_name_inc ( solution_file_name );

    solution_write ( node_num, u, solution_file_name, time );

  end

  fclose ( time_unit );

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

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