add frame_conventions to AnisotropicSolution

geodynamics · Mar 27, 2024 · 6ce3284 · 6ce3284
1 parent dbaa4c0
commit 6ce3284
Show file tree

Hide file tree

Showing 3 changed files with 10 additions and 8 deletions.
diff --git a/burnman/classes/anisotropicmineral.py b/burnman/classes/anisotropicmineral.py
@@ -372,11 +372,12 @@ def deformed_coordinate_frame(self):
         if self.orthotropic:
             return np.eye(3)
         else:
+            c = self.frame_convention
             M_T = self.unrotated_cell_vectors
             Q = np.empty((3, 3))
-            Q[0] = M_T[0] / np.linalg.norm(M_T[0])
-            Q[2] = np.cross(M_T[0], M_T[1]) / np.linalg.norm(np.cross(M_T[0], M_T[1]))
-            Q[1] = np.cross(Q[2], Q[0])
+            Q[c[0]] = M_T[c[0]] / np.linalg.norm(M_T[c[0]])
+            Q[c[2]] = np.cross(M_T[c[0]], M_T[c[1]]) / np.linalg.norm(np.cross(M_T[c[0]], M_T[c[1]]))
+            Q[c[1]] = np.cross(Q[c[2]], Q[c[0]])
             return Q
 
     @material_property
@@ -402,9 +403,10 @@ def cell_vectors(self):
         if self.orthotropic:
             return self.unrotated_cell_vectors
         else:
-            return np.einsum(
+            vectors = np.einsum(
                 "ij, jk->ik", self.unrotated_cell_vectors, self.rotation_matrix
             )
+            return vectors
 
     @material_property
     def cell_parameters(self):

diff --git a/burnman/classes/anisotropicsolution.py b/burnman/classes/anisotropicsolution.py
@@ -77,7 +77,7 @@ def __init__(
         self._logm_M0_mbr = np.einsum(
             "kij->ijk", np.array([logm(m[0].cell_vectors_0.T) for m in self.endmembers])
         )
-
+        self.frame_convention = self.endmembers[0][0].frame_convention
         self.anisotropic_params = anisotropic_parameters
         self.psi_excess_function = psi_excess_function
 

diff --git a/burnman/utils/unitcell.py b/burnman/utils/unitcell.py
@@ -100,9 +100,9 @@ def cell_vectors_to_parameters(vectors, frame_convention):
     """
 
     c = frame_convention
-    assert vectors[c[0], c[1]] == 0
-    assert vectors[c[0], c[2]] == 0
-    assert vectors[c[1], c[2]] == 0
+    assert np.abs(vectors[c[0], c[1]]) < np.finfo(float).eps
+    assert np.abs(vectors[c[0], c[2]]) < np.finfo(float).eps
+    assert np.abs(vectors[c[1], c[2]]) < np.finfo(float).eps
 
     lengths = np.empty(3)
     angles = np.empty(3)