function free_fem_poisson ( node_file_name, element_file_name )

%% MAIN is the main routine for FREE_FEM_POISSON.
%
%  Discussion:
%
%    This program solves the Poisson equation
%
%      -Laplacian U(X,Y) + K(X,Y) * U = F(X,Y)
%
%    in a triangulated region in the plane.
%
%    Along the boundary of the region, Dirichlet conditions
%    are imposed:
%
%      U(X,Y) = G(X,Y)
%
%    The code uses continuous piecewise linear basis functions on
%    triangles.
%
%  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 right hand side of the Dirichlet boundary
%    conditions by supplying a function
%
%      function node_bc = dirichlet_condition ( node_num, node_xy )
%
%    The user specifies the coefficient function K(X,Y) of the Poisson
%    equation by supplying a routine of the form
%
%      function node_k = k_coef ( node_num, node_xy )
%
%    The user specifies the right hand side of the Poisson equation
%    by supplying a routine of the form
%
%      function node_rhs = rhs ( node_num, node_xy )
%
%  Usage:
%
%    free_fem_poisson ( 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;
%    * "solution.txt", the value of the solution at every node.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Local parameters:
%
%    Local, real A(3*IB+1,NODE_NUM), the coefficient matrix.
%
%    Local, real EH1, the H1 seminorm error.
%
%    Local, real EL2, the L2 error.
%
%    Local, integer ELEMENT_NODE[3*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 element order.
%
%    Local, real F(NODE_NUM), the right hand side.
%
%    Local, integer IB, the half-bandwidth of the matrix.
%
%    Local, logical NODE_BOUNDARY(NODE_NUM), is TRUE if the node is
%    found to lie on the boundary of the region.
%
%    Local, integer NODE_CONDITION(NODE_NUM), 
%    indicates the condition used to determine the variable at a node.
%    0, there is no condition (and no variable) at this node.
%    1, a finite element equation is used;
%    2, a Dirichlet condition is used. 
%    3, a Neumann condition is used. 
%
%    Local, integer NODE_NUM, the number of nodes.
%
%    Local, real NODE_R(NODE_NUM), the residual error.
%
%    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_POISSON:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Finite element solution of the \n' );
  fprintf ( 1, '  steady Poisson equation on a triangulated region\n' );
  fprintf ( 1, '  in 2 dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  - Uxx - Uyy + K(x,y) * U = F(x,y) in the region\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                    U(x,y) = G(x,y) on the boundary.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The finite element method is used,\n' );
  fprintf ( 1, '  with triangular elements,\n' );
  fprintf ( 1, '  which must be a 3 node linear triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current status:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * Boundary conditions cannot be set easily.\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 ~= 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FREE_FEM_POISSON - Fatal error!\n' );
    fprintf ( 1, '  The input triangulation has order %d.\n', element_order );
    fprintf ( 1, '  However, a triangulation of order 3 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 ( 3, element_num, ...
    element_node, 1, 1, 3, 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_order3_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_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 )
      triangle_show = 2;
    else
      triangle_show = 1;
    end

    triangulation_order3_plot ( triangulation_eps_file_name, node_num, ...
      node_xy, element_num, element_node, node_show, triangle_show );

  end
%
%  Assemble the finite element coefficient matrix A and the right-hand side F.
%
  [ a, f ] = assemble_poisson ( node_num, node_xy, ...
    element_num, element_node, quad_num, ib );
%
%  Print a tiny portion of the matrix.
%
  dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, 10, 10, ...
    '  Initial block of finite element matrix A:' );

  r8vec_print_some ( node_num, f, 1, 10, '  Part of right hand side vector:' );
%
%  Adjust the linear system to account for Dirichlet boundary conditions.
%
  [ a, f ] = dirichlet_apply ( node_num, node_xy, node_condition, ib, a, f );

  dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, 10, 10, ...
    '  Matrix A after boundary condition adjustments:' );

  r8vec_print_some ( node_num, f, 1, 10, '  Part of right hand side vector:' );
%
%  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_POISSON - 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_POISSON - 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 );

  r8vec_print_some ( node_num, node_u, 1, 10, ...
    '  Part of the solution vector:' );

  node_r = residual_poisson ( node_num, node_xy, node_condition, ...
    element_num, element_node, quad_num, ib, a, f, node_u );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum absolute residual = %f\n',  ...
    max ( abs ( node_r(1:node_num) ) ) );
%
%  Write an ASCII file that can be read into MATLAB.
%
  solution_file_name = 'solution.txt';

  solution_write ( node_num, node_u, solution_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FREE_FEM_POISSON:\n' );
  fprintf ( 1, '  Wrote an ASCII file\n' );
  fprintf ( 1, '    "%s"\n', solution_file_name );
  fprintf ( 1, '  of the form\n' );
  fprintf ( 1, '    U ( X(I), Y(I) )\n' );
  fprintf ( 1, '  which can be used for plotting.\n' );

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

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

  return
end