function ffns_sparse (  node_file_name, element_file_name )

%% MAIN is the main routine of FFNS_SPARSE.
%
%  Discussion:
%
%    This program is a variant of FREE_FEM_NAVIER_STOKES.  That program is
%    particularly limited because of its use of banded matrix storage and
%    solving routines.
%
%    This program discards the banded approach.  Instead, it uses MATLAB's 
%    sparse matrix storage, factorization and solution facilities, 
%    which allow this program to solve larger problems faster.
%
%    A few routines needed to be changed:
%    * the main program FREE_FEM_NAVIER_STOKES is replaced by FFNS_SPARSE.
%    * the routine ASSEMBLE_STOKES is replaced by ASSEMBLE_STOKES_SPARSE.
%    * the routine DIRICHLET_APPLY is replaced by 
%      DIRICHLET_APPLY_SPARSE.
%    * the routine JACOBIAN_FEM is replaced by JACOBIAN_FEM_SPARSE.
%    * the routine JACOBIAN_ADJUST_DIRICHLET is replaced by 
%      JACOBIAN_ADJUST_DIRICHLET_SPARSE.
%
%    This also means that the following routines are NOT needed:
%    * BANDWIDTH
%    * DGB_FA
%    * DGB_SL.
%    * DGB_PRINT_SOME.
%
%
%    This program solves the steady incompressible Navier Stokes equations for 
%    velocity vector W and scalar pressure P:
%
%      -nu * Laplacian W(X,Y) + W dot Grad W + Grad P(X,Y) = F(X,Y)
%
%                                               Div W(X,Y) = G(X,Y)
%
%    in an arbitrary triangulated region in the plane.  
%
%    Let U and V denote the scalar components of the velocity vector W.
%
%    Along the boundary of the region, Dirichlet conditions 
%    are imposed:
%
%      U(X,Y) = U_BC(X,Y)
%      V(X,Y) = V_BC(X,Y)
%      P(X,Y) = P_BC(X,Y)
%
%    The code uses the finite element method.  The Taylor-Hood element
%    is used, in which a single reference triangle is used to define
%    both a piecewise quadratic representation of velocity, and a piecewise
%    linear representation of pressure.
%
%  Geometry 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.
%
%  Equation specification:
%
%    The user specifies 
%
%    * the kinematic viscosity NU;
%
%      nu = constants ( 'DUMMY' )
%
%    * the right hand sides of the Navier-Stokes equations
%      by supplying a routine of the form
%
%      function [ u_rhs, v_rhs, p_rhs ] = rhs ( xy )
%
%    * the right hand sides of the Dirichlet boundary
%      conditions by supplying a function
%
%      function [ u_bc, v_bc, p_bc ] = dirichlet_condition ( xy )
%
%  Usage:
%
%    ffns_sparse ( 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.
%
%    Graphics files created include:
%
%    * "nodes6.eps", an image of the nodes;
%    * "triangles6.eps", an image of the elements;
%
%    Data files created include:
%
%    * "nodes3.txt", the nodes associated with pressure;
%    * "triangles3.txt", the linear triangles associated with pressure;
%    * "stokes_pressure3.txt", the Stokes pressure at the pressure nodes;
%    * "stokes_velocity6.txt", the Stokes velocity at the velocity nodes.
%    * "pressure3.txt", the pressure at the pressure nodes;
%    * "velocity6.txt", the velocity at the velocity nodes;
%
%  Modified:
%
%    26 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILE_NAME, the name of the node file.  If
%    this argument is not supplied, it will be requested.
%
%    Input, string ELEMENT_FILE_NAME, the name of the element file.  If
%    this argument is not supplied, it will be requested.
%
%  Local parameters:
%
%    Local, real sparse A, the Stokes stiffness matrix or Navier Stokes
%    jacobian matrix, stored as a MATLAB sparse matrix.
%
%    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 element order.
%
%    Local, real F(VARIABLE_NUM), the right hand side.
%
%    Local, integer IT_MAX, the maximum number of Newton iterations allowed.
%    0 is a legal value.  It simply solves the Stokes problem, computes
%    the Navier-Stokes residual, and stops.
%    1 does the above, and then takes one Newton step, and so on.
%
%    Local, real NODE_C(VARIABLE_NUM), the finite element
%    coefficients of the current solution estimate.
%
%    Local, real NODE_C_DEL(VARIABLE_NUM), the correction to
%    the finite element coefficients of the current solution estimate.
%
%    Local, real NODE_C_OLD(VARIABLE_NUM), the finite element
%    coefficients computed by the prefious step of the iteration.
%
%    Local, integer NODE_NUM, the number of nodes.
%
%    Local, integer NODE_P_CONDITION(NODE_NUM), 
%    indicates the condition used to determine pressure at a node.
%    0, there is no condition 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_P_VARIABLE(NODE_NUM), 
%    is the index of the pressure variable associated with the node,
%    or -1 if there is no associated pressure variable.
%
%    Local, integer NODE_TYPE(NODE_NUM), determines if the node is a 
%    vertex or midside node.
%    1, the node is a vertex (P, U, V variables are associated with it).
%    2, the node is a midside node (only U and V variables are associated.)
%
%    Local, integer NODE_U_CONDITION(NODE_NUM), 
%    indicates the condition used to determine horizontal velocity at a node.
%    0, there is no condition 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_U_VARIABLE(NODE_NUM), 
%    is the index of the horizontal velocity variable associated with the node.
%
%    Local, integer NODE_V_CONDITION(NODE_NUM), 
%    indicates the condition used to determine vertical velocity at a node.
%    0, there is no condition 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_V_VARIABLE(NODE_NUM), 
%    is the index of the vertical velocity variable associated with the node.
%
%    Local, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Local, integer NODE3_NUM, the number of pressure nodes.
%
%    Local, integer NODE3_LABEL(NODE_NUM), contains the renumbered
%    labels of pressures nodes, and -1 for nodes that are not pressure nodes.
%
%    Local, real NU, the kinematic viscosity.
%
%    Local, integer NZ_NUM, the maximum number of nonzero entries that
%    might occur in the matrix.  This is only needed so that MATLAB can
%    efficiently allocate space for the matrix.
%
%    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.
%
%    Local, integer VARIABLE_NUM, the number of variables.
%
  debugging = 0;
  it_max = 5;
  quad_num = 7;
  res_tol = 0.0001;

  nu = constants ( 'DUMMY' );

  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FFNS_SPARSE\n' );
  fprintf ( 1, '  MATLAB version):\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This is a special version of FREE_FEM_NAVIER_STOKES\n' );
  fprintf ( 1, '  which uses MATLAB sparse matrix storage, factorization,\n' );
  fprintf ( 1, '  and solution methods, replacing the general banded approach.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Finite element solution of the \n' );
  fprintf ( 1, '  steady incompressible Navier Stokes equations\n' );
  fprintf ( 1, '  on a triangulated region in 2 dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  - nu * ( Uxx + Uyy ) + UUx + VUy + dPdx = F1(x,y)\n' );
  fprintf ( 1, '  - nu * ( Vxx + Vyy ) + UVx + VVy + dPdy = F2(x,y)\n' );
  fprintf ( 1, '                          Ux +  Vy        = F3(x,y).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Boundary conditions:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       U(x,y) = U_BC(x,y)\n' );
  fprintf ( 1, '       V(x,y) = V_BC(x,y)\n' );
  fprintf ( 1, '       P(x,y) = P_BC(x,y)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The finite element method uses Taylor-Hood\n' );
  fprintf ( 1, '  triangular elements which are linear for pressure\n' );
  fprintf ( 1, '  and quadratic for velocity.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum number of Newton iterations IT_MAX = %d:\n', it_max );
  fprintf ( 1, '  Quadrature order =       %d\n', quad_num );
  fprintf ( 1, '  Kinematic viscosity NU = %f\n', nu );
  fprintf ( 1, '\n' );
%
%  Get the file names, if not specified on command line.
%
  if ( nargin < 2 )
    element_file_name = [];
  end
  
  if ( nargin < 1 )
    node_file_name = [];
  end
  
  [ 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 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, 'FFNS_SPARSE - 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 ( 'FFNS_SPARSE - Fatal error!' );
  end

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

  i4mat_transpose_print_some ( 6, element_num, ...
    element_node, 1, 1, 6, 10, '  First 10 elements' );
%
%  Determine the node type (vertex or midside):
%
  node_type(1:node_num) = 1;

  for element = 1 : element_num
    for j = 4 : 6
      node = element_node(j,element);
      node_type(node) = 2;
    end
  end
%
%  Determine which nodes are boundary nodes.
%
  node_boundary = triangulation_order6_boundary_node ( node_num, ...
    element_num, element_node );

  if ( 0 == 1 )
    lvec_print ( node_num, node_boundary, '    Node  Boundary?' );
  end
%
%  Determine the node conditions:
%
%  All conditions begin as finite element conditions.
%
  node_p_condition(1:node_num) = 1;
  node_u_condition(1:node_num) = 1;
  node_v_condition(1:node_num) = 1;
%
%  Conditions on velocities associated with a boundary node are Dirichlet
%  conditions.
%
  for node = 1 : node_num
    if ( node_boundary(node) )
      node_u_condition(node) = 2;
      node_v_condition(node) = 2;
    end
  end
%
%  Midside nodes have no associated pressure variable.
%
  for node = 1 : node_num
    if ( node_type(node) == 2 )
      node_p_condition(node) = 0;
    end
  end
%
%  Replace a single finite element pressure condition by a Dirichlet
%  condition.
%
  p_node = -1;
  for node = 1 : node_num
    if ( node_p_condition(node) == 1 )
      node_p_condition(node) = 2;
      p_node = node;
      break
    end
  end

  if ( p_node == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FFNS_SPARSE - Fatal error!\n' );
    fprintf ( 1, '  Unable to find a finite element pressure condition\n' );
    fprintf ( 1, '  suitable for replacement by a Dirichlet condition.\n' );
    error ( 'FFNS_SPARSE - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Dirichlet boundary condition on pressure\n' );
  fprintf ( 1, '  will be applied at node %d\n', p_node );
%
%  Number the variables.
%
  variable_num = 0;

  for node = 1 : node_num

    variable_num = variable_num + 1;
    node_u_variable(node) = variable_num;

    variable_num = variable_num + 1;
    node_v_variable(node) = variable_num;

    if ( node_type(node) == 1 )
      variable_num = variable_num + 1;
      node_p_variable(node) = variable_num;
    else
      node_p_variable(node) = -1;
    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Total number of variables is %d\n', variable_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Variable indices per node:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Node  U_index  V_index  P_index\n' );
  fprintf ( 1, '\n' );
  for node = 1 : node_num
    if ( node <= 10 || node_num - 10 <= node )
      fprintf ( 1, '  %6d  %6d  %6d  %6d\n', ...
        node, node_u_variable(node), node_v_variable(node), ...
        node_p_variable(node) );
    end

    if ( node == 10 )
      fprintf ( 1, '(SKIPPING ENTRIES)\n' );
    end

  end
%
%  Plot the nodes.
%
  if ( node_num <= 100 )

    file_name = 'nodes6.eps';
    node_label = 1;

    points_plot ( file_name, node_num, node_xy, node_label );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Order 6 nodes plotted in "%s".\n', file_name );

  end
%
%  Plot the elements.
%
  if ( node_num <= 100 )

    file_name = 'triangles6.eps';
    node_show = 2;
    triangle_show = 2;

    triangulation_order6_plot ( file_name, node_num, ...
      node_xy, element_num, element_node, node_show, triangle_show );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Order 6 triangles plotted in "%s".\n', ...
      file_name );

  end
%
%  Determine the element neighbor array, just so we can estimate
%  the nonzeros.
%
  element_neighbor = triangulation_order6_neighbor_triangles ( ...
    element_num, element_node );
%
%  Determine the maximum number of nonzeros.
%
  nz_num = ns_adj_count ( node_num, element_num, variable_num, ...
    element_node, element_neighbor, node_u_variable, node_v_variable, ...
    node_p_variable );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NS_ADJ_COUNT returns NZ_NUM = %d\n', nz_num );
%
%  Get an initial condition, by assembling the Stokes coefficient matrix A 
%  and the right-hand side F, and solving.
%
  [ a, f ] = assemble_stokes_sparse ( node_num, element_num, quad_num, ...
    variable_num, node_xy, node_p_variable, node_u_variable, ...
    node_v_variable, element_node, nu, nz_num );

  if ( debugging )

    r8mat_print_some ( variable_num, variable_num, a, 1, 1, 10, 10, ...
      '  Part of Stokes matrix :' );

    r8vec_print_some ( variable_num, f, 1, 10, ...
      '  Part of Stokes right hand side:' );

  end
%
%  Adjust the linear system to account for Dirichlet boundary conditions.
%
  [ a, f ] = dirichlet_apply_sparse ( node_num, node_xy, node_p_variable, ...
    node_u_variable, node_v_variable, node_p_condition, ...
    node_u_condition, node_v_condition, variable_num, a, f );

  if ( debugging )

    r8mat_print_some ( variable_num, variable_num, a, 1, 1, 10, 10, ...
      '  Part of Stokes matrix, adjusted for Dirichlet BC:' );

    r8vec_print_some ( variable_num, f, 1, 10, ...
      '  Part of Stokes right hand side, adjusted for Dirichlet BC:' );

  end
%
%  Solve the linear system using a banded solver.
%
  node_c = a \ f';

  if ( debugging )
    r8vec_print_some ( variable_num, node_c, 1, 10, ...
      '  Part of the solution vector:' );
  end
%
%  Print the Stokes solution vector based at nodes.
%
  if ( debugging )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Solution to the STOKES equations:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '    Node        U        V        P\n' );
    fprintf ( 1, '\n' );

    for node = 1 : node_num

      if ( node <= 10 || node_num - 10 <= node )
 
        iu = node_u_variable(node);
        iv = node_v_variable(node);
        ip = node_p_variable(node);

        if ( 0 < ip )
          fprintf ( 1, '  %6d  %14f  %14f  %14f\n', ...
            node, node_c(iu), node_c(iv), node_c(ip) );
        else
          fprintf ( 1, '  %6d  %14f  %14f\n', ...
            node, node_c(iu), node_c(iv) );
        end

      end

      if ( node == 10 )
        fprintf ( 1, '(SKIPPING ENTRIES)\n' );
      end

    end 

  end
