function [A,RHS] = Navier_Stokes_BoundaryTreatment(mesh,Base,u1_D,u2_D,A,RHS)

% Treat the boundary condition for Navier_Stokes equation
%  Dim: the dimension for the rows of the velocity fem space


%-----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)];
%------get the information of the mesh-------------------------------------
Nr_node=size(mesh.node,1);
Nr_edge=size(mesh.node4edge,1);
Nr_elem=size(mesh.elem,1);
%----get the number of triangles-------------------------------------------
II3 = find(mesh.elem(:,4)<0);
Nr_elem3 = length(II3);
II4 = setdiff([1:Nr_elem]',II3);
Nr_elem4 = length(II4);
Dim = Nr_node*Base.DOF_3(1)+Nr_edge*Base.DOF_3(2)+Nr_elem3*Base.DOF_3(3)+Nr_elem4*Base.DOF_4(3);

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);
    dim2=size(bd,1);
    A(bd,:)=0;
    A(Dim+bd,:)=0;
    A(bd,bd)=eye(dim2);
    A(Dim+bd,Dim+bd)=eye(dim2);
    RHS(bd,:)=u1_D(mesh.node(bd,1:2))*ones(1,colu_RHS);
    RHS(Dim+bd,:)=u2_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);
    %----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(Dim+bd,:)=0;
    A(bd,bd)=eye(dim2);
    A(Dim+bd,Dim+bd) = eye(dim2);
    RHS(bd,:)=u1_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);
    RHS(Dim+bd,:)=u2_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);
end
% deal with the integral 0 for the pressure variable 
A(end,:) = 0;
A(end,end) = 1;
RHS(end,:) = 0;