function pk1_stress(F, global_u, global_l)
    
    [globaly, globalx] = size(global_u);
    P = zeros(4,4,globalx * globaly);
    
    for gp = 1:4
        Fe  = F{gp, 1};                    % 4 * m
        P(gp, : ,:)  = PK1_fast(global_u, global_l, Fe); % 4 * m
    end 
    Pm = squeeze(mean(P, 1));
    s1 = reshape(Pm(1, :), globaly, globalx);
    s2 = reshape(Pm(2, :), globaly, globalx);
    s3 = reshape(Pm(4, :), globaly, globalx);
    figure; imagesc(s1); colormap(jet); axis equal; axis off; %caxis([-0.1, 0.1]); 
    figure; imagesc(s2); colormap(jet); axis equal; axis off; %caxis([-0.15, 0.03]); 
    figure; imagesc(s3); colormap(jet); axis equal; axis off; %caxis([-0.1, 0.1]); 
end