function [ a, f ] = assemble_boundary ( node_num, node_xy, node_condition, ...
  ib, time, a, f )

%% ASSEMBLE_BOUNDARY modifies the linear system for the boundary conditions.
%
%  Discussion:
%
%    For now, we are only working with Dirichlet boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2007
%
%  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 IB, the half-bandwidth of the matrix.
%
%    Input, real TIME, the current time.
%
%    Input, real A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM coefficient 
%    matrix, stored in a compressed format.
%
%    Input, real F(NODE_NUM), the right hand side.
%
%    Output, real A(3*IB+1,NODE_NUM), the matrix has been adjusted for 
%    Dirichlet boundary conditions.
%
%    Output, real F(NODE_NUM), the right hand side has been adjusted for 
%    Dirichlet boundary conditions.
%
  DIRICHLET = 2;

  bc_value = dirichlet_condition ( node_num, node_xy, time );

  for node = 1 : node_num

    if ( node_condition(node) == DIRICHLET )

      column_low = max ( node - ib, 1 );
      column_high = min ( node + ib, node_num );

      for column = column_low : column_high
        a(node-column+2*ib+1,column) = 0.0;
      end
      a(2*ib+1,node) = 1.0;

      f(node) = bc_value(node);

    end

  end
 
  return
end