function ffs_sparse (  node_file_name, element_file_name )

%% MAIN is the main routine of FFS_SPARSE.
%
%  Discussion:
%
%    This program is a variant of FREE_FEM_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.
%
%    Only three routines needed to be changed:
%    * the main program FREE_FEM_STOKES is replaced by FFS_SPARSE.
%    * the routine ASSEMBLE_STOKES is replaced by ASSEMBLE_STOKES_SPARSE.
%    * the routine DIRICHLET_APPLY is replaced by 
%      DIRICHLET_APPLY_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 Stokes equations 
%    for velocity vector V and scalar pressure P:
%
%      - nu * Laplacian W(X,Y) + 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, the user controls the type of
%    boundary condition to be imposed, if any.  Currently, these
%    conditions may be of Dirichlet form:
%
%      U(X,Y) = U_BC(X,Y)
%      V(X,Y) = V_BC(X,Y)
%      P(X,Y) = P_BC(X,Y)
%
%    or Neumann form with ZERO right hand side:
%
%      dU/dn(X,Y) = 0
%      dV/dn(X,Y) = 0
%      dP/dn(X,Y) = 0
%
%    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;
%
%      function nu = constants ( 'DUMMY' )
%
%    * the right hand side of the Stokes equations:
%
%      function [ u_rhs, v_rhs, p_rhs ] = rhs ( node_num, node_xy )
%
%    * the type of boundary conditions imposed:
%
%      [ node_u_condition, node_v_condition, node_p_condition ] =
%        boundary_type ( node_num, node_xy, node_boundary, node_type,
%        node_u_condition, node_v_condition, node_p_condition )
%
%    * the right hand side of any Dirichlet boundary conditions:
%
%      function [ u_bc, v_bc, p_bc ] = dirichlet_condition ( node_num, node_xy )
%
%    * the right hand side of any Neumann boundary conditions:
%      (currently, nonzero values will be ignored)
%
%      function [ u_bc, v_bc, p_bc ] = neumann_condition ( node_num, node_xy )
%
%  Usage:
%
%    ffs_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;
%    * "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 finite element coefficient matrix, in
%    MATLAB's sparse matrix storage.
%
%    Local, integer ELEMENT_NODE(6,ELEMENT_NUM), the nodes that form each 
%    element.  Nodes 1, 2, and 3 are the vertices.  Node 4 is between 1 
%    and 2, and so on.
%
%    Local, integer ELEMENT_NUM, the number of elements.
%
%    Local, integer ELEMENT_ORDER, the element order.
%
%    Local, real F(VARIABLE_NUM), the finite element right hand side.
%
%    Local, logical NODE_BOUNDARY(NODE_NUM), is TRUE if the node is
%    found to lie on the boundary of the region.
%
%    Local, real NODE_C(NODE_NUM), the finite element coefficients.
%
%    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 (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_P_VARIABLE(NODE_NUM), the index of the pressure 
%    variable associated with a node, or -1 if there is none.
%
%    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 (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_U_VARIABLE(NODE_NUM), the index of the horizontal
%    velocity variable associated with a node, or -1 if there is none.
%
%    Local, integer NODE_V_CONDITION(NODE_NUM), 
%    indicates the condition used to determine vertical velocity 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_V_VARIABLE(NODE_NUM), the index of the vertical
%    velocity variable associated with a node, or -1 if there is none.
%
%    Local, real NODE_XY(2,NODE_NUM), the coordinates of the nodes.
%
%    Local, integer NODE3_NUM, the number of order 3 nodes.
%
%    Local, integer NODE3_LABEL(NODE_NUM), contains the renumbered
%    label of order3 nodes, and -1 for nodes that are not order3 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 unknowns.
%
  debugging = 1;
  quad_num = 7;

  nu = constants ( 'DUMMY' );

  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FFS_SPARSE:\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A version of FREE_FEM_STOKES using MATLAB''s \n' );
  fprintf ( 1, '  sparse matrix storage, factor and solve facilities.\n' ); 
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Finite element solution of the \n' );
  fprintf ( 1, '  steady incompressible Stokes equations\n' );
  fprintf ( 1, '  on a triangulated region in 2 dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  - nu * ( Uxx + Uyy ) + dPdx = F1(x,y)\n' );
  fprintf ( 1, '  - nu * ( Vxx + Vyy ) + dPdy = F2(x,y)\n' );
  fprintf ( 1, '                      Ux + Vy = F3(x,y).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Boundary conditions may be of Dirichlet type:\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, '  or of Neumann type with zero right hand side:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    dU/dn(x,y) = 0\n' );
  fprintf ( 1, '    dV/dn(x,y) = 0\n' );
  fprintf ( 1, '    dP/dn(x,y) = 0\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, '  Quadrature order =       %d\n', quad_num );
  fprintf ( 1, '  Kinematic viscosity NU = %f\n', nu );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current status:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * testing zero Neumann condition option.\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 );

  if ( debugging )

    r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, 2, 10, ...
      '  First 10 nodes' );

  end
%
%  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, 'FFS_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 ( 'FFS_SPARSE - Fatal error!' );
  end

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

  if ( debugging )

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

  end
%
%  Determine the "type" of each node.
%  A vertex node, of type 1, has U, V, and P variables.
%  A midside node, of type 2, has U and V only.
%
  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 )
    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, 'FFS_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 ( 'FFS_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 );
