Merge pull request #13 from pshriwise/gq-testing

Quadric classification testing

src/openmc_cad_adapter/gqs.py

@@ -1,6 +1,7 @@
 import numpy as np
 from numpy.linalg import matrix_rank
 
+UNKNOWN_QUADRIC = -1
 ELLIPSOID = 1
 ONE_SHEET_HYPERBOLOID = 2
 TWO_SHEET_HYPERBOLOID = 3
@@ -43,7 +44,8 @@ def characterize_general_quadratic( surface ): #s surface
     if np.abs( det_Ac ) < gq_tol:
         delta = 0
     else:
-        delta = -1 if det_Ac < 0 else -1
+        delta = -1 if det_Ac < 0 else 1
+
     eigen_results = np.linalg.eig(Aa)
     signs = np.array([ 0, 0, 0 ])
     for i in range( 0, 3 ):
@@ -65,7 +67,8 @@ def characterize_general_quadratic( surface ): #s surface
     dz = C[2]
 
     #Update the constant using the resulting translation
-    K_ = k + g/2*dx + h/2*dy + j/2*dz
+    K_ = k + (g/2)*dx + (h/2)*dy + (j/2)*dz
+    K_ = K_[0,0]
 
     if rank_Aa == 2 and rank_Ac == 3 and S == 1:
         delta = -1 if K_ * signs[0] else 1
@@ -74,30 +77,40 @@ def characterize_general_quadratic( surface ): #s surface
 
 
     def find_type( rAa, rAc, delta, S, D ):
+        quadric_type = UNKNOWN_QUADRIC
         if 3 == rAa and 4 == rAc and -1 == delta and 1 == S:
-            return ELLIPSOID
+            quadric_type = ELLIPSOID
         elif 3 == rAa and 4 == rAc and 1 == delta and -1 == S:
-            return ONE_SHEET_HYPERBOLOID
+            quadric_type = ONE_SHEET_HYPERBOLOID
         elif 3 == rAa and 4 == rAc and -1 == delta and -1 == S:
-            return TWO_SHEET_HYPERBOLOID
+            quadric_type = TWO_SHEET_HYPERBOLOID
         elif 3 == rAa and 3 == rAc and 0 == delta and -1 == S:
-            return ELLIPTIC_CONE
+            quadric_type = ELLIPTIC_CONE
         elif 2 == rAa and 4 == rAc and -1 == delta and 1 == S:
-            return ELLIPTIC_PARABOLOID
+            quadric_type = ELLIPTIC_PARABOLOID
         elif 2 == rAa and 4 == rAc and 1 == delta and -1 == S:
-            return HYPERBOLIC_PARABOLOID
+            quadric_type = HYPERBOLIC_PARABOLOID
         elif 2 == rAa and 3 == rAc and -1 == delta and 1 == S:
-            return ELLIPTIC_CYLINDER
+            quadric_type = ELLIPTIC_CYLINDER
         elif 2 == rAa and 3 == rAc and 0 == delta and -1 == S:
-            return HYPERBOLIC_CYLINDER
-        elif 2 == rAa and 3 == rAc and 0 == delta and 1 == S:
-            return PARABOLIC_CYLINDER
-        elif 2 == rAa and 3 == rAc and 1 == S and D != 0 :
-            return find_type( rAa, rAc, D, S, 0 )
+            quadric_type = HYPERBOLIC_CYLINDER
+        elif 1 == rAa and 3 == rAc and 0 == delta and 1 == S:
+            quadric_type = PARABOLIC_CYLINDER
         else:
-            raise "UNKNOWN QUADRATIC"
+            quadric_type = UNKNOWN_QUADRIC
+
+        # special case, replace delta with D
+        if 2 == rAa and 3 == rAc and 1 == S and D != 0 :
+            quadric_type = find_type( rAa, rAc, D, S, 0 )
+
 
-    gq_type = find_type( rank_Aa, rank_Ac, delta, S, D )
+        if quadric_type == UNKNOWN_QUADRIC:
+            msg = f'UNKNOWN QUADRIC: rAa={rAa}, rAc={rAc}, delta={delta}, S={S}, D={D}'
+            raise ValueError(msg)
+
+        return quadric_type
+
+    gq_type = find_type(rank_Aa, rank_Ac, delta, S, D)
 
     #set the translation
     translation = C
@@ -114,16 +127,19 @@ def find_type( rAa, rAc, delta, S, D ):
     C_ = eigenvalues[2];
     D_ = 0; E_ = 0; F_ = 0;
     G_ = 0; H_ = 0; J_ = 0;
+
+    # alter type and coefficients for special cases
+    # where coefficients are near-zero
     if gq_type == ONE_SHEET_HYPERBOLOID:
-        if abs( K_) < equivalence_tol:
+        if abs(K_) < equivalence_tol:
            K_ = 0
            gq_type = ELLIPTIC_CONE
     if gq_type == TWO_SHEET_HYPERBOLOID:
-        if abs( K_) < equivalence_tol:
+        if abs(K_) < equivalence_tol:
            K_ = 0
            gq_type = ELLIPTIC_CONE
     if gq_type == ELLIPSOID:
