for step=1:Nstep,
    for k =1:Nbel,
        F11=DEF_GRADIENT(k,step);
        F22=DEF_GRADIENT(Nbel+k,step);
        F12=DEF_GRADIENT(2*Nbel+k,step);
        F21=DEF_GRADIENT(3*Nbel+k,step);
        F33=1/(F11*F22-F12*F21);
        F=[F11 F12 0; F21 F22 0; 0 0 F33];
        sigma=2*F*F'-2*F33^2*eye(3);
        sigma11(k,step)=sigma(1,1);
        sigma22(k,step)=sigma(2,2);
        sigma12(k,step)=sigma(1,2);
        EV11(k,step)=EPSILON_VIRT(k,step);
        EV22(k,step)=EPSILON_VIRT(Nbel+k,step);
        EV12(k,step)=EPSILON_VIRT(2*Nbel+k,step);   
        Wint(k,step)=(EV11(k,step)*sigma11(k,step)+...
            EV22(k,step)*sigma22(k,step)+...
            2*EV12(k,step)*sigma12(k,step))*surf_fin(k,step)*F33;
    end
end
Wint_tot=sum(Wint);
Wext_tot=P'.*sum(surf_fin);

C10=10;
figure;plot(C10*Wint_tot(25:end));hold on;plot(Wext_tot(25:end));
% 
% Myplot(xGfin,yGfin,zGfin,sigma11(:,end))
% Myplot(xGfin,yGfin,zGfin,sigma22(:,end))
% Myplot(xGfin,yGfin,zGfin,sigma12(:,end))