diff --git a/src/openmc_cad_adapter/gqs.py b/src/openmc_cad_adapter/gqs.py index 7b0a510..6eeb775 100644 --- a/src/openmc_cad_adapter/gqs.py +++ b/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"] \ No newline at end of file + "PARABOLIC_CYLINDER"] diff --git a/test/test_local.py b/test/test_local.py index 31479f7..095fddd 100644 --- a/test/test_local.py +++ b/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') diff --git a/test/test_quadric.py b/test/test_quadric.py new file mode 100644 index 0000000..9b58717 --- /dev/null +++ b/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 \ No newline at end of file