%
%  Compute a renumbering of the pressure nodes.
%
  node3_num = 0;

  for node = 1 : node_num
    if ( node_type(node) == 1 )
      node3_num = node3_num + 1;
    end
  end

  node3_num = 0;

  for node = 1 : node_num
    if ( node_type(node) == 1 )
      node3_num = node3_num + 1;
      node3_label(node) = node3_num;
    else
      node3_label(node) = -1;
    end
  end
%
%  Write the pressure nodes to a file.
%
  file_name = 'nodes3.txt';

  nodes3_write ( file_name, node_num, node_xy, node_type );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pressure nodes written to "%s".\n', file_name );
%
%  Write the pressure triangles to a file.
%
  file_name = 'triangles3.txt';

  triangles3_write ( file_name, element_num, element_node, ...
    node_num, node3_label );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pressure triangles written to "%s".\n', file_name );
%
%  Write the pressures to a file.
%
  if ( 0 )

    file_name = 'stokes_pressure3.txt';

    pressure3_write ( file_name, node_num, node_p_variable, ...
      variable_num, node_c );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Stokes pressures written to "%s".\n', file_name );

  end
%
%  Write an ASCII file that can be read into MATLAB.
%
  if ( 0 )

    file_name = 'stokes_velocity6.txt';

    velocity6_write ( file_name, node_num, node_u_variable, ...
      node_v_variable, variable_num, node_c );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Stokes velocities written to "%s".\n', file_name );

  end
