function Ke = computeKE(a,b,DH)
    GaussNodes = [-1/sqrt(3); 1/sqrt(3)]; GaussWeigh = [1 1]; 
    L = [1 0 0 0; 0 0 0 1; 0 1/2 1/2 0]; 
    Ke = zeros(8,8); 

    for i = 1:length(GaussNodes) 
        for j = 1:length(GaussNodes)      
            GN_x = a*GaussNodes(i); 
            GN_y = b*GaussNodes(j);            
            
            dN_x = 1/4*[-(1-GN_x)  (1-GN_x) (1+GN_x) -(1+GN_x)];            
            dN_y = 1/4*[-(1-GN_y) -(1+GN_y) (1+GN_y)  (1-GN_y)];   

            J = [dN_x; dN_y]*[ -a  a  a  -a;  
                                -b  -b  b  b]';            
            G = [inv(J) zeros(size(J)); 
                zeros(size(J)) inv(J)];      

            dN(1,1:2:8) = dN_x; dN(2,1:2:8) = dN_y;            
            dN(3,2:2:8) = dN_x; dN(4,2:2:8) = dN_y;   

            Be = L*G*dN;        
            Ke = Ke + GaussWeigh(i)*GaussWeigh(j)*det(J)*Be'*DH*Be; 
        end  
    end  
    
% D = [d11 d12 d13
%      d12 d22 d23
%      d13 d23 d33];
% d11 = DH(1,1);
% d12 = DH(1,2);
% d13 = DH(1,3);
% d22 = DH(2,2);
% d23 = DH(2,3);
% d33 = DH(3,3);
% Ke = ...
%  [               0.3333*d11 + 0.25*d13 + 0.08333*d33,     0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33,                          0.04167*d33 - 0.3333*d11,  0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33,               - 0.1667*d11 - 0.25*d13 - 0.04167*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33,                          0.1667*d11 - 0.08333*d33,  0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33;
%      0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33,                 0.3333*d22 + 0.25*d23 + 0.08333*d33,  0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33,                          0.1667*d22 - 0.08333*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33,               - 0.1667*d22 - 0.25*d23 - 0.04167*d33,  0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33,                          0.04167*d33 - 0.3333*d22;
%                             0.04167*d33 - 0.3333*d11,    0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33,               0.3333*d11 - 0.25*d13 + 0.08333*d33,   0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33,                            0.1667*d11 - 0.08333*d33,    0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33,               0.25*d13 - 0.1667*d11 - 0.04167*d33, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33;
%     0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33,                            0.1667*d22 - 0.08333*d33,   0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33,               0.3333*d22 - 0.25*d23 + 0.08333*d33,    0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33,                            0.04167*d33 - 0.3333*d22, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33,               0.25*d23 - 0.1667*d22 - 0.04167*d33;
%                - 0.1667*d11 - 0.25*d13 - 0.04167*d33, - 0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33,                          0.1667*d11 - 0.08333*d33,  0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33,                 0.3333*d11 + 0.25*d13 + 0.08333*d33,     0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33,                          0.04167*d33 - 0.3333*d11,  0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33;
%   -0.25*d12 - 0.08333*d13 - 0.08333*d23 - 0.0625*d33,               - 0.1667*d22 - 0.25*d23 - 0.04167*d33,  0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33,                          0.04167*d33 - 0.3333*d22,     0.25*d12 + 0.1667*d13 + 0.1667*d23 + 0.0625*d33,                 0.3333*d22 + 0.25*d23 + 0.08333*d33,  0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33,                          0.1667*d22 - 0.08333*d33;
%                             0.1667*d11 - 0.08333*d33,    0.25*d12 + 0.08333*d13 - 0.1667*d23 - 0.0625*d33,               0.25*d13 - 0.1667*d11 - 0.04167*d33, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33,                            0.04167*d33 - 0.3333*d11,    0.08333*d23 - 0.1667*d13 - 0.25*d12 + 0.0625*d33,               0.3333*d11 - 0.25*d13 + 0.08333*d33,   0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33;
%     0.08333*d13 - 0.25*d12 - 0.1667*d23 + 0.0625*d33,                            0.04167*d33 - 0.3333*d22, 0.25*d12 - 0.08333*d13 - 0.08333*d23 + 0.0625*d33,               0.25*d23 - 0.1667*d22 - 0.04167*d33,    0.25*d12 - 0.1667*d13 + 0.08333*d23 - 0.0625*d33,                            0.1667*d22 - 0.08333*d33,   0.1667*d13 - 0.25*d12 + 0.1667*d23 - 0.0625*d33,               0.3333*d22 - 0.25*d23 + 0.08333*d33;];

end