function [A,RHS,BDDOF] = Poission_BoundaryTreatment(mesh,Base,u_D,A,RHS)

%-------------------------------------------------------------------------%
% Treat the boundary condition for Poission equation                      %
%-------------------------------------------------------------------------%

%-----get the local coordination in the reference edge [0 1)---------------
coord = [Base.nodal_3(3*[0:Base.DOF_3(1)-1]+1,1);...
         Base.nodal_3(3*(Base.DOF_3(1)+[0:Base.DOF_3(2)-1])+3,1)];
colu_RHS = size(RHS,2);
N_bd = size(coord,1);
%----deal with the boundary condition--------------------------------------
if Base.DOF_3(1) == 1
    %---we need to treat the boundary nodals    
    bd = unique(mesh.Dirichlet_node);
    dim1 = size(bd,1);
    A(bd,:) = 0;
    A(bd,bd) = eye(dim1);
    %RHS(bd,:) = 0;
    RHS(bd,:) = u_D(mesh.node(bd,1:2))*ones(1,colu_RHS);    
end
if Base.DOF_3(2) > 0
    Nr_edge = size(mesh.node4edge,1);
    bd_edge = unique(mesh.Dirichlet_edge);  
    dim2 = size(bd_edge,1);
end

if nargout == 3
    if Base.DOF_3(2) == 0
        BDDOF = bd;
    else
        BDDOF = zeros(dim1*Base.DOF_3(1)+dim2*Base.DOF_3(2),1);
        BDDOF(1:dim1,1) = bd;
        Pos = dim1;
    end
end
%---deal with the boundary condition on the Dirichlet edge-----------------
if  Base.DOF_3(2)>0
    Nr_edge = size(mesh.node4edge,1);
    bd_edge = unique(mesh.Dirichlet_edge);    
    %----start the Free_deg for boundary edge
    if Base.DOF_3(1)>0
        Nr_node = size(mesh.node,1);
        bd = Nr_node*Base.DOF_3(1) - Nr_edge + bd_edge;
    else
        bd = bd_edge - Nr_edge;
    end
    dim2 = size(bd,1);
end
%----treat boundary edges--------------------------------------------------
for i = Base.DOF_3(1) + 1 : N_bd
    bd = Nr_edge + bd;
    A(bd,:) = 0;
    A(bd,bd) = eye(dim2);
    RHS(bd,:) = u_D(mesh.node(mesh.node4edge(bd_edge,1),1:2)*(1-coord(i))+...
                mesh.node(mesh.node4edge(bd_edge,2),1:2)*coord(i))*...
                ones(1,colu_RHS);
     if nargout == 3
         BDDOF(Pos+1:Pos+dim2,1) = bd;
         Pos = Pos + dim2;
     end
end