function [mesh,solu,Prolong] = refinement_multi(mesh,solu,Base)

%---------refine the mesh by the regular way-------------------------------
tic;

fprintf('----------Refine the mesh by regular way---------------------\n')

%--------------------------------------------------------------------------
Nr_edge=size(mesh.node4edge,1);
Nr_node=size(mesh.node,1);
Nr_elem=size(mesh.elem,1);


numnew = 3*Nr_elem;
if (nargin == 3 & nargout >= 2) 
    [FreeDeg31,Dim1]=DOF_Mangement(mesh,Base.DOF_3);
    dim3_basis = size(Base.b_3([0,0]));
end

%---Chech if to do the multilevel correction 
Ind_Multilevel = isfield(mesh,'owner');
if Ind_Multilevel
    owner = mesh.owner;
    elemlocalnode = mesh.elemlocalnode;
end

%--------------get the new points(mid points)------------------------------
mesh.node=[mesh.node;
           [(mesh.node(mesh.node4edge(:,1),1:2)+...
            mesh.node(mesh.node4edge(:,2),1:2))/2]];    
%------------renew the mesh construct--------------------------------------
mesh.elem =  [mesh.elem(:,1),Nr_node+mesh.edge(:,[3,2]);
             Nr_node+mesh.edge(:,3),mesh.elem(:,2),Nr_node+mesh.edge(:,1);
             Nr_node+mesh.edge(:,2),Nr_node+mesh.edge(:,1),mesh.elem(:,3);
             Nr_node+mesh.edge(:,1:3)];

%------------get the two level relationship--------------------------------
% Rel: save the coarse mesh inclue which elements in the fine mesh
if nargout >= 2 & nargin >= 2 
    Rel(1:Nr_elem,:) = [1:Nr_elem]'*ones(1,4)+Nr_elem*ones(Nr_elem,1)*[0 1 2 3]; 
end                        
%--------------------------------------------------------------------------
if size(mesh.Dirichlet_node,1)>0
    mesh.Dirichlet_node = [mesh.Dirichlet_node(:,1),Nr_node+mesh.Dirichlet_edge;
                       Nr_node+mesh.Dirichlet_edge, mesh.Dirichlet_node(:,2)];
end
%--------------------------------------------------------------------------
if size(mesh.Neumann_node,1)>0              
    mesh.Neumann_node = [mesh.Neumann_node(:,1),Nr_node+mesh.Neumann_edge;
                         Nr_node+mesh.Neumann_edge,mesh.Neumann_node(:,2)];
end
                
%--------------------------------------------------------------------------
%------produce the new mesh------------------------------------------------
mesh = Initial_FEM_matlab(mesh.node,mesh.elem,mesh.Dirichlet_node,mesh.Neumann_node);
if Ind_Multilevel 
    %----------renew the owner information for the new mesh----------------
    %----------owner: save the owner information to save the original element
    %--comes from the coarse elements--------------------------------------
    mesh.owner = [owner; owner; owner; owner];
    %----------------------------------------------------------------------
    mesh.elemlocalnode = sparse(4*Nr_elem,9);                                                
    mesh.elemlocalnode(1:Nr_elem,:) = [elemlocalnode(:,1:3),...
                                (elemlocalnode(:,1:3)+elemlocalnode(:,4:6))/2,...
                                (elemlocalnode(:,1:3)+elemlocalnode(:,7:9))/2];
    mesh.elemlocalnode(Nr_elem+1:2*Nr_elem,:) = [(elemlocalnode(:,1:3)+elemlocalnode(:,4:6))/2,...
                                 elemlocalnode(:,4:6),...
                                (elemlocalnode(:,4:6)+elemlocalnode(:,7:9))/2];
    mesh.elemlocalnode(2*Nr_elem+1:3*Nr_elem,:) = [(elemlocalnode(:,1:3)+elemlocalnode(:,7:9))/2,...
                                (elemlocalnode(:,4:6)+elemlocalnode(:,7:9))/2,...
                                 elemlocalnode(:,7:9),];
    mesh.elemlocalnode(3*Nr_elem+1:4*Nr_elem,:) = [(elemlocalnode(:,4:6)+elemlocalnode(:,7:9))/2,...
                                (elemlocalnode(:,1:3)+elemlocalnode(:,7:9))/2,...
                                (elemlocalnode(:,1:3)+elemlocalnode(:,4:6))/2];
