%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-11-26 18:33:42
 % @LastEditTime: 2022-04-29 15:27:29
%--------------------------
function [psi_cell] = preparation_cvt(e1, e0, penal, nbnodes_cell, phi_cell, dofid_cell, ...
    rho_cell, faces_cell, c_nodes, sK_cell, c_edofMat, cbns)

    ncvt = size(rho_cell, 1);

    psi_cell = cell(ncvt, 1);
    c_NH = cell(ncvt, 1);

    % plot ele = 44
    parfor ele = 1 : ncvt
        
        faces = faces_cell{ele};
        nele  = size(faces, 1);
        [iK, jK, edofMat, ~, ~] = forAssemble_trimesh(faces);
        order1 = reshape(edofMat', [], 1);

        rho = rho_cell{ele}(:); % rho(:) = 1;
        phi = phi_cell{ele};
        sK0 = sK_cell{ele};
        dofid  = dofid_cell{ele};

        nedges = length(nbnodes_cell{ele});
        nbdofs = (sum(nbnodes_cell{ele}) - nedges) * 2;

        k  = assemble_micro_k_hetero(dofid, rho, e0, e1, penal, sK0, iK, jK);
        Me = build_M_simply(nbdofs, k, phi);

        Re = [phi; Me]; % [2n, 2v] 
        R1 = Re(dofid(order1), :); % [6m, 2v]  
        R2 = reshape(R1, 6, nele, []); % [6, m, 2v]
        R3 = permute(R2, [1,3,2]);

%         edof = c_edofMat{ele};
%         u = reshape(cbns', [], 1);
%         ue = u(edof);
%         umi = Re * ue;
%         fineNodes = reshape(umi(dofid), 2, [])';
%         figure; scatter(fineNodes(:,1), fineNodes(:,2)); hold on;
%         re = reshape(ue, 2, [])';
%         scatter(re(:,1), re(:,2), 'filled'); 
% 
%         dt = triangulation(faces, fineNodes);
%         triplot(dt);
        
        psi_cell{ele} = R3;   

        %% surface of NH
%         nodes = c_nodes{ele};
%         Nh = init_Nh_triangle(nodes, faces); % [m,3]
% 
%         n_vdofs = size(Re, 2);
%         
%         R3 = Re(dofid, :); % [2n, 2v]
%         NHe = NH_alpha(edofMat, n_vdofs/2, Nh, R3, size(nodes, 1));
%         NH{ele, 1} = NHe; % [n,2,2,v]
% 
%         plot_NH_cvt(nodes, faces, NHe);
    end 
end

function plot_NH_cvt(nodes, faces, NH)

    % [n,2,2,v]
    x = nodes(:, 1); % [n, 1]
    % x = x(faces);    % [m, 3]
    y = nodes(:, 2); % [n, 1]
    % y = y(faces);    % [m, 3]
    
    for i = 1 :size(NH,4)
        figure;
        % for j = 1:3
            cv = NH(:,2,2,i) * 10;           

            surfir(x,y,cv,0);
%             scatter3(x, y, cv); 
%             colormap(jet); caxis([0,1]);
%             s.FaceColor = 'w';
%             s.EdgeColor = [0.2,0.2,0.2];

%             hold on;
            grid off; axis equal; axis off;
            set(gca,'FontSize',14,'FontWeight','bold');
        % end
    end
end