%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% matlab example 01 - tensor calculus %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
% based on a give deformation gradient, F = grad_X phi, evaluate 
% constitutive equations, calculate relevant quantities, play around ;-)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
ndim = 3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prescribe deformation gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(ndim,ndim);
alpha = 0.01; gamma = 0.2; theta = pi/4;

F=[1,0,0; 0,1,0; 0,0,1];                                      % unit map
F=[1.001,0,0.02; 0,1.0004,0; 0,0,1];                 % small deformation
F=[alpha,0,0; 0,sqrt(alpha),0; 0,0,sqrt(alpha)];  % uniax stretch incomp
F=[1,gamma,0; 0,1,0; 0,0,1];                                     % shear
%F=[cos(theta),sin(theta),0; -sin(theta),cos(theta),0; 0,0,1]; % rotation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material parameter E, nue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
emod = 1.0;
nue  = 0.2;

lam = 0.577;
mue = 0.385;
kap = emod / (3.0*(1.0-2.0*nue));
lam = emod * nue / (1.0+nue)/ ( 1.0 - 2.0 * nue );
mue = emod / (2.0*(1.0+    nue));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% jacobi determinant                  det_F=det(F)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
det_F = F(1,1)*F(2,2)*F(3,3)+F(2,1)*F(3,2)*F(1,3)+F(3,1)*F(1,2)*F(2,3) ...
	  - F(1,1)*F(2,3)*F(3,2)-F(2,2)*F(3,1)*F(1,3)-F(3,3)*F(1,2)*F(2,1);

det_F = det(F); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% inverses                            F_inv = F^-1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I=zeros(ndim,ndim);
F_inv = inv(F);
for i=1:3
  for j=1:3
    for k=1:3
	  I(i,j) = I(i,j) + F(i,k) * F_inv(k,j);
	end
  end
end  	  

I = F*F_inv;
I = eye(ndim);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% left and right cauchy green tensor     C = F^t * F      b = F * F^t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(ndim,ndim);
b = zeros(ndim,ndim);
for i=1:ndim
  for j=1:ndim
	for k=1:ndim
	   C(i,j) = C(i,j) + F(k,i) * F(k,j);
	   b(i,j) = b(i,j) + F(i,k) * F(j,k);
	end
  end
end

