%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-12-09 11:54:24
 % @LastEditTime: 2021-12-15 11:00:41

 % pb for blending
%--------------------------
function [rho_cell, dH_cell] = projection_mmc_cvt(nodes_cell, faces_cell, MMCs, epsilon, nd, pb)

%     nd = 3;

    ncvt = size(nodes_cell, 1);

    rho_cell = cell(ncvt, 1);
    dH_cell  = cell(ncvt, 1);

    % figure; clf
    for ele = 1 : ncvt
        nodes = nodes_cell{ele};
        faces = faces_cell{ele};
        
        mmc = MMCs{ele};
        nc  = size(mmc, 1);
        
        phi_max = repmat(-1e5, size(nodes, 1), 1);
%         phi0 = 0;
        phi_cell = cell(nc, 1);
        for i = 1 : nc
            phi = tPhi(mmc(i, :), nodes(:,1), nodes(:,2));
            phi_max = max(phi_max, phi);
            phi_cell{i} = phi;
%             phi0 = phi0 + exp(pb * phi);
        end
%         phi_max = log(phi0) / pb;

        dx = 4 * epsilon;
        H = Heaviside(phi_max, dx); 

        rho_cell{ele} = sum(H(faces), 2) / 4;

        % patch('Faces',faces,'Vertices',nodes,'FaceVertexCData',1 - H','FaceColor','interp','EdgeColor','interp'); hold on;

        %% sensitivity
        step = repmat(epsilon / 2, 1, nd);
        dH   = zeros(nc, nd, size(faces, 1));
        for i = 1 : nc
            for di = 1 : nd
                mmc_plus = mmc(i, :); 
                mmc_plus(di) = mmc(di) + step(di);
                tmpPhiD1 = tPhi(mmc_plus, nodes(:,1), nodes(:,2)); 
                phi_plus = tmpPhiD1; 
%                 phi_plus = exp(pb * tmpPhiD1);
                for ik = 1:i-1 
%                     phi_plus = phi_plus + exp(pb * phi_cell{ik});
                    phi_plus = max(phi_plus, phi_cell{ik}); 
                end
                for ik = i+1:nc
%                     phi_plus = phi_plus + exp(pb * phi_cell{ik});
                    phi_plus = max(phi_plus, phi_cell{ik}); 
                end
%                 phi_plus_blend = log(phi_plus) / pb;

                mmc_minus = mmc(i, :); 
                mmc_minus(di) = mmc(di) - step(di); 
                tmpPhiD2 = tPhi(mmc_minus, nodes(:,1), nodes(:,2)); 
                phi_minus=tmpPhiD2; 
%                 phi_minus = exp(pb * tmpPhiD2);
                for ik = 1:i-1 
%                     phi_minus = phi_minus + exp(pb * phi_cell{ik});
                    phi_minus = max(phi_minus, phi_cell{ik}); 
                end
                for ik = i+1:nc
%                     phi_minus = phi_minus + exp(pb * phi_cell{ik});
                    phi_minus = max(phi_minus, phi_cell{ik}); 
                end
%                 phi_minus_blend = log(phi_minus) / pb;

%                 HD1 = Heaviside(phi_plus_blend,  dx); 
%                 HD2 = Heaviside(phi_minus_blend, dx); 
                HD1 = Heaviside(phi_plus,  dx); 
                HD2 = Heaviside(phi_minus, dx); 
                rho1 = sum(HD1(faces), 2) / 4;
                rho2 = sum(HD2(faces), 2) / 4;
                dH(i, di, :) = (rho1 - rho2) / (2 * step(di));     
            end
        end
        dH_cell{ele} = dH;
    end
    % colormap("gray"); axis equal; axis off; pause(1e-6);
end

function H = Heaviside(phi, epsilon) 
    alpha=1e-3; % parameter alpha in the Heaviside function 

    num_all=[1 : length(phi)]'; 
    num1=find(phi>epsilon); 
    H(num1)=1; 
    num2=find(phi<-epsilon); 
    H(num2)=alpha; 
    num3=setdiff(num_all,[num1;num2]); 
    H(num3)=3*(1-alpha)/4*(phi(num3)/epsilon-phi(num3).^3/(3*(epsilon)^3))+ (1+alpha)/2; 
end