%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-06-16 15:55:30
 % @LastEditTime: 2021-07-01 16:16:54

 % [cx, cy] is the coordinate of the target coarse-node
%--------------------------
function plot_NH_unitp(cx, cy, lx, microx, nelx, nely, addNum, NH)
    
    % find node-id based on [cx, cy]
    [ref_nodes, ~, ~] = multi_generate_mesh(lx, nelx, nely, addNum);
    nid = project_force(cx, cy, ref_nodes);
    
    eid_global = nodes_eles_id(nelx, nely, 2, microx);
    eid = sort(eid_global{nid}); % acsending

    dx = lx / microx;
    [gridx, gridy] = meshgrid(0:dx:2*lx, 0:dx:2*lx);
    gridy = flipud(gridy);
    gridx = fliplr(gridx);

    for ii = 2
        for jj = 2
            
            figure;
            
            for i = 1:4
                e0 = eid(i);
                NHe = NH{e0,1}; % [(microx+1)^2, 2, 2, nnodes]
                
                switch i
                case 3
%                     NHi = NHe(:,:,:,2);
                    NHi = NHe(:,:,:,4);
                    gxi = gridx(1:microx+1, 1:microx+1);
                    gyi = gridy(1:microx+1, 1:microx+1);

                    gx_b = gridx(1:microx+1, microx+1);
                    gy_b = gridy(1:microx+1, 1);
                    cv_i = microx+1;
                    cv_j = fliplr(1:microx+1);

                    % to plot full-scale bilinear surf
                    cv = reshape(NHi(:,ii,jj), microx+1, microx+1);
                    cv(:) = 0;
                    cv(1, end) = 1;
                    cv = flipud(cv);
                    
                case 4
                    NHi = NHe(:,:,:,7);
                    gxi = gridx(1:microx+1, 1:microx+1);
                    gyi = gridy(microx+1:end, microx+1:end);
                   
                    gx_b = gxi(1, :);
                    gy_b = gyi(1, :);
                    cv_i = 1;
                    cv_j = fliplr(1:microx+1);

                    cv = reshape(NHi(:,ii,jj), microx+1, microx+1);
                    cv(:) = 0;
                    cv(end, end) = 1;
                    cv = flipud(cv);
                    
                case 1
                    NHi = NHe(:,:,:,1);
                    gxi = gridx(microx+1:end, microx+1:end);
                    gyi = gridy(1:microx+1, 1:microx+1);
                                       
                    gx_b = gxi(end, :);
                    gy_b = gyi(end, :);
                    cv_i = microx+1;
                    cv_j = fliplr(1:microx+1);

                    cv = reshape(NHi(:,ii,jj), microx+1, microx+1);
                    cv(:) = 0;
                    cv(1, 1) = 1;
                    cv = flipud(cv);
                    
                case 2
                    NHi = NHe(:,:,:,10);
                    gxi = gridx(microx+1:end, microx+1:end);
                    gyi = gridy(microx+1:end, microx+1:end);  
                    
                    gx_b = gxi(:, 1);
                    gy_b = gyi(:, 1);
                    cv_i = 1:microx+1;
                    cv_j = 1;

                    cv = reshape(NHi(:,ii,jj), microx+1, microx+1);
                    cv(:) = 0;
                    cv(end, 1) = 1;
                    cv = flipud(cv);
                end
                
%                 cv = reshape(NHi(:,ii,jj), microx+1, microx+1);
                s = surf(gxi, gyi, cv); 
                colormap(jet); caxis([0,1]);
                s.FaceColor = 'interp';
                s.EdgeColor = [0.2,0.2,0.2];
%                 s.EdgeColor = 'none';
                
                hold on;
                
                bsurf2 = cv(cv_i,cv_j);
                s1 = plot3(gx_b, gy_b, fliplr(bsurf2),'--');
                s1.LineWidth = 2;
                s1.Color = 'k';
            end            
        end
    end
    grid off; axis equal; axis off;
    set(gca,'FontSize',14,'FontWeight','bold');
end

% if vx / vy < 0, means taking all indices
function nodes2 = project_force(vx, vz, nodes)
    nodes(:,3) = 1:size(nodes,1);
    if vx > -1e-3 % valid
        v = min(abs(nodes(:,1) - vx));
        [ix, ~] = find(abs(nodes(:,1)-v-vx) < 1e-3);
        nodes1 = nodes(ix, :);
    else
        nodes1 = nodes;
    end

    if vz > -1e-3 % valid
        v = min(abs(nodes1(:,2) - vz));
        [iz, ~] = find(abs(nodes1(:,2)-v-vz) < 1e-3);
        nodes2 = nodes1(iz, 3);
    else
        nodes2 = nodes1(:,3);
    end
end