%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-11-28 13:49:21
 % @LastEditTime: 2022-02-01 12:02:45
%--------------------------
function K = stiffness_cvt(rho_cell, sK_cell, psi_cell, c_edofMat, nvar, e0, e1, penal, c_vol)

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

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

        sK0 = sK_cell{ele}; % [12^2, m]
        rho = rho_cell{ele};

        sK1 = sK0 .* (e0 + rho(:)'.^penal * (e1 - e0));
        sK2 = reshape(sK1, 12,12,[]); % [12,12,m]
        
        Ke0 = mtimesx(mtimesx(psi2, sK2), psi1); % [3v,3v,m]
        
%         vol0 = c_vol{ele}; % [m,1]
%         vol1 = zeros(1,1,length(vol0));
%         vol1(1,1,:) = vol0;
%         Ke0 = Ke0 .* vol1;
        
        Ke1 = squeeze(sum(Ke0, 3)); % [3v,3v], each ele is different

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

    % add small value to diagonal elements
    % in case of disconnected projected rho caused by thin components
    S = spdiags(repmat(1e-5, nvar, 1),0,nvar,nvar);
    K = K + S;

    if sum(isinf(K(:))) > 0 || sum(isnan(K(:))) > 0
        error("here");
    end
end