-        if abs( A_) < equivalence_tol:
+        if abs(A_) < equivalence_tol:
            A_ = 0
            gq_type = ELLIPTIC_CYLINDER
         elif abs( B_) < equivalence_tol:
@@ -145,4 +161,4 @@ def find_type( rAa, rAc, delta, S, D ):
            "HYPERBOLIC_PARABOLOID",
            "ELLIPTIC_CYLINDER",
            "HYPERBOLIC_CYLINDER",
-           "PARABOLIC_CYLINDER"]
+           "PARABOLIC_CYLINDER"]

test/test_local.py

@@ -124,8 +124,8 @@ def test_general_cone(request, run_in_tmpdir):
         to_cubit_journal(g, world=(500, 500, 500), filename='cone.jou')
 
 @reset_openmc_ids
-def test_gq_ellipsoid(request):
-    ellipsoid = openmc.Quadric(1, 2, 3, k=1)
+def test_gq_ellipsoid(request, run_in_tmpdir):
+    ellipsoid = openmc.Quadric(1, 2, 3, k=-1)
     g = openmc.Geometry([openmc.Cell(region=-ellipsoid)])
     to_cubit_journal(g, world=(500, 500, 500), filename='ellipsoid.jou')
     diff_gold_file('ellipsoid.jou')

test/test_quadric.py

@@ -0,0 +1,106 @@
+
+import pytest
+
+import openmc
+
+from openmc_cad_adapter import to_cubit_journal
+from openmc_cad_adapter.gqs import *
+
+
+def test_ellipsoid_classification():
+    # ELLIPSOID
+    testEllip = openmc.Quadric(1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testEllip)
+    assert quadric_type == ELLIPSOID
+
+
+def test_one_sheet_hyperboloid_classification():
+    # ONE_SHEET_HYPERBOLOID
+    testOneSheet = openmc.Quadric(1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testOneSheet)
+    assert quadric_type == ONE_SHEET_HYPERBOLOID
+
+
+def test_two_sheet_hyperboloid_classification():
+    # TWO_SHEET_HYPERBOLOID
+    testTwoSheet = openmc.Quadric(-1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testTwoSheet)
+    assert quadric_type == TWO_SHEET_HYPERBOLOID
+
+
+def test_elliptic_cone_classification():
+    # ELLIPTIC_CONE
+    testEllCone = openmc.Quadric(1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testEllCone)
+    assert quadric_type == ELLIPTIC_CONE
+
+
+def test_elliptic_paraboloid_classification():
+    # ELLIPTIC_PARABOLOID
+    testEllPara = openmc.Quadric(1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testEllPara)
+    assert quadric_type == ELLIPTIC_PARABOLOID
+
+
+def test_hyperbolic_paraboloid_classification():
+    # HYPERBOLIC_PARABOLOID
+    testHypPara = openmc.Quadric(1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testHypPara)
+    assert quadric_type == HYPERBOLIC_PARABOLOID
+
+
+def test_elliptic_cyl_classification():
+    # ELLIPTIC_CYL
+    testEllCyl = openmc.Quadric(1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testEllCyl)
+    assert quadric_type == ELLIPTIC_CYLINDER
+
+
+def test_hyperbolic_cyl_classification():
+    # HYPERBOLIC_CYL
+    testHypCyl = openmc.Quadric(1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testHypCyl)
+    assert quadric_type == HYPERBOLIC_CYLINDER
+
+
+def test_parabolic_cyl_classification():
+    # PARABOLIC_CYL
+    testParaCyl = openmc.Quadric(1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testParaCyl)
+    assert quadric_type == PARABOLIC_CYLINDER
+
+
+# Transformation Tests
+def test_ellipsoid_classification():
+    # ELLIPSOID
+    testRotEllip = openmc.Quadric(103, 125, 66, -48, -12, -60, 0, 0, 0, -294)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testRotEllip)
+    assert quadric_type == ELLIPSOID
+
+
+def test_elliptic_cone_classification():
+    # ELLIPTIC_CONE
+    testRotCone = openmc.Quadric(3, 3, -1, 2,  0, 0, 0, 0, 0, 0)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testRotCone)
+    assert quadric_type == ELLIPTIC_CONE
+
+
+def test_elliptic_paraboloid_classification():
+    # ELLIPTIC_PARABOLOID
+    testRotEllParab = openmc.Quadric(1, 3, 1, 2, 2, 2, -2, 4, 2, 12)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testRotEllParab)
+    assert quadric_type == ELLIPTIC_PARABOLOID
+
+
+def test_elliptic_cylinder_classification():
+    # ELLIPTIC_CYLINDER
+    testRotEllCyl = openmc.Quadric(5, 2, 5, -4, -4, -2, 6, -12, 18, -3)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testRotEllCyl)
+    assert quadric_type == ELLIPTIC_CYLINDER
+
+
+def test_parabolic_cylinder_classification():
+    # PARABOLIC CYLINDER
+    testRotParaCyl = openmc.Quadric(9, 36, 4, -36, -24, 12, -16, -24, -48, 56)
+    quadric_type, A, B, C, K, _, _ = characterize_general_quadratic(testRotParaCyl)
+    assert quadric_type == PARABOLIC_CYLINDER