%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-05-02 20:28:27
 % @LastEditTime: 2022-05-03 14:46:16
%--------------------------
% function pure_CVT()
    nelx = 4;
    nely = 2; 
    nelz = 1;
    microx = 64;
    dx = 1;

    %% init
    lx = dx * microx;
    globalx = nelx * microx;
    globaly = nely * microx;
    globalz = nelz * microx;
    maxx = globalx / 2; minx = -maxx;
    maxy = globaly / 2; miny = -maxy;
    maxz = globalz / 2; minz = -maxz;

%     [gx, gy, gz] = meshgrid(linspace(minx,maxx,8), linspace(miny,maxy,4), linspace(miny,maxy,2));
    gx = randi([minx, maxx],64,1);
    gy = randi([miny, maxy],64,1);
    gz = randi([minz, maxz],64,1);
    seeds = [gx(:), gy(:), gz(:)];

    npnts = size(seeds, 1);
    T = repmat(8, npnts, 1);
    [xval, xold1, xold2, xmin, xmax, low, upp, m_mma, c_mma, d_mma, a0_mma, ...
            a_mma] = MMA_paras(seeds, maxx, minx, maxy, miny, maxz, minz);

    loop = 0;
    while loop < 100
        loop = loop + 1;

        indices = 1 : int32(npnts);
        [c_edges, cenP, vedges, ~, c_nnbSeedIDX] = intersection_info(seeds, ...
            indices, nelx, nely, nelz, microx, 'cantilever'); 

        MMCs = update_MMCs(vedges, T, c_nnbSeedIDX);
        plot_mmc(MMCs, minx, maxx, miny, maxy, minz, maxz);

        [E, dE] = cvt_energy(seeds, cenP);

        f0val= E; 
        df0dx= dE; 
        fval = 0;
        dfdx = repmat(1, 1, size(seeds, 1)*3);
        
        [xmma, ~, ~, ~, ~, ~, ~, ~, ~, low, upp] = mmasub_mmc(1, npnts * 3, ...
            loop, xval, xmin, xmax, xold1, xold2, f0val, df0dx, fval, dfdx,...
            low, upp, a0_mma, a_mma, c_mma, d_mma);
        xold2 = xold1; xold1 = xval; 
        xval = xmma;
            
        seeds = reshape(xmma, 3, [])';
    end
    
    write_cin(seeds, 'box'); 
% end

function [xval, xold1, xold2, xmin, xmax, low, upp, m_mma, c_mma, d_mma, a0_mma, ...
    a_mma] = MMA_paras(seeds, maxx, minx, maxy, miny, maxz, minz)

    np = size(seeds, 1);

    xval = reshape([seeds]', [], 1); 
    xold1= xval; xold2 = xval; 

    xmin = repmat([minx, miny, minz], np, 1)'; 
    xmax = repmat([maxx, maxy, maxz], np, 1)'; 
    xmin = xmin(:);
    xmax = xmax(:);

    low  = xmin; upp = xmax; 

    m_mma = 1; % number of constrain
    c_mma = 1000*ones(m_mma,1); d_mma = zeros(m_mma,1); a0_mma = 1; a_mma = zeros(m_mma,1); 
end


function [E, dE] = cvt_energy(seeds, cenP)

    npnts   = size(seeds, 1);
    E0  = sum((seeds - cenP) .^ 2, 2); % [n,1]
    dE0 = 2 * (seeds - cenP);          % [n,2]
    
    E  = sum(E0);
    dE = dE0;
    dE = reshape(dE', [], 1);
end