function r = residual_navier_stokes ( node_num, node_xy, triangle_num, ...
  triangle_node, quad_num, node_u_variable, node_v_variable, ...
  node_p_variable, variable_num, reynolds, c )

%% RESIDUAL_NAVIER_STOKES evaluates the Navier Stokes residual.
%
%  Discussion:
%
%    The residual terms due to boundary conditions must be applied separately.
%
%    The Navier Stokes equations in weak form are:
%
%      Integral ( ( dBdx(I) * dUdx + dBdy(I) * dUdy ) 
%        + B(I) * ( R * ( U * dUdx + V * dUdy ) + dPdx - U_RHS ) ) = 0
%
%      Integral ( ( dBdx(I) * dVdx + dBdy(I) * dVdy ) 
%        + B(I) * ( R * ( U * dVdx + V * dVdy ) + dPdy - V_RHS ) ) = 0
%
%      Integral ( Q(I) * ( dUdx + dVdy - P_RHS ) ) = 0
%
%  Modified:
%
%    19 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the nodes.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(6,TRIANGLE_NUM), the nodes of each triangle.
%
%    Input, integer QUAD_NUM, the order of the quadrature rule.
%
%    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 VARIABLE_NUM, the number of unknowns.
%
%    Input, real REYNOLDS, the value of the Reynolds number.
%
%    Input, real C(VARIABLE_NUM), the finite element coefficients of an 
%    approximate solution of the Navier Stokes equations.
%
%    Output, real R(VARIABLE_NUM), the Navier Stokes residual.
%

%
%  Initialize the residual to 0.
%
  r(1:variable_num) = 0.0;
%
%  Get the quadrature weights and nodes.
%
  [ quad_w, quad_xy ] = quad_rule ( quad_num );
%
%  Consider all quantities associated with a given TRIANGLE.
%
  for triangle = 1 : triangle_num
%
%  Make a copy of the triangle.
%
    t3(1:2,1:3) = node_xy(1:2,triangle_node(1:3,triangle));
    t6(1:2,1:6) = node_xy(1:2,triangle_node(1:6,triangle));
%
%  Map the quadrature points QUAD_XY to points XY in the physical triangle.
%
    xy(1:2,1:quad_num) = reference_to_physical_t3 ( t3, quad_num, quad_xy );

    w(1:quad_num) = quad_w(1:quad_num) * triangle_area ( t3 );
%
%  Evaluate the basis functions at the quadrature points.
%
    [ b, dbdx, dbdy ] = basis_t6 ( t6, quad_num, xy );
    [ q, dqdx, dqdy ] = basis_t3 ( t3, quad_num, xy );
%
%  Extract the indices of the finite element coefficients for this triangle.
%
    iu(1:6) = node_u_variable(triangle_node(1:6,triangle));
    iv(1:6) = node_v_variable(triangle_node(1:6,triangle));
    ip(1:3) = node_p_variable(triangle_node(1:3,triangle));
%
%  Extract the finite element coefficients for this triangle.
%
    cu(1:6) = c(iu(1:6));
    cv(1:6) = c(iv(1:6));
    cp(1:3) = c(ip(1:3));
%
%  Evaluate the flowfield at each quadrature point.
%
     u  (1:quad_num) = cu(1:6) *  b  (1:6,1:quad_num);
    dudx(1:quad_num) = cu(1:6) * dbdx(1:6,1:quad_num);
    dudy(1:quad_num) = cu(1:6) * dbdx(1:6,1:quad_num);

     v  (1:quad_num) = cv(1:6) *  b  (1:6,1:quad_num);
    dvdx(1:quad_num) = cv(1:6) * dbdx(1:6,1:quad_num);
    dvdy(1:quad_num) = cv(1:6) * dbdx(1:6,1:quad_num);

     p  (1:quad_num) = cp(1:3) *  q  (1:3,1:quad_num);
    dpdx(1:quad_num) = cp(1:3) * dqdx(1:3,1:quad_num);
    dpdy(1:quad_num) = cp(1:3) * dqdx(1:3,1:quad_num);
%
%  Evaluate the right hand side at each quadrature point.
%
    [ u_rhs, v_rhs, p_rhs ] = rhs ( quad_num, xy );
%
%  Consider the contribution to the integrals from the QUAD-th quadrature point.
%
    for quad = 1 : quad_num

      for i = 1 : 6

        r(iu(i)) = r(iu(i)) + w(quad) * ( ...
          ( dbdx(i,quad) * dudx(quad) + dbdy(i,quad) * dudy(quad) ) ...
          + b(i,quad) * ( ...
          reynolds * ( u(quad) * dudx(quad) + v(quad) * dudy(quad) ) ...
          + dpdx(quad) - u_rhs(quad) ) );

        r(iv(i)) = r(iv(i)) + w(quad) * ( ...
          ( dbdx(i,quad) * dvdx(quad) + dbdy(i,quad) * dvdy(quad) ) ...
          + b(i,quad) * ( ...
          reynolds * ( u(quad) * dvdx(quad) + v(quad) * dvdy(quad) ) ...
          + dpdy(quad) - v_rhs(quad) ) );

      end

      for i = 1 : 3

        r(ip(i)) = r(ip(i)) + w(quad) * q(i,quad) * ( ...
          dudx(quad) + dvdy(quad) - p_rhs(quad) );

      end

    end

  end
