Jbjacob94 patch 1

ImageD11/eps_sig_solver_NEW.py → ImageD11/stress.py

@@ -40,14 +40,13 @@ def __init__(self,
         -----------
         cij (float)        : elastic constants of the crystal
         Cij_symmetry (str) : symmetry considered for the Stiffness and Compliance matrices. Should be one of the following:
-                               'cubic', 'trigonal_high', 'trigonal_low', 'tetragonal', 'hexagonal','orthorombic', 'monoclinic', 'triclinic'
+                               'cubic', 'trigonal_high', 'trigonal_low', 'tetragonal_high', 'tetragonal_low', 'hexagonal', 'orthorhombic', 'monoclinic', 'triclinic'
 
-        dzero_unitcell (array_like)   : Unstrained unit cell parameters [a, b, c, alpha,beta, gamma]
+        unitcell (array_like)   : Unstrained unit cell parameters [a, b, c, alpha,beta, gamma]
         UBI_list (list of 3x3 arrays) : List of real-space unit cell vectors (ubi in ImageD11).
         """
 
-        assert symmetry in ['cubic', 'trigonal_high','trigonal_low', 'tetragonal',
-                            'hexagonal','orthorombic', 'monoclinic', 'triclinic'], 'symmetry not recognized!'
+        assert symmetry in Cij_symmetry.keys(), 'symmetry not recognized!'
 
         self.phase_name = name
         self.unitcell = unitcell
@@ -80,16 +79,17 @@ def __init__(self,
         self.stress_unit = 'GPa'
         self.Cij_symmetry = Cij_symmetry[self.symmetry]
         self.Cij = self.form_stiffness_tensor()
-        self.Sij = np.linalg.inv(self.Cij)
+        if not np.alltrue(self.Cij == 0):
+            self.Sij = np.linalg.inv(self.Cij)
 
         self.UBIs = UBI_list
         self.F_list = None
 
 
 
-    def __str__(self):
-        return f"EpsSigSolver:\n phase name: {self.phase_name}\n reference unitcell: {self.unitcell}\n symmetry:" +\
-               f"{self.symmetry}\n unit:{self.stress_unit}\n Stiffness:\n {self.Cij}\n  n ubis: {len(self.UBIs)}"
+ #   def __str__(self):
+ #       return f"EpsSigSolver:\n phase name: {self.phase_name}\n reference unitcell: {self.unitcell}\n symmetry:" +\
+ #              f"{self.symmetry}\n unit:{self.stress_unit}\n Stiffness:\n {self.Cij}\n  n ubis: {len(self.UBIs)}"
 
     # Load / save / update functions for parameters (from former eps_sig_solver.py)
     ########################################
@@ -138,14 +138,18 @@ def form_stiffness_tensor(self):
         Cij : 6x6 matrix containing the elastic components
         """
         Cij = np.zeros((6,6))
-        pattern = self.Cij_symmetry  # pattern for the stiffness matrix
+        pattern = self.Cij_symmetry  # pattern for the stiffness matrix. From Mouhat & Coudert (2014). Necessary and sufficient elastic stability conditions in various crystal systems. Physical Review
+
 
         parlist = 'c11,c22,c33,c44,c55,c66,c12,c13,c14,c15,c16,c23,c24,c25,c26,c34,c35,c36,c45,c46,c56'.split(',')
 
         for i,parname in enumerate(parlist):
             v = self.parameterobj.parameters[parname]
             if v is None:
                 continue
+            if self.symmetry in ['hexagonal','trigonal_high','trigonal_low'] and parname == 'c66':
+                v = 0.5 * (self.parameterobj.parameters['c11'] - self.parameterobj.parameters['c12'])
+
             c_ij = np.where( np.abs(pattern) == i+1, np.sign(pattern) * v, 0 )
             Cij += c_ij
         return Cij 
@@ -274,7 +278,7 @@ def convert_tensors_to_vecs(self, output_format = 'voigt', debug=0):
         # select all strain and stress tensors list
         dnames = [attr for attr in dir(self) if any([attr.startswith(s) for s in ['eps','sigma']]) ]
          # filter out all data that begins with 'eps' or 'sigma' but are not strain or stress tensors  
-        dnames = [d for d in dnames if any([d.endswith(s) for s in ['Lab','Ref','d']]) ] 
+        dnames = [d for d in dnames if any([d.endswith(s) for s in ['Lab','Ref','_d']]) ] 
 
         # stop if no strain / stress tensor list fond
         if len(dnames) == 0:
@@ -581,29 +585,8 @@ def invariants(T):
     I3 = np.linalg.det(T)
 
     return I1, I2, I3
-
-
-def invariants_quantities(T):
-    """
-    compute relevant invariant parameters of strain / stress tensor T
-
-    Returns
-    --------
-    P (float)     : -I1/3 : hydrostatic pressure / volumetric strain
-    τ_oct (float) : √(2*J2/3) : octahedral shear stress / strain
-    """
 
-    T_dev = deviatoric(T)
-    I1, I2, I3 = invariants(T)
-    J1, J2, J3 = invariants(T_dev)
-
-    return -I1/3, np.sqrt(2*J2/3)
-
-
-
-
-
-# Cij_symmetry for stiffness tensor (copied from matscipy.elasticity : https://github.com/libAtoms/matscipy/blob/master/matscipy/elasticity.py)
+# Cij_symmetry for stiffness tensor (from Mouhat & Coudert 2014)
 ########################
 
 Cij_symmetry = {
@@ -614,12 +597,12 @@ def invariants_quantities(T):
                                 [0, 0, 0, 0, 4, 0],
                                 [0, 0, 0, 0, 0, 4]]),
 
-   'trigonal_high':   np.array([[1, 7, 8, 9, 10, 0],
-                                [7, 1, 8, 0,-9, 0],
-                                [8, 8, 3, 0, 0, 0],
+   'trigonal_high':   np.array([[1, 7, 8, 9,  0, 0],
+                                [7, 1, 8, -9, 0, 0],
+                                [8, 8, 3, 0,  0, 0],
                                 [9, -9, 0, 4, 0, 0],
-                                [10, 0, 0, 0, 4, 0],
-                                [0, 0, 0, 0, 0, 6]]),
+                                [0,  0, 0, 0, 4, 9],
+                                [0, 0, 0, 0,  9, 6]]),
 
    'trigonal_low':    np.array([[1,  7,  8,  9,  10,  0 ],
                                 [7,  1,  8, -9, -10,  0 ],
@@ -665,5 +648,6 @@ def invariants_quantities(T):
    }
 
 
-Cij_symmetry['hexagonal'] = Cij_symmetry['trigonal_high']
+Cij_symmetry['hexagonal'] = Cij_symmetry['tetragonal_high']
+Cij_symmetry['tetragonal'] = Cij_symmetry['tetragonal_high']
 Cij_symmetry[None] = Cij_symmetry['triclinic']