function [sigma, d] = model(F)
%
% Let W = 1/2 * K * (J-1)^2 + mu/2 (tr C - 3) [terrible model!]
%

% Hardwire
K = 999;
mu = 11;

b = F*F';
J = det(F);

up = K*(J-1);
upp = K;

% Cauchy stress
sigma(1) = up + (1/J)*mu*b(1,1);
sigma(2) = up + (1/J)*mu*b(2,2);
sigma(3) = up + (1/J)*mu*b(3,3);
sigma(4) =      (1/J)*mu*b(1,2);
sigma(5) =      (1/J)*mu*b(2,3);
sigma(6) =      (1/J)*mu*b(3,1);


ioi =  [1 1 1 0 0 0;
        1 1 1 0 0 0;
        1 1 1 0 0 0;
        0 0 0 0 0 0;
        0 0 0 0 0 0;
        0 0 0 0 0 0];

Isym = [ 1 0 0 0 0 0;
         0 1 0 0 0 0;
         0 0 1 0 0 0;
         0 0 0 1/2 0 0;
         0 0 0 0 1/2 0;
         0 0 0 0 0 1/2];

% spatial moduli
d = (J*upp + up)*ioi - up*2*Isym;

