%
% Defining isotropic compliance and stiffness tensors using
% complianceTensor  and  stiffnessTensor introduced in MTEX 5
%
% David Mainprice 9/03/2018
%
% Data from
% J.Fortin, Y.Gueguen, and A.Schubnel (2007)
% Effects of pore collapse and grain crushing on ultrasonic velocities and
% Vp/Vs, J. Geophys. Res., 112, B08207, doi:10.1029/2005JB004005.
%

% sandstone constants
Eo = 42.13; % Youngs Modulus in GPa
Ko = 21.3;  % Bulk Modulus in GPa
Go = 18.0;  % Shear Modulus in GPa
vo=(3*Ko-2*Go)/(2*(3*Ko+Go)); % Poissons ratio

Isotropic tensor compliance Sij defined by Eo and vo

S11=(1/Eo); S12=(-vo/Eo); S44=2*(S11-S12);
% define matrix M
 M =...
 [[  S11     S12    S12    0.0     0.0    0.0];...
 [   S12     S11    S12    0.0     0.0    0.0];...
 [   S12     S12    S11    0.0     0.0    0.0];...
 [   0.0     0.0    0.0    S44     0.0    0.0];...
 [   0.0     0.0    0.0    0.0     S44    0.0];...
 [   0.0     0.0    0.0    0.0     0.0    S44]];

CS_Tensor_Iso = crystalSymmetry('1',[1 1 1],...
 [ 90.0 90.0 90.0]*degree,'X||a*','Y||b',...
 'mineral','Isotropic tensor')

S_iso = complianceTensor(M,CS_Tensor_Iso)
 
CS_Tensor_Iso = crystalSymmetry  
 
  mineral           : Isotropic tensor  
  symmetry          : 1                 
  a, b, c           : 1, 1, 1           
  alpha, beta, gamma: 90°, 90°, 90°     
  reference frame   : X||a, Y||b*, Z||c*
 
 
S_iso = complianceTensor  
  unit            : 1/GPa                                   
  rank            : 4 (3 x 3 x 3 x 3)                       
  doubleConvention: true                                    
  mineral         : Isotropic tensor (1, X||a, Y||b*, Z||c*)
 
  tensor in Voigt matrix representation: *10^-3
 23.736 -4.043 -4.043      0      0      0
 -4.043 23.736 -4.043      0      0      0
 -4.043 -4.043 23.736      0      0      0
      0      0      0 55.558      0      0
      0      0      0      0 55.558      0
      0      0      0      0      0 55.558

Isotropic tensor stiffness Cij defined by Ko and Go

C11 = Ko+(4/3)*Go ; C12=C11-2*Go; C44=(C11-C12)/2;
% define matrix M
 M =...
 [[  C11     C12    C12    0.0     0.0    0.0];...
 [   C12     C11    C12    0.0     0.0    0.0];...
 [   C12     C12    C11    0.0     0.0    0.0];...
 [   0.0     0.0    0.0    C44     0.0    0.0];...
 [   0.0     0.0    0.0    0.0     C44    0.0];...
 [   0.0     0.0    0.0    0.0     0.0    C44]];
CS_Tensor_Iso = crystalSymmetry('1',[1 1 1],...
 [ 90.0 90.0 90.0]*degree,'X||a*','Y||b',...
 'mineral','Isotropic tensor')
C_iso_solid = stiffnessTensor(M,CS_Tensor_Iso)
% Calculate  Compliance tensor by inversion
S_iso_solid = inv(C_iso_solid)
 
CS_Tensor_Iso = crystalSymmetry  
 
  mineral           : Isotropic tensor  
  symmetry          : 1                 
  a, b, c           : 1, 1, 1           
  alpha, beta, gamma: 90°, 90°, 90°     
  reference frame   : X||a, Y||b*, Z||c*
 
 
C_iso_solid = stiffnessTensor  
  unit   : GPa                                     
  rank   : 4 (3 x 3 x 3 x 3)                       
  mineral: Isotropic tensor (1, X||a, Y||b*, Z||c*)
 
  tensor in Voigt matrix representation:
 45.3  9.3  9.3    0    0    0
  9.3 45.3  9.3    0    0    0
  9.3  9.3 45.3    0    0    0
    0    0    0   18    0    0
    0    0    0    0   18    0
    0    0    0    0    0   18
 
S_iso_solid = complianceTensor  
  unit            : 1/GPa                                   
  rank            : 4                                       
  rank            : 4 (3 x 3 x 3 x 3)                       
  doubleConvention: true                                    
  mineral         : Isotropic tensor (1, X||a, Y||b*, Z||c*)
 
  tensor in Voigt matrix representation: *10^-3
 23.735 -4.043 -4.043      0      0      0
 -4.043 23.735 -4.043      0      0      0
 -4.043 -4.043 23.735      0      0      0
      0      0      0 55.556      0      0
      0      0      0      0 55.556      0
      0      0      0      0      0 55.556