% J2 plasticity material subroutine with isotropic hardening.
function [sv, C] = get_stress(sv, Ce, H, deps)
    % sv is a structure w/ two fields
    % sv.sig is a 4x1 vector of the stress components  [sxx, syy, sxy, szz].
    % sv.k is the current value of the yield stress.

    % Compute trial stress (Ce is given an extra row to update szz).
    sv.sig = sv.sig + [Ce; Ce(1,2), Ce(1,2), 0]*deps;
    % Compute pressure and deviatoric stress.
    p = -(sv.sig(1) + sv.sig(2) + sv.sig(4))/3.0;
    s_tr = sv.sig + p*[1;1;0;1];

    % Deviatoric stress norm; include sxy twice due to voigt notation.
    a = norm([s_tr; s_tr(3)]);
    if a < sv.k,
        C = Ce;
    else
        % Solve for plastic increment.
        nhat = s_tr/a;
        mu = Ce(3,3);
        dlambda = (a-sv.k) / (2.0*mu+H);

        % Conversion from engineering strain.
        s = diag([2,2,1,2]);
        % Update stress and yield strength.
        sv.sig = sv.sig - s*mu*dlambda*nhat;
        sv.k = sv.k + H*dlambda;

        assert(abs(nhat(1)+nhat(2)+nhat(4)) < 1e-14, 'nhat not deviatoric');

        NN = nhat(1:3)*nhat(1:3)';
        s = s(1:3,1:3);
        % Compute algorithmic tangent modulus.
        Cep = Ce - s*mu*NN / (1.0+H/(2.0*mu));
        gamma = s*mu*dlambda/a;
        D = eye(3) - [1;1;0]*[1,1,0]/3.0 - NN;
        C = Cep - s*mu*gamma*D;
        %C = Ce;
    end
end