% Many thanks to Sanjay for setting me straight (again!) on
% tensors.
%  From: Sanjay Govindjee <sanjay@gough.ce.berkeley.edu>
%  Ok here is my run for E110 which does work out: -sg.

% stiffness matrix
c = [ 165 64 64 0 0 0;
64 165 64 0 0 0 ;
64 64 165 0 0 0 ;
0 0 0 80 0 0 ;
0 0 0 0 80 0;
0 0 0 0 0 80]

% (Define loading direction)
dir = [1;1;0]

% get a unit vector in that direction
v = dir/sqrt(dir'*dir);

% (Get rank 2 Tensor for for the stress)

s = v*v';

% s = (for 110 direction)
%   0.5000    0.5000         0
%   0.5000    0.5000         0
%        0         0         0

%(Map rank 2 tensor into a 6x1 vector ie. Voigt Notation)
sv = [s(1,1) ; s(2,2) ; s(3,3); s(1,2) ; s(2,3) ; s(3,1)];

% sv = (for 110 direction)
%   0.5000
%   0.5000
%        0
%   0.5000
%        0
%        0

% (Get strain as a 6x1 vector)

inv(c)*sv;

%ans = (for 110 direction)
%    0.0028
%    0.0028
%   -0.0022
%    0.0063
%         0
%         0

% (Map 6x1 Voigt strain vector into rank 2 tensor form; make sure you
% map engineering shear strains to tensorial shear strains)

eps = [ ans(1) ans(4)/2.0 ans(6)/2.0;
ans(4)/2.0 ans(2) ans(5)/2.0;
ans(6)/2.0 ans(5)/2.0 ans(3)];

%eps =
%    0.0028    0.0031         0
%    0.0031    0.0028         0
%         0         0   -0.0022

% (extract relevant strain)

ans = v'*eps*v;

% (Divide into stress to get modulus)

1/ans

% ans = (for 110 direction)
%   169.1240