function A = stiff_matrix(Base1,Base2,coord3,mesh,coeff)
%function A = stiff_matrix(base3_1,base3_2,base4_1,base4_2,coord3,coord4,mesh,coeff,DOF13,DOF14,DOF23,DOF24)

% Assemble the matrix corresponding the mesh, Base1,Base2 and the
% coifficient matrix coeff
%
% INPUT:  Base1, Base2, the left and right hand side base functions, coord3
% and coord4 is the quadrature matrix on the triangle and quadrilateral,
% mesh is the corresponding the triangulation, coeff is the coefficient
% matrix
%
% OUTPUT: A: is the matrix
%
% HEHU XIE,    02-10-2010



%--------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);
%----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);

A = sparse(Dim1,Dim2);
%-------------produce the stiffness matrix---------------------------------
if Nr_elem > 0
    m = size(Base1.gb_3([0,0]),2);
    n = size(Base2.gb_3([0,0]),2);
    [area,inv_transf] = Transform_information(mesh,1,1);
    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.gb_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
        B_2 = Base2.gb_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
                
               vv(:,1) = B_1(1,i)*inv_transf(:,1,1)+B_1(2,i)*inv_transf(:,2,1);
               vv(:,2) = B_1(1,i)*inv_transf(:,1,2)+B_1(2,i)*inv_transf(:,2,2);
        
               vv      = [vv(:,1).*coeff_matrix(:,1,1)+vv(:,2).*coeff_matrix(:,2,1),...
                          vv(:,1).*coeff_matrix(:,1,2)+vv(:,2).*coeff_matrix(:,2,2)];
               
               vv      = [vv(:,1).*inv_transf(:,1,1)+vv(:,2).*inv_transf(:,1,2),...
                          vv(:,1).*inv_transf(:,2,1)+vv(:,2).*inv_transf(:,2,2)];
               
               vv(:,1) = coord3(kk,3)*[B_2(1,j)*vv(:,1)+B_2(2,j)*vv(:,2)].*area; 
        
               A       = A + sparse(Free_deg31(:,i),Free_deg32(:,j),vv(:,1),Dim1,Dim2);
        
            end
        
        end
    end 
    
end