function err = difference_estimate(mesh,Base1,solu1,Base0,solu0,coord3,type)

% err=difference_estimate(mesh,Base1,solu1,Base0,solu0,coord3,coord4,type)
% Compute the square of the error between the approximation solution1 and
% another finite element solution solution0 
% INPUT: mesh, the finite element mesh,
%        Base1 is the base for the finite element solution solution1,
%        solution1 is the first solution
%        Base0 is the base for the another finite element solution solution0,
%        coord3 and coord4 are the quadrature informations for triangle and
%        quadrilaterals
%       type: means compute the H_1 error or mass error
%       1: means compute the H_1 error 
%       0:  means compute the L_2 error
% OUTPUT: error is the desired 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_deg1_3,Dim1] = DOF_Mangement(mesh,Base1.DOF_3);
[Free_deg0_3,Dim0] = DOF_Mangement(mesh,Base0.DOF_3);
%----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);
          
err=0;
if type == 1
    %------compute the error for the gradient approximation----------------        
    if Nr_elem > 0        
        m1 = size(Base1.gb_3([0,0]),2);
        m0 = size(Base0.gb_3([0,0]),2);
        [area,inv_transf] = Transform_information(mesh,1,1);
        for kk = 1:Nr_int3            
           B_1 = Base1.gb_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));
           B_0 = Base0.gb_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));
           %---------compute the gradient of the approximation-------------         
            v = zeros(Nr_elem,2);               
            for i=1:m1                
                v = v + [(inv_transf(:,1,1)*B_1(1,i)+inv_transf(:,2,1)*B_1(2,i)).*solu1(Free_deg1_3(:,i)),...
                         (inv_transf(:,1,2)*B_1(1,i)+inv_transf(:,2,2)*B_1(2,i)).*solu1(Free_deg1_3(:,i))];                    
            end            
            for i=1:m0                
                v = v - [(inv_transf(:,1,1)*B_0(1,i)+inv_transf(:,2,1)*B_0(2,i)).*solu0(Free_deg0_3(:,i)),...
                         (inv_transf(:,1,2)*B_0(1,i)+inv_transf(:,2,2)*B_0(2,i)).*solu0(Free_deg0_3(:,i))];                 
            end            
            err = err + coord3(kk,3)*sum(sum(v.^2,2).*area);
        end        
    end
      
   %return    
else
    err=0;
    %----------compute the error for function value approximation----------
    if Nr_elem > 0        
        [area]=Transform_information(mesh,1,0);        
        for kk = 1: Nr_int3            
            B_1 = Base1.b_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));
            B_0 = Base0.b_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));            
            %----------compute the error-----------------------------------                   
            err = err + coord3(kk,3)*sum((solu1(Free_deg1_3)*B_1' - solu0(Free_deg0_3)*B_0').^2.*area);
        end        
    end
end
err =sqrt(err);