%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-12-07 09:13:42
 % @LastEditTime: 2021-12-31 09:54:48
%--------------------------
function [dc, dv] = sensitivity_cvt_align(c_dH, c_interEles, area_cell, rho_cell, psi_cell, sK_cell, ...
    edofMat, U, e0, e1, penal, area0)

    ncvt = size(rho_cell, 1);
    dc_local = cell(ncvt, 1);

    % pn = 16; % p-norm for volume fraction
    % dpnorm = 1/ncvt * (mean(v_local.^pn))^(1/pn-1) * (v_local.^(pn-1)); % [N, 1]

    parfor ele = 1 : ncvt
        %% compliance sensitivity
        psi1 = psi_cell{ele};          % [6, 2v, m]
        psi2 = permute(psi1, [2,1,3]); % [2v, 6, m]
        nvdofs = size(psi2, 1);

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

        sK1 = sK0 .* (penal * rho(:)'.^(penal-1) * (e1 - e0));
        sK2 = reshape(sK1, 6,6,[]); % [6,6,m]

        Ke0 = mtimesx(mtimesx(psi2, sK2), psi1); % [2v,2v,m]
        Ke1 = reshape(Ke0, nvdofs, []); % [2v, 2v*m]

        edof = edofMat{ele};
        ue = U(edof); % [2v,1]

        dc1 = ue' * Ke1; % [1, 2v*m]
        dc2 = reshape(dc1, nvdofs, [])'; % [m,2v]
        dc_local{ele} = dc2 * ue; % [m, 1]
    end

    dcdp = zeros(ncvt, 3);
    dvdp = zeros(ncvt, 3);

    parfor i = 1 : ncvt
        for j = 1:3

            dH_ele = c_dH{i,j};
            interEles = c_interEles{i,j};
            nnb = length(interEles);

            for k = 1 : nnb
                ele = interEles(k);
                dcdp(i, j) = dcdp(i, j) + sum(dc_local{ele} .* dH_ele{k}); 
                dvdp(i, j) = dvdp(i, j) + sum(area_cell{ele} / area0 .* dH_ele{k}); 
            end
        end
    end
    dc = reshape(dcdp', [], 1); % TODO: NOT coorespond
    dv = reshape(dvdp', [], 1);
end