end
%-------renew the solution-------------------------------------------------
%  !!!! So far only for the linear function--------------------------------
if (nargin == 3  & nargout >= 2)     
    [FreeDeg32,Dim2]=DOF_Mangement(mesh,Base.DOF_3);
    dim_solu = size(solu,2);
    solu(Dim1+1:Dim2,:) = sparse(Dim2-Dim1,dim_solu);
    solu_tmp = zeros(Dim2,1);
    if(nargout ==3 )
        %----define the prolong matrix------
        Prolong = sparse(Dim2,Dim1);
    end
    if Nr_elem > 0
        %--define the left basis and right hand basis--------------------
        mesh_left = struct('node',[0 0;1/2 0;0 1/2], 'elem',[1 2 3],...
            'edge',[1 2 3],'node4edge',[2 3; 3 1; 1 2]);
        left_basis = Base.nodal_funct_3(Base.b_3, mesh_left)';
        %----------------------------------------------------------------
        mesh_right = struct('node',[1/2 0; 1 0; 1/2 1/2], 'elem',[1 2 3],...
            'edge',[1 2 3],'node4edge',[2 3; 3 1; 1 2]);
        right_basis = Base.nodal_funct_3(Base.b_3, mesh_right)';
        %--define the left basis and right hand basis--------------------
        mesh_up = struct('node',[0 1/2;1/2 1/2;0 1], 'elem',[1 2 3],...
            'edge',[1 2 3],'node4edge',[2 3; 3 1; 1 2]);
        up_basis = Base.nodal_funct_3(Base.b_3, mesh_up)';
        %----------------------------------------------------------------
        mesh_mid = struct('node',[1/2 1/2;0 1/2;1/2 0], 'elem',[1 2 3],...
            'edge',[1 2 3],'node4edge',[2 3; 3 1; 1 2]);
        mid_basis = Base.nodal_funct_3(Base.b_3, mesh_mid)';
        %----------------------------------------------------------------
    end        
        %------------------------------------------------------------------        
        if(nargout==3)
            if Nr_elem >= 1               
                for i=1:Nr_elem
                    Prolong(FreeDeg32(Rel(i,1),:),FreeDeg31(i,:)) = left_basis';
                    Prolong(FreeDeg32(Rel(i,2),:),FreeDeg31(i,:)) = right_basis';
                    Prolong(FreeDeg32(Rel(i,3),:),FreeDeg31(i,:)) = up_basis';
                    Prolong(FreeDeg32(Rel(i,4),:),FreeDeg31(i,:)) = mid_basis';
                end
            end 
            %Prolong = Prolong';
        end
        %-------Prolong the solution to the new mesh-----------------------
        for k=1:dim_solu
            solu_init = solu(:,k);
            if Nr_elem > 1 
                solu_tmp(FreeDeg32(Rel(:,1),:)) = solu_init(FreeDeg31)*left_basis;
                solu_tmp(FreeDeg32(Rel(:,2),:)) = solu_init(FreeDeg31)*right_basis;
                solu_tmp(FreeDeg32(Rel(:,3),:)) = solu_init(FreeDeg31)*up_basis;
                solu_tmp(FreeDeg32(Rel(:,4),:)) = solu_init(FreeDeg31)*mid_basis;
            end

            if Nr_elem == 1
                solu_tmp(FreeDeg32(Rel(:,1),:)) = solu_init(FreeDeg31)'*left_basis;
                solu_tmp(FreeDeg32(Rel(:,2),:)) = solu_init(FreeDeg31)'*right_basis;
                solu_tmp(FreeDeg32(Rel(:,3),:)) = solu_init(FreeDeg31)'*up_basis;
                solu_tmp(FreeDeg32(Rel(:,4),:)) = solu_init(FreeDeg31)'*mid_basis;
            end
            solu(:,k) = solu_tmp;  
        end
end
toc