function 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 )

%% NEUMANN_APPLY_SPARSE accounts for Neumann boundary conditions.
%
%  Discussion:
%
%    At the moment, this program only allows Neumann boundary conditions
%    of the form
%
%      dU/dn = 0
%      dV/dn = 0
%      dP/dn = 0
%
%    For such conditions, there is NO change necessary to the linear system.
%    So this routine actually does nothing.  It is here as preparation
%    for later treatment of nonzero Neumann conditions.
%
%  Modified:
%
%    08 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Input, 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.
%
%    Input, integer NODE_U_VARIABLE(NODE_NUM),
%    is the index of the horizontal velocity variable associated with the node.
%
%    Input, integer NODE_V_VARIABLE(NODE_NUM),
%    is the index of the vertical velocity variable associated with the node.
%
%    Input, 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.
%
%    Input, 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.
%
%    Input, 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.
%
%    Input, integer VARIABLE_NUM, the number of variables.
%
%    Input, real F(VARIABLE_NUM), the right hand side.
%
%    Output, real F(VARIABLE_NUM), the right hand side has been adjusted
%    for Dirichlet boundary conditions.
%

%
%  The user routine supplies a right hand side value for a possible
%  Neumann condition at EVERY node.
%
  [ u_bc, v_bc, p_bc ] = neumann_condition ( node_num, node_xy );

  NEUMANN = 3;

  for node = 1 : node_num

    iu = node_u_variable(node);
    iv = node_v_variable(node);
    ip = node_p_variable(node);

    if ( node_u_condition(node) == NEUMANN )

%     f(iu) = f(iu) + line integral;

    end

    if ( node_v_condition(node) == NEUMANN )

%     f(iv) = f(iv) + line integral;

    end

    if ( 0 < ip )

      if ( node_p_condition(node) == NEUMANN )

%       f(ip) = f(ip) + line integral;

      end

    end

  end