function node_r = residual_poisson ( node_num, node_xy, node_condition, ...
  element_num, element_node, quad_num, ib, a, f, node_u )

%% RESIDUAL_POISSON evaluates the residual for the Poisson equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 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_CONDITION(NODE_NUM), reports the condition
%    used to set the unknown associated with the node.
%    0, unknown.
%    1, finite element equation.
%    2, Dirichlet condition;
%    3, Neumann condition.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(3,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer QUAD_NUM, the number of quadrature points used in assembly.
%
%    Input, integer INDX(NODE_NUM), gives the index of the unknown quantity
%    associated with the given node.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Workspace, real A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM
%    coefficient matrix, stored in a compressed format.
%
%    Workspace, real F(NODE_NUM), the right hand side.
%
%    Input, real NODE_U(NODE_NUM), the value of the solution
%    at each node.
%
%    Output, real NODE_R(NODE_NUM), the finite element
%    residual at each node.
%
%  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.
%

%
%  Initialize the arrays to zero.
%
  f(1:node_num) = 0.0;
  a(1:3*ib+1,1:node_num) = 0.0;
%
%  Get the quadrature weights and nodes.
%
  [ quad_w, quad_xy ] = quad_rule ( quad_num );
%
%  The actual values of A and F are determined by summing up
%  contributions from all the elements.
%
  for element = 1 : element_num
%
%  Make a copy of the element.
%
    t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element));
%
%  Map the quadrature points QUAD_XY to points XY in the physical element.
%
    xy(1:2,1:quad_num) = reference_to_physical_t3 ( t3, quad_num, quad_xy );
    area = abs ( triangle_area_2d ( t3 ) );
    w(1:quad_num) = area * quad_w(1:quad_num);

    quad_rhs = rhs ( quad_num, quad_xy );

    quad_k = k_coef ( quad_num, quad_xy );
%
%  Consider a quadrature point QUAD, with coordinates (X,Y).
%
    for quad = 1 : quad_num
%
%  Consider one of the basis functions, which will play the
%  role of test function in the integral.
%
%  We generate an integral for every node associated with an unknown.
%  But if a node is associated with a boundary condition, we do nothing.
%
      for test = 1 : 3

        i = element_node(test,element);

        [ bi, dbidx, dbidy ] = basis_11_t3 ( t3, test, xy(1:2,quad) );

        f(i) = f(i) + w(quad) * quad_rhs(quad) * bi;
%
%  Consider another basis function, which is used to form the
%  value of the solution function.
%
        for basis = 1 : 3

          j = element_node(basis,element);

          [ bj, dbjdx, dbjdy ] = basis_11_t3 ( t3, basis, xy(1:2,quad) );

          a(i-j+2*ib+1,j) = a(i-j+2*ib+1,j) + w(quad) * ( ...
            dbidx * dbjdx + dbidy * dbjdy + quad_k(quad) * bj * bi );

        end

      end

    end

  end
%
%  Apply boundary conditions.
%
  [ a, f ] = dirichlet_apply ( node_num, node_xy, node_condition, ib, a, f );
%
%  Compute A*U.
%
  node_r = dgb_mxv ( node_num, node_num, ib, ib, a, node_u );
%
%  Set RES = A * U - F.
%
  node_r(1:node_num) = node_r(1:node_num) - f(1:node_num);

  return
end