%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-11-28 13:49:21
 % @LastEditTime: 2021-12-06 16:52:55
%--------------------------
function K = stiffness_cvt(rho_cell, sK_cell, psi_cell, edofMat, nvar, e0, e1, penal)

    ncvt = size(psi_cell, 1);
    K = sparse(nvar, nvar);

    for ele = 1 : ncvt
        psi1 = psi_cell{ele};          % [6, 2v, m]
        psi2 = permute(psi1, [2,1,3]); % [2v, 6, m]

        sK0 = sK_cell{ele}; % [36, m]
        rho = rho_cell{ele};

        sK1 = sK0 .* (e0 + rho(:)'.^penal * (e1 - e0));
        sK2 = reshape(sK1, 6,6,[]); % [6,6,m]
        
        Ke0 = mtimesx(mtimesx(psi2, sK2), psi1); % [2v,2v,m]
        Ke1 = squeeze(sum(Ke0, 3)); % [2v,2v], each ele is different

        edof = edofMat{ele};
        K(edof, edof) = K(edof, edof) + Ke1;
    end
end