function M  = mass_matrix(Base1,Base2,coord3,mesh,coeff)

%  build the maxx matrix 
%--------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_deg31,Dim1] = DOF_Mangement(mesh,Base1.DOF_3);
[Free_deg32,Dim2] = DOF_Mangement(mesh,Base2.DOF_3);
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);


M = sparse(Dim1,Dim2);
%-------------produce the mass matrix--------------------------------------
if Nr_elem > 0
    m = size(Base1.b_3([0,0]),2);
    n = size(Base2.b_3([0,0]),2);
    [area] = Transform_information(mesh,1,0);
    for kk=1:Nr_int3
        %----------Compute the coefficient matrix at the integral point----
        
        coeff_matrix = coeff((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));
                   
        B_1 = Base1.b_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
        B_2 = Base2.b_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2)); 
                   
        for i=1:m            
            for j=1:n                
              
               M       = M + coord3(kk,3)*sparse(Free_deg31(:,i),Free_deg32(:,j),...
                   B_1(i)*coeff_matrix*B_2(j).*area,Dim1,Dim2);
            end        
        end
    end     
end