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

%% ASSEMBLE_STOKES assembles the finite element system for steady 2D Stokes flow.
%
%  Discussion:
%
%    The matrix is known to be banded.  A special matrix storage format
%    is used to reduce the space required.  Details of this format are
%    discussed in the routine DGB_FA.
%
%    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:
%
%    16 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 IB, the matrix half-bandwidth.
%
%    Output, real A(3*IB+1,VARIABLE_NUM), the VARIABLE_NUM by VARIABLE_NUM
%    finite element coefficient matrix.
%
%    Output, real F(VARIABLE_NUM), the finite element right hand side.
%
  f(1:variable_num) = 0.0;
  a(1:3*ib+1,1:variable_num) = 0.0;
%
%  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(1:2,1:quad_num) = 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 );
%
%  Evaluate the basis functions at the quadrature points.
%
    [ b, dbdx, dbdy ] = basis_mn_t6 ( t6, quad_num, xy );
    [ q, dqdx, dqdy ] = basis_mn_t3 ( t3, quad_num, xy );
%
%  Extract the indices of the finite element coefficients for this element.
%
    iu(1:6) = node_u_variable(element_node(1:6,element));
    iv(1:6) = node_v_variable(element_node(1:6,element));
    ip(1:3) = node_p_variable(element_node(1:3,element));
%
%  The horizontal momentum equation.
%
    for i = 1 : 6

      f(iu(i)) = f(iu(i)) + sum ...
      ( w(1:quad_num) .* u_rhs(1:quad_num) .* b(i,1:quad_num) );

      for j = 1 : 6

        a(iu(i)-iu(j)+2*ib+1,iu(j)) = a(iu(i)-iu(j)+2*ib+1,iu(j)) + sum ...
        ( w(1:quad_num) .*                                               ...
          (                                                             ...
            nu * ( dbdx(j,1:quad_num) .* dbdx(i,1:quad_num)              ...
                 + dbdy(j,1:quad_num) .* dbdy(i,1:quad_num) )            ...
          )                                                             ...
        );

      end

      for j = 1 : 3
        a(iu(i)-ip(j)+2*ib+1,ip(j)) = a(iu(i)-ip(j)+2*ib+1,ip(j)) + sum ...
          ( w(1:quad_num) .* dqdx(j,1:quad_num) .* b(i,1:quad_num) );
      end

    end
%
%  The vertical momentum equation.
%
    for i = 1 : 6

      f(iv(i)) = f(iv(i)) + sum ...
      ( w(1:quad_num) .* b(i,1:quad_num) .* v_rhs(1:quad_num) );

      for j = 1 : 6

        a(iv(i)-iv(j)+2*ib+1,iv(j)) = a(iv(i)-iv(j)+2*ib+1,iv(j)) + sum ...
        ( w(1:quad_num) .*                                               ...
          (                                                             ...
            nu * ( dbdx(j,1:quad_num) .* dbdx(i,1:quad_num)              ...
                 + dbdy(j,1:quad_num) .* dbdy(i,1:quad_num) )            ...
          )                                                             ...
        );

      end

      for j = 1 : 3
        a(iv(i)-ip(j)+2*ib+1,ip(j)) = a(iv(i)-ip(j)+2*ib+1,ip(j)) + sum ...
        ( w(1:quad_num) .* dqdy(j,1:quad_num) .* b(i,1:quad_num) );
      end

    end
%
%  The pressure equation.
%
    for i = 1 : 3

      f(ip(i)) = f(ip(i)) + sum ...
      ( w(1:quad_num) .* q(i,1:quad_num) .* p_rhs(1:quad_num) );

      for j = 1 : 6

        a(ip(i)-iu(j)+2*ib+1,iu(j)) = a(ip(i)-iu(j)+2*ib+1,iu(j)) + sum ...
        ( w(1:quad_num) .* dbdx(j,1:quad_num) .* q(i,1:quad_num) );

        a(ip(i)-iv(j)+2*ib+1,iv(j)) = a(ip(i)-iv(j)+2*ib+1,iv(j)) + sum ...
        ( w(1:quad_num) .* dbdy(j,1:quad_num) .* q(i,1:quad_num) );

      end

    end

  end