%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-04-15 15:10:00
 % @LastEditTime: 2021-04-15 16:12:43
%--------------------------
function plot_stress(datapath, U, D0, D1, dx, globalx, globaly, deform_nodes, global_rho)

    [~, dNh, ~] = shape_gradient_ele(dx/2, dx/2, 4);
    [~, ~, edofMat] = forAssemble(globalx, globaly);

    defn = reshape(deform_nodes', [], 1);
    nele = globalx * globaly;
    global_cx = zeros(4, nele);
    global_cy = zeros(4, nele);
    
    % global_info = zeros(4, nele);
    def_info1 = zeros(4, nele);
    def_info2 = zeros(4, nele);
    def_info3 = zeros(4, nele);

    stress_field1 = zeros(globaly, globalx);
    stress_field2 = zeros(globaly, globalx);
    stress_field3 = zeros(globaly, globalx);
    
    strain_field1 = zeros(globaly, globalx);
    strain_field2 = zeros(globaly, globalx);
    strain_field3 = zeros(globaly, globalx);

    for ele = 1 : globalx * globaly
        edof = edofMat(ele, :);
        ue = U(edof);

        if global_rho(ele)==1
            D = D1;
        else
            D = D0;
        end

        s = zeros(4,3);
        se = zeros(4,3);
        for gp = 1:4
            dNh_g = dNh{gp}; % 3,8
            s(gp, :)  = D * dNh_g * ue;
            se(gp, :) = dNh_g * ue;
        end

        coor = reshape(defn(edof), 2, [])'; % [4,2]        
        global_cx(:, ele) = coor(:,1);
        global_cy(:, ele) = coor(:,2);

        def_info1(:, ele) = se(:, 1);
        def_info2(:, ele) = se(:, 2);
        def_info3(:, ele) = se(:, 3);

        stress_field1(ele) = mean(s(:, 1));
        stress_field2(ele) = mean(s(:, 2));
        stress_field3(ele) = mean(s(:, 3));
        
        strain_field1(ele) = mean(se(:, 1));
        strain_field2(ele) = mean(se(:, 2));
        strain_field3(ele) = mean(se(:, 3));
    end
%     figure; imagesc(stress_field1); colormap(jet); caxis([-10, 10]); axis equal; axis off;
%     figure; imagesc(stress_field2); colormap(jet); caxis([-10, 10]); axis equal; axis off;    
%     figure; imagesc(stress_field3); colormap(jet); caxis([-5, 5]); axis equal; axis off;

    figure; imagesc(strain_field1); colormap(jet); axis equal; axis off; %caxis([-0.1, 0.1]); 
    figure; imagesc(strain_field2); colormap(jet); axis equal; axis off; %caxis([-0.05, 0.05]); 
    figure; imagesc(strain_field3); colormap(jet); axis equal; axis off; %caxis([-0.1, 0.1]); 
    
%     figure;
%     h = fill(global_cx, global_cy, def_info1); axis equal; axis off; 
%     set(h, {'LineStyle'}, {'none'});
%     colormap('jet'); caxis([-0.1, 0.1]);
%     
%     h = gcf;
%     set(h,'Units','Inches');
%     pos = get(h,'Position');
%     set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
%     print(h, datapath,'-dpdf','-r0')

%     figure;
%     h = fill(global_cx, global_cy, def_info2); axis equal; 
%     set(h, {'LineStyle'}, {'none'});
%     colormap('jet'); caxis([-0.05, 0.05]);
% 
%     figure
%     h = fill(global_cx, global_cy, def_info3); axis equal; 
%     set(h, {'LineStyle'}, {'none'});
%     colormap('jet'); caxis([-0.1, 0.1]); 
end