+ handle rotation matrix ambiguities

src/openmc_cad_adapter/gqs.py

@@ -13,19 +13,70 @@
 PARABOLIC_CYLINDER = 9
 
 
+def _reorder_columns_for_orthogonality(A):
+    # Get the size of the matrix
+    n = A.shape[0]
+
+    # List to track used rows for the diagonal
+    used_rows = set()
+
+    # List to store the new column order
+    column_order = []
+
+    # Tolerance for ambiguity
+    tolerance = 1e-8
+
+    # Determine the best column for each diagonal position
+    for i in range(n):  # For each diagonal position
+        max_value = -np.inf
+        best_col = -1
+        ambiguous_cols = []
+
+        for j in range(n):  # Check each column
+            if j not in column_order:  # Only consider unused columns
+                value = abs(A[i, j])
+                if value > max_value + tolerance:
+                    max_value = value
+                    best_col = j
+                    ambiguous_cols = []  # Clear ambiguities
+                elif abs(value - max_value) < tolerance:  # Handle ambiguity
+                    ambiguous_cols.append(j)
+
+        if best_col != -1:
+            column_order.append(best_col)
+            used_rows.add(i)
+        elif ambiguous_cols:  # Handle ambiguous columns
+            ambiguous_cols.append(best_col)
+            # Check for orthogonality with already selected columns
+            for col1, col2 in zip(ambiguous_cols, ambiguous_cols[1:]):
+                v1 = A[:, col1]
+                v2 = A[:, col2]
+                if abs(np.dot(v1, v2)) < tolerance:  # Prefer orthogonal vectors
+                    best_col = col1 if col1 not in column_order else col2
+                    column_order.append(best_col)
+                    break
+
+    # Reorder the columns of the matrix
+    A_reordered = A[:, column_order]
+
+    return A_reordered, column_order
+
+
 def characterize_general_quadratic( surface ): #s surface
     gq_tol = 1e-6
     equivalence_tol = 1e-8
-    a = surface.coefficients['a']
-    b = surface.coefficients['b']
-    c = surface.coefficients['c']
-    d = surface.coefficients['d']
-    e = surface.coefficients['e']
-    f = surface.coefficients['f']
-    g = surface.coefficients['g']
-    h = surface.coefficients['h']
-    j = surface.coefficients['j']
-    k = surface.coefficients['k']
+
+    a = np.round(surface.coefficients['a'], 8)
+    b = np.round(surface.coefficients['b'], 8)
+    c = np.round(surface.coefficients['c'], 8)
+    d = np.round(surface.coefficients['d'], 8)
+    e = np.round(surface.coefficients['e'], 8)
+    f = np.round(surface.coefficients['f'], 8)
+    g = np.round(surface.coefficients['g'], 8)
+    h = np.round(surface.coefficients['h'], 8)
+    j = np.round(surface.coefficients['j'], 8)
+    k = np.round(surface.coefficients['k'], 8)
+
     #coefficient matrix
     Aa = np.asarray([[a, d/2, f/2],
                      [d/2, b, e/2],
@@ -45,6 +96,10 @@ def characterize_general_quadratic( surface ): #s surface
         delta = -1 if det_Ac < 0 else 1
 
     eigenvalues, eigenvectors = np.linalg.eig(Aa)
+
+    eigenvectors, order = _reorder_columns_for_orthogonality(eigenvectors)
+    eigenvalues = eigenvalues[order] # Reorder eigenvalues
+
     signs = np.array([0, 0, 0])
     for i in range(0, 3):
         if eigenvalues[i] > -1 * gq_tol:

src/openmc_cad_adapter/to_cubit_journal.py

@@ -134,26 +134,26 @@ def rotation_to_axis_angle( mat ):
 
                         if abs(theta) <= np.finfo(np.float64).eps:
                             return ( np.array([ 0, 0, 0 ]), 0 )
-                        elif abs( theta - math.pi ) <= np.finfo(np.float64).eps:
-                          # theta is pi (180 degrees) or extremely close to it
-                          # find the column of mat with the largest diagonal
-                          col = 0
-                          if mat[1,1] > mat[col,col]: col = 1
-                          if mat[2,2] > mat[col,col]: col = 2
-
-                          axis = np.array([ 0, 0, 0 ])
-
-                          axis[col] = math.sqrt( (mat[col,col]+1)/2 )
-                          denom = 2*axis[col]
-                          axis[(col+1)%3] = mat[col,(col+1)%3] / denom
-                          axis[(col+2)%3] = mat[col,(col+2)%3] / denom
-                          return ( axis, theta )
+                        if abs( theta - math.pi ) <= np.finfo(np.float64).eps:
+                            # theta is pi (180 degrees) or extremely close to it
+                            # find the column of mat with the largest diagonal
+                            col = 0
+                            if mat[1,1] > mat[col,col]: col = 1
+                            if mat[2,2] > mat[col,col]: col = 2
+
+                            axis = np.array([ 0, 0, 0 ])
+
+                            axis[col] = math.sqrt( (mat[col,col]+1)/2 )
+                            denom = 2*axis[col]
+                            axis[(col+1)%3] = mat[col,(col+1)%3] / denom
+                            axis[(col+2)%3] = mat[col,(col+2)%3] / denom
+                            return ( axis, theta )
                         else:
-                          axis = np.array([ x/r, y/r, z/r ])
-                          return ( axis, theta )
+                            axis = np.array([ x/r, y/r, z/r ])
+                        return ( axis, theta )
                     (r_axis, r_theta ) = rotation_to_axis_angle( rotation_matrix )
                     #compensate for cubits insertion of a negative
-                    r_degs = - math.degrees( r_theta )
+                    r_degs = math.degrees( r_theta )
                     print( r_axis, math.degrees( r_theta ), r_degs )
                     if gq_type == ELLIPSOID : #1
                             r1 = math.sqrt( abs( -K_/A_ ) )