function free_fem_stokes (  node_file_name, element_file_name )

%*****************************************************************************80
%
%% MAIN is the main routine of FREE_FEM_STOKES.
%
%  Discussion:
%
%    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:
%
%    free_fem_stokes ( 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.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  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 A(3*IB+1,VARIABLE_NUM), the VARIABLE_NUM by VARIABLE_NUM
%    finite element coefficient matrix, stored in a banded format.
%
%    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, 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, 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 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, 'FREE_FEM_STOKES\n' );
  fprintf ( 1, '  MATLAB version:\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 );
%
%  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, 'FREE_FEM_STOKES - 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_STOKES - 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 "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, 'FREE_FEM_STOKES - 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 ( 'FREE_FEM_STOKES - 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
    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', ...
      node, node_u_condition(node), node_v_condition(node), ...
      node_p_condition(node) );
  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
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', node, node_u_variable(node), ...
      node_v_variable(node), node_p_variable(node) );
  end
%
%  Determine the bandwidth of the coefficient matrix.
%
  ib = bandwidth ( element_num, element_node, ...
    node_num, node_p_variable, node_u_variable, node_v_variable );

  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 );
%
%  Plot the nodes.
%
  if ( node_num <= 100 )

    file_name = 'free_fem_stokes_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 = 'free_fem_stokes_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
%
%  Assemble the finite element coefficient matrix A and the right-hand side F.
%
  [ a, f ] = assemble_stokes ( node_num, element_num, quad_num, ...
    variable_num, node_xy, node_p_variable, node_u_variable, ...
    node_v_variable, element_node, nu, ib );

  if ( debugging )

    dgb_print_some ( variable_num, variable_num, ib, ib, a, 1, 1, 10, 10, ...
      '  Initial block of finite element 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 ( node_num, node_xy, node_p_variable, ...
    node_u_variable, node_v_variable, node_p_condition, ...
    node_u_condition, node_v_condition, variable_num, ib, a, f );

  if ( debugging )

    dgb_print_some ( variable_num, variable_num, ib, ib, 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 adjustment:' );

  end
%
%  Adjust the linear system to account for Neumann boundary conditions.
%
  f = neumann_apply ( 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 )

    dgb_print_some ( variable_num, variable_num, ib, ib, 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 adjustment:' );

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

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FREE_FEM_STOKES - 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_STOKES - Fatal error!' );
  end

  job = 0;
  node_c = dgb_sl ( variable_num, ib, ib, a_lu, pivot, f, job );

  r8vec_print_some ( variable_num, node_c, 1, 10, ...
    '  Part of the solution vector:' );
%
%  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

    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
%
%  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, 'FREE_FEM_STOKES:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end