%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-09-27 15:20:17
 % @LastEditTime: 2021-09-29 15:23:19
%--------------------------
function [H, H2, dH] = mmc_project2H(lx, nelx, nely, microx, mmc_descriptor, MMCs, nodes_in_ele, coorx, coory, e1)
    dx = lx / microx;
    globalx = microx * nelx;
    globaly = microx * nely;

    epsilon = 4*dx;
    
    phi_max = repmat(-1e5, globaly+1, 1+globalx);
    phi_cell = cell(nely*nelx*20, 1); idx = 1;
    
    for i = 1 : nelx
        for j = 1 : nely
            mmc = mmc_descriptor{j,i};
            for ci = 1:20
                cx = mmc(ci, 1);
                cy = mmc(ci, 2);
                l  = mmc(ci, 3);
                sint = mmc(ci, 4);
                t = MMCs(ci,j,i);

                phi = phi_func(coorx, coory, sint, cx, cy, t, l);                  
                phi_max = max(phi_max, phi);

                phi_cell{idx} = phi;
                idx = idx + 1;
            end
            % phi_cell{j,i} = phi_max((j-1)*(microx+1)+1:j*(microx+1), (i-1)*(microx+1)+1:i*(microx+1));
        end
    end
    heav = Heaviside(phi_max, epsilon);

    H  = reshape(sum(heav(nodes_in_ele), 2) /4, globaly, globalx);
    H2 = reshape(sum(heav(nodes_in_ele).^2, 2)/4, globaly, globalx) * e1;    

    dH = sensitivies_dt(phi_cell, heav, lx, nelx, nely, microx, mmc_descriptor, MMCs, nodes_in_ele, coorx, coory);

end

function dH = sensitivies_dt(phi_cell, heav, lx, nelx, nely, microx, mmc_descriptor, MMCs, nodes_in_ele, coorx, coory)
    % compute dH2/dt = dH2 * dHdt
    dx = lx / microx;
    globalx = microx * nelx;
    globaly = microx * nely;

    epsilon = 4*dx;
    step = max(2*dx, 5e-3);

    dH = cell(nely, nelx, 20);

    for i = 1:nelx
        for j = 1:nely
            mmc = mmc_descriptor{j,i};           
            for ci = 1:20
                global_ci = ((i-1)*nely + j - 1)*20 + ci;

                cx = mmc(ci, 1);
                cy = mmc(ci, 2);
                l  = mmc(ci, 3);
                sint = mmc(ci, 4);

                t_plus = MMCs(ci,j,i) + step;
                phi_plus = phi_func(coorx, coory, sint, cx, cy, t_plus, l);
                for ck1 = 1 : global_ci-1
                    phi_max_plus = max(phi_cell{ck1}, phi_plus);
                end
                for ck2 = global_ci+1 : nelx*nely*20
                    phi_max_plus = max(phi_cell{ck2}, phi_plus);
                end

                t_minus = max(MMCs(ci,j,i) - step, 1e-3); 
                phi_minus = phi_func(coorx, coory, sint, cx, cy, t_minus, l);
                for ck1 = 1 : global_ci-1
                    phi_max_minus = max(phi_cell{ck1}, phi_minus);
                end
                for ck2 = global_ci+1 : nelx*nely*20
                    phi_max_minus = max(phi_cell{ck2}, phi_minus);
                end

                heav_plus  = Heaviside(phi_max_plus, epsilon);
                heav_minus = Heaviside(phi_max_minus, epsilon);
                dHdt = (heav_plus - heav_minus) / 2*step;

                dH_i = 2*heav .* dHdt;
                dH{j, i, ci} = reshape(sum(dH_i(nodes_in_ele), 2)/4, globaly, globalx);
            end
        end
    end
end

function H = Heaviside(phi, epsilon)

    alpha = 1e-3; 
    
    [y,x] = size(phi);
    H = zeros(x*y, 1);
    num_all=[1 : x*y]';
    
    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;
    
    H = reshape(H, y, x);
end

function phi = phi_func(coorx, coory, sint, cx, cy, t, l)
    cost = sqrt(abs(1-sint^2));
    
    % center coordinate of ellipse is (cx, cy)
    x1 =  cost*(coorx-cx) + sint*(coory-cy);
    y1 = -sint*(coorx-cx) + cost*(coory-cy);
    
    % NOTE: negative sign!!!
    phi = -((x1/l).^6 + (y1/t).^6 - 1);
end