function [ a, f ] = dirichlet_apply ( node_num, node_xy, node_p_variable, ...
  node_u_variable, node_v_variable, node_p_condition, ...
  node_u_condition, node_v_condition, variable_num, ib, a, f )

%% DIRICHLET_APPLY accounts for Dirichlet boundary conditions.
%
%  Modified:
%
%    21 June 2005
%
%  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_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 NODE_P_CONDITION(NODE_NUM),
%    indicates the condition used to determine pressure at a node.
%    0, there is no condition at this node.
%    1, a finite element equation is used;
%    2, a Dirichlet condition is used.
%    3, a Neumann condition is used.
%
%    Input, integer NODE_U_CONDITION(NODE_NUM),
%    indicates the condition used to determine horizontal velocity at a node.
%    0, there is no condition at this node.
%    1, a finite element equation is used;
%    2, a Dirichlet condition is used.
%    3, a Neumann condition is used.
%
%    Input, integer NODE_V_CONDITION(NODE_NUM),
%    indicates the condition used to determine vertical velocity at a node.
%    0, there is no condition at this node.
%    1, a finite element equation is used;
%    2, a Dirichlet condition is used.
%    3, a Neumann condition is used.
%
%    Input, integer VARIABLE_NUM, the number of variables.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Input, real A(3*IB+1,VARIABLE_NUM), the VARIABLE_NUM by
%    VARIABLE_NUM coefficient matrix, stored in a banded format.
%
%    Input, real F(VARIABLE_NUM), the right hand side.
%
%    Output, real A(3*IB+1,VARIABLE_NUM), the matrix has been adjusted for 
%    Dirichlet boundary conditions.
%
%    Output, real F(VARIABLE_NUM), the right hand side has been adjusted
%    for Dirichlet boundary conditions.
%
  DIRICHLET = 2;

  for node = 1 : node_num

    [ u_bc, v_bc, p_bc ] = dirichlet_condition ( 1, node_xy(1:2,node) );

    iu = node_u_variable(node);
    iv = node_v_variable(node);
    ip = node_p_variable(node);

    if ( node_u_condition(node) == DIRICHLET )

      column_low = max ( iu - ib, 1 );
      column_high = min ( iu + ib, variable_num );

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

      f(iu) = u_bc;

    end

    if ( node_v_condition(node) == DIRICHLET )

      column_low = max ( iv - ib, 1 );
      column_high = min ( iv + ib, variable_num );

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

      f(iv) = v_bc;

    end

    if ( 0 < ip )

      if ( node_p_condition(node) == DIRICHLET )

        column_low = max ( ip - ib, 1 );
        column_high = min ( ip + ib, variable_num );

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

        f(ip) = p_bc;

      end

    end

  end