%
%  The Stokes solution provides an initial guess for the Newton iteration
%  applied to the Navier Stokes equations.
%
  it_num = 0;

  while ( 1 )

    node_r = residual_fem ( node_num, node_xy, element_num, ...
      element_node, quad_num, node_u_variable, ...
      node_v_variable, node_p_variable, variable_num, nu, node_c );

    if ( debugging )
      r8vec_print_some ( variable_num, node_r, 1, 10, ...
        '  Part of Navier-Stokes FEM residual:' );
    end

    res_norm = sqrt ( sum ( node_r(1:node_num).^2 ) );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  l2-norm of FEM residual = %f\n', res_norm );

    node_r = residual_adjust_dirichlet ( node_num, node_xy, node_p_variable, ...
      node_u_variable, node_v_variable, node_p_condition, ...
      node_u_condition, node_v_condition, variable_num, node_c, node_r );

    if ( debugging )
      r8vec_print_some ( variable_num, node_r, 1, 10, ...
        '  Part of Navier-Stokes FEM residual adjusted for BC:' );
    end

    node_r_norm = sqrt ( sum ( node_r(1:node_num).^2 ) );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  l2-norm of adjusted FEM residual = %f\n', node_r_norm );

    if ( node_r_norm <= res_tol )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Convergence.\n' );
      break;
    end

    if ( it_max <= it_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  The iteration limit has been reached.\n' );
      break
    end

    it_num = it_num + 1;
 
    a = jacobian_fem_sparse ( node_num, node_xy, element_num, ...
      element_node, quad_num, node_u_variable, ...
      node_v_variable, node_p_variable, variable_num, nu, node_c, nz_num ); 
%
%  Adjust the jacobian for boundary conditions.
%
    a = jacobian_adjust_dirichlet_sparse ( node_num, node_xy, ...
      node_p_variable, node_u_variable, node_v_variable, ...
      node_p_condition, node_u_condition, node_v_condition, ...
      variable_num, a );
%
%  Set up and solve the Newton system J * dX = - F(x)
%
    node_c_del = - a \ node_r';

    node_c_old(1:variable_num) = node_c(1:variable_num);
    node_c(1:variable_num) = node_c(1:variable_num) + node_c_del(1:variable_num);

    if ( debugging )
      r8vec_print_some ( variable_num, node_c_del, 1, 10, ...
        '  Part of Newton correction vector:' );
    end

    node_c_del_norm = sqrt ( sum ( node_c_del(1:variable_num).^2 ) );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  l2-norm of Newton correction = %f\n', node_c_del_norm );

  end
%
%  Print the Navier Stokes solution vector based at nodes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution to the NAVIER STOKES equations:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Node        U        V        P\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num

    if ( node <= 10 || node_num - 10 <= node )
 
      iu = node_u_variable(node);
      iv = node_v_variable(node);
      ip = node_p_variable(node);

      if ( 0 < ip )
        fprintf ( 1, '  %6d  %14f  %14f  %14f\n', ...
         node, node_c(iu), node_c(iv), node_c(ip) );
      else
        fprintf ( 1, '  %6d  %14f  %14f\n', ...
         node, node_c(iu), node_c(iv) );
      end

    end

    if ( node == 10 )
      fprintf ( 1, '(SKIPPING ENTRIES)\n' );
    end

  end
%
%  Write the pressures to a file.
%
  file_name = 'pressure3.txt';
  
  pressure3_write ( file_name, node_num, node_p_variable, ...
    variable_num, node_c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Navier Stokes pressures written to "%s".\n', file_name );
%
%  Write the velocities to a file.
%
  file_name = 'velocity6.txt';

  velocity6_write ( file_name, node_num, node_u_variable, ...
    node_v_variable, variable_num, node_c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Navier Stokes velocities written to "%s".\n', file_name );

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

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