function err = error_estimate(mesh,Base,coord3,u_h,u,type)

% err=error_estimate(mesh,Base,coord3,coord4,u_h,u,type)
% Compute the square of the error between the approximation u_h and the 
% exact funtion u mesh describe the mesh structure
% base is the bases function
% coord3 are the quadrature informations
% u_h is the finite element approximation
% u is the exact function
% type: means compute the H_1 error or mass error
%       gradient: means compute the H_1 error 
%       else:  means compute the L_2 error
%
%
%  Hehu Xie   23-05-2009




%--------get the element matrix size---------------------------------------
Nr_int3 = size(coord3,1);
%--------compute the base value at the integral points---------------------
Nr_elem=size(mesh.elem,1);
%----------compute the freedom degree vector-------------------------------
[Free_deg3,Dim] = DOF_Mangement(mesh,Base.DOF_3,Base.DOF_4);
%----get the number of triangles-------------------------------------------

          
err=0;
if type=='gradient'
    %------compute the error for the gradient approximation----------------
       
    if Nr_elem > 0
        m=size(Base.gb_3([0,0]),2);

        [area,inv_transf]=Transform_information(mesh,1,1);
        for kk=1:Nr_int3
            
           B = Base.gb_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
           
           grad_u = u((1-coord3(kk,1)-coord3(kk,2))*mesh.node(mesh.elem(:,1),1:2)+...
                       coord3(kk,1)*mesh.node(mesh.elem(:,2),1:2)+...
                       coord3(kk,2)*mesh.node(mesh.elem(:,3),1:2));
    
            %---------compute the gradient of the approximation------------         
            v=zeros(Nr_elem,2);
    
            for i=1:m
                
                v = v + [(inv_transf(:,1,1)*B(1,i)+inv_transf(:,2,1)*B(2,i)).*u_h(Free_deg3(:,i)),...
                            (inv_transf(:,1,2)*B(1,i)+inv_transf(:,2,2)*B(2,i)).*u_h(Free_deg3(:,i))];
            end
            err = err + coord3(kk,3)*sum(sum((grad_u-v).^2,2).*area);
        end        
    end
    
else
   
    %----------compute the error for function value approximation----------
    if Nr_elem > 0
        
        [area]=Transform_information(mesh,1,0);
        
        for kk = 1: Nr_int3            
            B=Base.b_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
            
            %----------compute the error-----------------------------------
            
            u_exact = u((1-coord3(kk,1)-coord3(kk,2))*mesh.node(mesh.elem(:,1),1:2)+...
                                coord3(kk,1)*mesh.node(mesh.elem(:,2),1:2)+...
                                coord3(kk,2)*mesh.node(mesh.elem(:,3),1:2));
            err = err + coord3(kk,3)*sum((u_h(Free_deg3)*B' - u_exact).^2.*area);
        end        
    end
end
err =sqrt(err);