function a = jacobian_fem ( node_num, node_xy, element_num, ...
  element_node, quad_num, node_u_variable, node_v_variable, ...
  node_p_variable, variable_num, nu, c, ib )

%% JACOBIAN_FEM evaluates the Navier Stokes jacobian.
%
%  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 Navier Stokes equations in weak form are:
%
%      Integral ( nu * ( dBdx(I) * dUdx + dBdy(I) * dUdy ) 
%        + B(I) * ( ( U * dUdx + V * dUdy ) + dPdx - U_RHS ) ) = 0
%
%      Integral ( nu * ( dBdx(I) * dVdx + dBdy(I) * dVdy ) 
%        + B(I) * ( ( U * dVdx + V * dVdy ) + dPdy - V_RHS ) ) = 0
%
%      Integral ( Q(I) * ( dUdx + dVdy - P_RHS ) ) = 0
%
%    The residual terms due to boundary conditions must be applied separately.
%
%  Modified:
%
%    09 September 2006
%
%  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 ELEMENT_NUM, the number of element.
%
%    Input, integer ELEMENT_NODE(6,ELEMENT_NUM), the nodes of each element.
%
%    Input, integer QUAD_NUM, the order of the quadrature rule.
%
%    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_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 VARIABLE_NUM, the number of unknowns.
%
%    Input, real NU, the kinematic viscosity.
%
%    Input, real C(VARIABLE_NUM), the finite element coefficients of an 
%    approximate solution of the Navier Stokes equations.
%
%    Input, integer IB, the bandwidth of the jacobian.
%
%    Output, real A(3*IB+1,VARIABLE_NUM), the VARIABLE_NUM by 
%    VARIABLE_NUM Navier Stokes jacobian matrix, stored in a banded format.
%

%
%  Initialize the jacobian to 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 );
%
%  Consider all quantities associated with a given ELEMENT.
%
  for element = 1 : element_num
%
%  Make a copy of the triangle.
%
    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);
%
%  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 triangle.
%
    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));
%
%  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) * dbdy(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) * dbdy(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) * dqdy(1:3,1:quad_num);
%
%  dUeqn/dUcof,
%  dUeqn/dVcof,
%  dUeqn/dPcof.
%
    for i = 1 : 6
      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) ) ...
              + ...
              ( b(j,1:quad_num) .* dudx(1:quad_num) ...
              + u(1:quad_num) .* dbdx(j,1:quad_num) ...
              + v(1:quad_num) .* dbdy(j,1:quad_num) ) .* b(i,1:quad_num) ...
            ) ...
          );

        a(iu(i)-iv(j)+2*ib+1,iv(j)) = a(iu(i)-iv(j)+2*ib+1,iv(j)) + sum ...
        ( w(1:quad_num) .* b(j,1:quad_num) .* dudy(1:quad_num) ...
        .* b(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
%
%  dVeqn/dUcof,
%  dVeqn/dVcof,
%  dVeqn/dPcof.
%
    for i = 1 : 6
      for j = 1 : 6 

        a(iv(i)-iu(j)+2*ib+1,iu(j)) = a(iv(i)-iu(j)+2*ib+1,iu(j)) + sum ...
        ( w(1:quad_num) .* b(j,1:quad_num) .* dvdx(1:quad_num)            ...
          .* b(i,1:quad_num) );

        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) )            ...
            +                                                           ...
            ( u(1:quad_num) .* dbdx(j,1:quad_num)                        ...
            + b(j,1:quad_num) .* dvdy(1:quad_num)                        ...
            + v(1:quad_num) .* dbdy(j,1:quad_num) )                      ...
            .* b(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
%
%  dPeqn/dUcof,
%  dPeqn/dVcof,
%
    for i = 1 : 3
      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