%
%  Allow the user to examine and modify the tentative boundary conditions.
%
  [ node_u_condition, node_v_condition, node_p_condition ] = ...
    boundary_type ( node_num, node_xy, node_boundary, node_type, ...
    node_u_condition, node_v_condition, node_p_condition );

  neumann_num = 0;

  for node = 1 : node_num

    if ( node_u_condition(node) == 3 )
      neumann_num = neumann_num + 1;
    end

    if ( node_v_condition(node) == 3 )
      neumann_num = neumann_num + 1;
    end

    if ( node_p_condition(node) == 3 )
      neumann_num = neumann_num + 1;
    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of Neumann conditions added = \n', neumann_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Boundary conditions per node:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Node    U_cond	V_cond    P_cond\n' );
  fprintf ( 1, '\n' );
  for node = 1 : node_num

    if ( node <= 10 || node_num - 10 <= node )
      fprintf ( 1, '  %8d  %8d  %8d  %8d\n', ...
        node, node_u_condition(node), node_v_condition(node), ...
        node_p_condition(node) );
    end

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

  end
%
%  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 = 'ffs_sparse_nodes.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 = 'ffs_sparse_elements.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 );
%
%  Assemble the finite element coefficient matrix A and the right-hand side F.
%
  [ 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, ...
      '  Initial block of Stokes stiffness matrix A:' );

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

  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, ...
      '  Matrix A after Dirichlet BC adjustments:' );

    r8vec_print_some ( variable_num, f, 1, 10, ...
      '  Part of right hand side after Dirichlet BC adjustments:' );

  end
%
%  Adjust the linear system to account for Neumman boundary conditions.
%
  f = neumann_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, f );

  if ( 0 ) 

    r8mat_print_some ( variable_num, variable_num, a, 1, 1, 10, 10, ...
      '  Matrix A after Neumann BC adjustments:' );

    r8vec_print_some ( variable_num, f, 1, 10, ...
      '  Part of right hand side after Neumann BC adjustments:' );

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

  if ( debugging )

    r8vec_print_some ( variable_num, node_c, 1, 10, ...
      '  Part of the solution vector:' );

  end
%
%  Print the solution vector based at nodes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Variable indices per node:\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
%
%  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.
%
  file_name = 'pressure3.txt';

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

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pressures written to "%s".\n', file_name );
%
%  Write an ASCII file that can be read into MATLAB.
%
  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, 'FFS_SPARSE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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