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

%% ASSEMBLE_STOKES_SPARSE assembles the Stokes stiffness SPARSE matrix.
%
%  Discussion:
%
%    The finite element coefficient matrix is set up as a MATLAB
%    sparse matrix.
%
%    The Stokes equations in weak form are:
%
%      Integral ( nu * ( dBdx(I) * dUdx + dBdy(I) * dUdy ) 
%        + B(I) * ( dPdx - U_RHS ) ) = 0
%
%      Integral ( nu * ( dBdx(I) * dVdx + dBdy(I) * dVdy ) 
%        + B(I) * (  dPdy - V_RHS ) ) = 0
%
%      Integral ( Q(I) * ( dUdx + dVdy - P_RHS ) ) = 0
%
%    Once the basic finite element system is set up by this routine, another 
%    routine adjusts the system to account for boundary conditions.
%
%  Modified:
%
%    25 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer QUAD_NUM, the number of quadrature points in an element.
%
%    Input, integer VARIABLE_NUM, the number of unknowns.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the nodes.
%
%    Input, integer NODE_P_VARIABLE(NODE_NUM), the index of the pressure 
%    variable associated with a node, or -1 if there is none.
%
%    Input, integer NODE_U_VARIABLE(NODE_NUM), the index of the horizontal 
%    velocity variable associated with a node, or -1 if there is none.
%
%    Input, integer NODE_V_VARIABLE(NODE_NUM), the index of the vertical 
%    velocity variable associated with a node, or -1 if there is none.
%
%    Input, 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.
%
%    Input, real NU, the kinematic viscosity.
%
%    Input, integer NZ_NUM, the (maximum) number of nonzeros in the matrix.
%    If set to 0 on input, we hope MATLAB's sparse utility will be able
%    to take over the task of reallocating space as necessary.
%
%    Output, real sparse A, the finite element coefficient matrix, stored
%    as a MATLAB sparse matrix.
%
%    Output, real F(VARIABLE_NUM), the finite element right hand side.
%
%  Local parameters:
%
%    Local, real BI, DBIDX, DBIDY, the value of some basis function
%    and its first derivatives at a quadrature point.
%
%    Local, real BJ, DBJDX, DBJDY, the value of another basis 
%    function and its first derivatives at a quadrature point.
%
  f(1:variable_num) = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASSEMBLE_STOKES_SPARSE:\n' );
  fprintf ( 1, '  Setting up sparse Stokes matrix with NZ_NUM = %d\n', nz_num );

  a = sparse ( [], [], [], variable_num, variable_num, nz_num );
%
%  Get the quadrature weights and nodes.
%
  [ quad_w, quad_xy ] = quad_rule ( quad_num );

  for element = 1 : element_num
%
%  Extract the nodes of the linear and quadratic triangles.
%
    t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element));
    t6(1:2,1:6) = node_xy(1:2,element_node(1:6,element));
%
%  Map the quadrature points QUAD_XY to points XY in the physical triangle.
%
    xy = reference_to_physical_t6 ( t6, quad_num, quad_xy );
    area = abs ( triangle_area_2d ( t3 ) );
    w(1:quad_num) = area * quad_w(1:quad_num);
    [ u_rhs, v_rhs, p_rhs ] = rhs ( quad_num, xy );
%
%  Consider the QUAD-th quadrature point.
%
    for quad = 1 : quad_num

      point = xy(1:2,quad);
%
%  Consider the test functions.
%
      [ bi, dbidx, dbidy ] = basis_mn_t6 ( t6, 1, point );
      [ qi, dqidx, dqidy ] = basis_mn_t3 ( t3, 1, point );

      for test = 1 : 6

        test_node = element_node(test,element);

        iu = node_u_variable(test_node);
        iv = node_v_variable(test_node);
        ip = node_p_variable(test_node);
%
%  Compute the source terms for the right hand side.
%
        f(iu) = f(iu) + w(quad) * u_rhs(quad) * bi(test);
        f(iv) = f(iv) + w(quad) * v_rhs(quad) * bi(test);
        if ( 0 < ip )
          f(ip) = f(ip) + w(quad) * p_rhs(quad) * qi(test);
        end
%
%  Consider the basis functions.
%
        [ bj, dbjdx, dbjdy ] = basis_mn_t6 ( t6, 1, point );
        [ qj, dqjdx, dqjdy ] = basis_mn_t3 ( t3, 1, point );

        for basis = 1 : 6

          basis_node = element_node(basis,element);

          ju = node_u_variable(basis_node);
          jv = node_v_variable(basis_node);
          jp = node_p_variable(basis_node);
%
%  Add terms to the horizonal momentum equation.
%
          a(iu,ju) = a(iu,ju) + w(quad) * nu ...
            * ( dbidx(test) * dbjdx(basis) + dbidy(test) * dbjdy(basis) );

          if ( 0 < jp ) 
            a(iu,jp) = a(iu,jp) + w(quad) * bi(test) * dqjdx(basis);
          end
%
%  Add terms to the vertical momentum equation.
%
          a(iv,jv) = a(iv,jv) + w(quad) * nu ...
            * ( dbidx(test) * dbjdx(basis) + dbidy(test) * dbjdy(basis) );

          if ( 0 < jp )
            a(iv,jp) = a(iv,jp) + w(quad) * bi(test) * dqjdy(basis);
          end
%
%  Add terms to the continuity equation.
%
          if ( 0 < ip )
            a(ip,ju) = a(ip,ju) + w(quad) * qi(test) * dbjdx(basis);
            a(ip,jv) = a(ip,jv) + w(quad) * qi(test) * dbjdy(basis);
          end

        end

      end

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASSEMBLE_STOKES_SPARSE:\n' );
  fprintf ( 1, '  Sparse Stokes matrix used NZ_NUM = %d\n', nz_num );