function [umicro, dudU, c] = microFEA_FE2(edofMat, nelx, nely, KE, iK, jK, freedofs, fixeddofs, ...
                                            ufixed, xPhys, E0, Emin)

    sK = reshape(KE(:) * (Emin + xPhys(:)' * (E0 - Emin)), 64 * nelx * nely, 1);
    K = sparse(iK, jK, sK); 
    K = (K + K') / 2;
    
    umicro = ufixed;
    
    k21 = K(freedofs, fixeddofs);
    k22 = K(freedofs, freedofs);

    Fr = -k21*ufixed(fixeddofs);
    umicro(freedofs) = k22\Fr;

    % dudU includes two parts: du_free/dU, dUdU = 1
    dudU = ones(2*(nelx+1)*(nely+1), length(fixeddofs));
    dudU(freedofs, :) = -k22\k21;

    % energy
    ce = sum((umicro(edofMat)*KE) .* umicro(edofMat), 2);
    c = sum((Emin + xPhys(:) * (E0 - Emin)) .* ce);
end

% % regard macro_U as micro BC
% function [Ue] = microFEA_FE2(ke, iK, jK, ufixed, freedofs, fixeddofs, xPhys, e, emin)    

%     if nargin>9
%         sK = ke(:) * (emin + xPhys(:)'*(e-emin)); 
%     else
%         sK = ke(:) * xPhys(:)'; % xPhys = sum(H(De)^2)/4, thus not need penalty
%     end
%     K = sparse(iK, jK, sK(:)); 
%     K = (K + K') / 2;  

%     %% FE2 FEA
%     fr = - K(freedofs, fixeddofs) * ufixed(fixeddofs, 1);
%     Ue(fixeddofs, 1) = ufixed(fixeddofs, 1);
%     Ue(freedofs, 1) = K(freedofs, freedofs) \ fr;
% end