%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-10-06 14:49:15
 % @LastEditTime: 2021-10-06 14:51:23
%--------------------------
function [den, energy_nod, comp] = simulation_MMC(edofMat, EleNodesID, iEner, H, iK, jK, KE, freedofs, U, F)
    % FEA 
    denk= sum(H(EleNodesID).^2, 2) / 4; 
    den = sum(H(EleNodesID), 2) / 4; 

    sK= KE(:)*denk(:)'; 
    K = sparse(iK(:),jK(:),sK(:)); K=(K+K')/2; 
    U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); 
    
    %Energy of element 
    energy = sum((U(edofMat)*KE).*U(edofMat), 2); 
    sEner = ones(4,1)*energy'/4; 
    energy_nod = sparse(iEner(:),1,sEner(:)); 
    comp = F'*U; 
end