%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-06-16 15:55:30
 % @LastEditTime: 2021-06-29 19:09:26

 % [cx, cy] is the coordinate of the target coarse-node
%--------------------------
function plot_bezier_surf(lx, microx, bezier_B, nelx, nely)
   
    dx = lx / microx;
    [gridx, gridy] = meshgrid(0:dx:2*lx, 0:dx:2*lx);
    gridy = flipud(gridy);
    % gridx = fliplr(gridx);
    
    [ref_nodes, ~, ~] = multi_generate_mesh(lx, 1, 1, 2);
    [~, ~, edofMat,  ~] = macro_edofMat_bezier(nelx, nely, 2); 
    nodes = reshape(ref_nodes', [], 1);
%     nodes = nodes(edofMat(1, :));
    boundary_nodes = bezier_B * nodes;
    
    figure;
    for i = [1:4]
%         nodei = nodes(edofMat(i, :));
        
        
        switch i 
        case 3 % left-up area, bottom boundary
            % in bezier_B
            coarse_nid = 4;            
            boundary_nid = 1 : microx+1;
            
            gx = gridx(microx+1, 1 : microx+1);
            gy = gridy(microx+1, 1 : microx+1);
            
        case 1 % left-bottom area, right boundary
            coarse_nid = 7;
            boundary_nid = microx+1 : 2*microx+1;
            
            gx = gridx(1:microx+1, microx+1);
            gy = gridy(1:microx+1, microx+1);
            
        case 2 % right-up area, left boundary
            coarse_nid = 1;
            boundary_nid = fliplr([3*microx+1 : 4*microx, 1]);

            gx = gridx(microx+1:2*microx+1, microx+1);
            gy = gridy(microx+1:2*microx+1, microx+1);
            
        case 4 % right-bottom area, up boundary
            coarse_nid = 10;
            boundary_nid = fliplr(2*microx+1 : 3*microx+1);

            gx = gridx(microx+1, microx+1:2*microx+1);
            gy = gridy(microx+1, microx+1:2*microx+1);
        end

%         bsurf2 = bezier_B(boundary_nid(:)*2, coarse_nid*2);
        bsurf2 = boundary_nodes(boundary_nid(:)*2);
%         scatter3(gx, gy, bsurf2, 3, bsurf2, 'filled'); colormap(jet);
        bsurf_dense = interp1(1:length(gx), bsurf2, 1:0.25:length(gx));
        gx_dense = interp1(1:length(gx), gx, 1:0.25:length(gx));
        gy_dense = interp1(1:length(gx), gy, 1:0.25:length(gx));
        scatter3(gx_dense, gy_dense, bsurf_dense, 3, bsurf_dense, 'filled'); colormap(jet);
%         s1 = plot3(gx, gy, bsurf2);
%         s1.Color = bsurf2;
%         s1.LineWidth = 2;
        hold on;
    end
%     scatter3(gx_dense, gy_dense, bsurf_dense, 2, bsurf_dense); colormap(jet);
    O = [0,0,0];
    L = [2*lx, 2*lx, 2*lx];
    s1 = plot3(O, 0:2*lx, O); hold on;
    s1.Color = [0         0    0.5156];
    s1.LineWidth = 2;
    s1 = plot3(L, 0:2*lx, O); hold on;
    s1.Color = [0         0    0.5156];
    s1.LineWidth = 2;
    s1 = plot3(0:2*lx, O, O); hold on;
    s1.Color = [0         0    0.5156];
    s1.LineWidth = 2;
    s1 = plot3(0:2*lx, L, O); hold on;
    s1.Color = [0         0    0.5156];
    s1.LineWidth = 2;
    
    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