C = F'*F;
b = F*F';
C_inv = inv(C);
b_inv = inv(b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% green lagrange strain tensor        E = 1/2 [C -I]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(ndim,ndim);
for i=1:ndim
  for j=1:ndim
	   E(i,j) = E(i,j) + 1/2 *(C(i,j) - I(i,j));
  end
end

E =  (C-I)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% invariants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C2 = C*C;
b2 = b*b;

tr_C  =  C(1,1) + C(2,2) + C(3,3);
tr_b  =  b(1,1) + b(2,2) + b(3,3);
tr_C  = trace(C);
tr_b  = trace(b);

I1 = trace(C);
I2 = 1/2*((trace(C))^2 - trace(C2));
I3 = det(C);

I1 = trace(b);
I2 = 1/2*((trace(b))^2 - trace(b2));
I3 = det(b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fourth order unit tensors           I_sym, I_skw, I_vol, I_dev
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_sym = zeros(ndim,ndim,ndim,ndim);
I_skw = zeros(ndim,ndim,ndim,ndim);
I_vol = zeros(ndim,ndim,ndim,ndim);
I_dev = zeros(ndim,ndim,ndim,ndim);

for i=1:ndim
  for j=1:ndim
    for k=1:ndim
	  for l=1:ndim
	    I_sym(i,j,k,l) = 1/2 * (I(i,k) * I(j,l) + I(i,l) * I(j,k));
        I_skw(i,j,k,l) = 1/2 * (I(i,k) * I(j,l) - I(i,l) * I(j,k));
        I_vol(i,j,k,l) = 1/3 *  I(i,j) * I(k,l);
        I_dev(i,j,k,l) =-1/3 *  I(i,j) * I(k,l) ...
		                +1/2 * (I(i,k) * I(j,l) + I(i,l) * I(j,k));
	  end
	end
  end
end  	  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% projections of E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E_sym = zeros(ndim,ndim);
E_skw = zeros(ndim,ndim);
E_vol = zeros(ndim,ndim);
E_dev = zeros(ndim,ndim);

for i=1:ndim
  for j=1:ndim
    for k=1:ndim
	  for l=1:ndim
		E_sym(i,j) = E_sym(i,j) + I_sym(i,j,k,l) * E(k,l);
		E_skw(i,j) = E_skw(i,j) + I_skw(i,j,k,l) * E(k,l);
		E_vol(i,j) = E_vol(i,j) + I_vol(i,j,k,l) * E(k,l);
		E_dev(i,j) = E_dev(i,j) + I_dev(i,j,k,l) * E(k,l);
	  end
	end
  end
end  	  

E_sym = (E+E')/2;
E_skw = (E-E')/2;
E_skw =  E-E_sym;
E_vol =  I*trace(E)/3;
E_dev =  E-E_vol;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% energy for neo hooke material
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NRG = 1/2 * lam * log(det_F)* log(det_F) ...
    + 1/2 * mue * (trace(C)-ndim-2*log(det_F));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stress tensors P and S for neo hooke material 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P = zeros(ndim,ndim);
S = zeros(ndim,ndim);

for i=1:ndim
  for j=1:ndim
      P(i,j) = (lam*log(det_F)-mue)*F_inv(j,i)+mue*F(i,j);
	  S(i,j) = (lam*log(det_F)-mue)*C_inv(j,i)+mue*I(i,j); 
	sig(i,j) =((lam*log(det_F)-mue)*    I(j,i)+mue*b(i,j))/(det_F); 
  end
end	  

S = F_inv*P;

	  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% material operators A_mat & C_mat for neo hooke material
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A_mat = zeros(ndim,ndim,ndim,ndim);
C_mat = zeros(ndim,ndim,ndim,ndim);
	  
for i=1:ndim
  for j=1:ndim
    for k=1:ndim
	  for l=1:ndim
		A_mat(i,j,k,l) =     lam            * F_inv(j,i)*F_inv(l,k) ...
		               + mue                *     I(i,k)*    I(j,l) ...
					   +(mue-lam*log(det_F))*(F_inv(l,i)*    F(j,k));
		C_mat(i,j,k,l) =     lam            * C_inv(i,j)*C_inv(k,l) ...
					   +(mue-lam*log(det_F))*(C_inv(i,k)*C_inv(j,l) ...
											 +C_inv(i,l)*C_inv(j,k)); 
	  end
	end
  end
end	  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stress tensor S for neo hooke material @ small strains
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S = zeros(ndim,ndim);
for i=1:ndim
  for j=1:ndim
    for k=1:ndim
	  for l=1:ndim
	    S(i,j)=S(i,j)+C_mat(i,j,k,l)*E(k,l);
	  end
	end
  end
end	  	  
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% transformation strain tensor intofour_ten voigt notation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
strain_ten = E;

strain_vgt(1) =   strain_ten(1,1);
strain_vgt(2) =   strain_ten(2,2);
strain_vgt(3) =   strain_ten(3,3);
strain_vgt(4) = 2*strain_ten(1,2);
strain_vgt(5) = 2*strain_ten(2,3);
strain_vgt(6) = 2*strain_ten(3,1);	  

strain_vgt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% transformation 4th order tensor into voigt notation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
four_ten = I_vol;

four_vgt(1,1) = four_ten(1,1,1,1); four_vgt(1,2) = four_ten(1,1,2,2); 
four_vgt(1,3) = four_ten(1,1,3,3); four_vgt(1,4) = four_ten(1,1,1,2); 
four_vgt(1,5) = four_ten(1,1,2,3); four_vgt(1,6) = four_ten(1,1,3,1);

four_vgt(2,1) = four_ten(2,2,1,1); four_vgt(2,2) = four_ten(2,2,2,2); 
four_vgt(2,3) = four_ten(2,2,3,3); four_vgt(2,4) = four_ten(2,2,1,2); 
four_vgt(2,5) = four_ten(2,2,2,3); four_vgt(2,6) = four_ten(2,2,3,1);

four_vgt(3,1) = four_ten(3,3,1,1); four_vgt(3,2) = four_ten(3,3,2,2); 
four_vgt(3,3) = four_ten(3,3,3,3); four_vgt(3,4) = four_ten(3,3,1,2); 
four_vgt(3,5) = four_ten(3,3,2,3); four_vgt(3,6) = four_ten(3,3,3,1);

four_vgt(4,1) = four_ten(1,2,1,1); four_vgt(4,2) = four_ten(1,2,2,2); 
four_vgt(4,3) = four_ten(1,2,3,3); four_vgt(4,4) = four_ten(1,2,1,2); 
four_vgt(4,5) = four_ten(1,2,2,3); four_vgt(4,6) = four_ten(1,2,3,1);

four_vgt(5,1) = four_ten(2,3,1,1); four_vgt(5,2) = four_ten(2,3,2,2); 
four_vgt(5,3) = four_ten(2,3,3,3); four_vgt(5,4) = four_ten(2,3,1,2); 
four_vgt(5,5) = four_ten(2,3,2,3); four_vgt(5,6) = four_ten(2,3,3,1);

four_vgt(6,1) = four_ten(3,1,1,1); four_vgt(6,2) = four_ten(3,1,2,2); 
four_vgt(6,3) = four_ten(3,1,3,3); four_vgt(6,4) = four_ten(3,1,1,2); 
four_vgt(6,5) = four_ten(3,1,2,3); four_vgt(6,6) = four_ten(3,1,3,1);

four_vgt;
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% double scalar product in voigt notation 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
strain_vol = four_vgt * strain_vgt';