From 08e68dbb3734822f2eb9ed2cf8958dd40f391b59 Mon Sep 17 00:00:00 2001 From: Stefano Segantin Date: Fri, 22 Nov 2024 10:23:02 -0500 Subject: [PATCH 1/2] first commit --- src/openmc_cad_adapter/to_cubit_journal.py | 372 ++++++++++----------- 1 file changed, 186 insertions(+), 186 deletions(-) diff --git a/src/openmc_cad_adapter/to_cubit_journal.py b/src/openmc_cad_adapter/to_cubit_journal.py index eff69ad..5aac483 100644 --- a/src/openmc_cad_adapter/to_cubit_journal.py +++ b/src/openmc_cad_adapter/to_cubit_journal.py @@ -108,200 +108,200 @@ def surface_to_cubit_journal(node, w, indent = 0, inner_world = None, hex = Fals def ind(): return ' ' * (2*indent) if isinstance(node, Halfspace): - seen.add( node.surface ) - surface = node.surface - - nonlocal cmds - - def reverse(): - return "reverse" if node.side == '-' else "" - - if cad_surface := _CAD_SURFACE_DICTIONARY.get(surface._type): - cad_surface = cad_surface.from_openmc_surface(surface) - ids, cad_cmds = cad_surface.to_cubit_surface(ent_type, node, w, inner_world, hex) - cmds += cad_cmds - return ids - elif surface._type == "quadric": - (gq_type, A_, B_, C_, K_, translation, rotation_matrix) = characterize_general_quadratic(surface) - - def rotation_to_axis_angle( mat ): - x = mat[2, 1]-mat[1, 2] - y = mat[0, 2]-mat[2, 0] - z = mat[1, 0]-mat[0, 1] - r = math.hypot( x, math.hypot( y,z )) - t = mat[0,0] + mat[1,1] + mat[2,2] - theta = math.atan2(r,t-1) - - 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 ) - else: - 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 ) - print( r_axis, math.degrees( r_theta ), r_degs ) - if gq_type == ELLIPSOID : #1 - r1 = math.sqrt( abs( -K_/A_ ) ) - r2 = math.sqrt( abs( -K_/B_ ) ) - r3 = math.sqrt( abs( -K_/C_ ) ) - cmds.append( f"sphere radius 1") - ids = emit_get_last_id( ent_type , cmds) - cmds.append( f"body {{ { ids } }} scale x { r1 } y { r2 } z { r3 }") - move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) - elif gq_type == ELLIPTIC_CYLINDER : #7 - if A_ == 0: - print( "X", gq_type, A_, B_, C_, K_, r_axis, r_degs ) - h = inner_world[0] if inner_world else w[0] - r1 = math.sqrt( abs( K_/C_ ) ) - r2 = math.sqrt( abs( K_/B_ ) ) - cmds.append( f"cylinder height {h} Major Radius {r1} Minor Radius {r2}") - ids = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ { ids } }} about y angle 90") - if node.side != '-': - wid = 0 - if inner_world: - if hex: - cmds.append( f"create prism height {inner_world[2]} sides 6 radius { inner_world[0] / 2 } " ) - wid = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ {wid} }} about z angle 30" ) - cmds.append( f"rotate body {{ {wid} }} about y angle 90") - else: - cmds.append( f"brick x {inner_world[0]} y {inner_world[1]} z {inner_world[2]}" ) - wid = emit_get_last_id( ent_type , cmds) + seen.add( node.surface ) + surface = node.surface + + nonlocal cmds + + def reverse(): + return "reverse" if node.side == '-' else "" + + if cad_surface := _CAD_SURFACE_DICTIONARY.get(surface._type): + cad_surface = cad_surface.from_openmc_surface(surface) + ids, cad_cmds = cad_surface.to_cubit_surface(ent_type, node, w, inner_world, hex) + cmds += cad_cmds + return ids + elif surface._type == "quadric": + (gq_type, A_, B_, C_, K_, translation, rotation_matrix) = characterize_general_quadratic(surface) + + def rotation_to_axis_angle( mat ): + x = mat[2, 1]-mat[1, 2] + y = mat[0, 2]-mat[2, 0] + z = mat[1, 0]-mat[0, 1] + r = math.hypot( x, math.hypot( y,z )) + t = mat[0,0] + mat[1,1] + mat[2,2] + theta = math.atan2(r,t-1) + + 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 ) + else: + 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 ) + print( r_axis, math.degrees( r_theta ), r_degs ) + if gq_type == ELLIPSOID : #1 + r1 = math.sqrt( abs( -K_/A_ ) ) + r2 = math.sqrt( abs( -K_/B_ ) ) + r3 = math.sqrt( abs( -K_/C_ ) ) + cmds.append( f"sphere radius 1") + ids = emit_get_last_id( ent_type , cmds) + cmds.append( f"body {{ { ids } }} scale x { r1 } y { r2 } z { r3 }") + move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) + elif gq_type == ELLIPTIC_CYLINDER : #7 + if A_ == 0: + print( "X", gq_type, A_, B_, C_, K_, r_axis, r_degs ) + h = inner_world[0] if inner_world else w[0] + r1 = math.sqrt( abs( K_/C_ ) ) + r2 = math.sqrt( abs( K_/B_ ) ) + cmds.append( f"cylinder height {h} Major Radius {r1} Minor Radius {r2}") + ids = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ { ids } }} about y angle 90") + if node.side != '-': + wid = 0 + if inner_world: + if hex: + cmds.append( f"create prism height {inner_world[2]} sides 6 radius { inner_world[0] / 2 } " ) + wid = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ {wid} }} about z angle 30" ) + cmds.append( f"rotate body {{ {wid} }} about y angle 90") else: - cmds.append( f"brick x {w[0]} y {w[1]} z {w[2]}" ) + cmds.append( f"brick x {inner_world[0]} y {inner_world[1]} z {inner_world[2]}" ) wid = emit_get_last_id( ent_type , cmds) - cmds.append( f"subtract body {{ { ids } }} from body {{ { wid } }}" ) - cmds.append( f"Rotate body {{ {wid } }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( wid, translation[0,0], translation[1,0], translation[2,0], cmds) - return wid - cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) - return ids - if B_ == 0: - print( "Y", gq_type, A_, B_, C_, K_ ) - h = inner_world[1] if inner_world else w[1] - r1 = math.sqrt( abs( K_/A_ ) ) - r2 = math.sqrt( abs( K_/C_ ) ) - cmds.append( f"cylinder height {h} Major Radius {r1} Minor Radius {r2}") - ids = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ { ids } }} about x angle 90") - if node.side != '-': - wid = 0 - if inner_world: - if hex: - cmds.append( f"create prism height {inner_world[2]} sides 6 radius { inner_world[0] / 2 } " ) - wid = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ {wid} }} about z angle 30" ) - cmds.append( f"rotate body {{ {wid} }} about y angle 90") - else: - cmds.append( f"brick x {inner_world[0]} y {inner_world[1]} z {inner_world[2]}" ) - wid = emit_get_last_id( ent_type , cmds) + else: + cmds.append( f"brick x {w[0]} y {w[1]} z {w[2]}" ) + wid = emit_get_last_id( ent_type , cmds) + cmds.append( f"subtract body {{ { ids } }} from body {{ { wid } }}" ) + cmds.append( f"Rotate body {{ {wid } }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( wid, translation[0,0], translation[1,0], translation[2,0], cmds) + return wid + cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) + return ids + if B_ == 0: + print( "Y", gq_type, A_, B_, C_, K_ ) + h = inner_world[1] if inner_world else w[1] + r1 = math.sqrt( abs( K_/A_ ) ) + r2 = math.sqrt( abs( K_/C_ ) ) + cmds.append( f"cylinder height {h} Major Radius {r1} Minor Radius {r2}") + ids = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ { ids } }} about x angle 90") + if node.side != '-': + wid = 0 + if inner_world: + if hex: + cmds.append( f"create prism height {inner_world[2]} sides 6 radius { inner_world[0] / 2 } " ) + wid = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ {wid} }} about z angle 30" ) + cmds.append( f"rotate body {{ {wid} }} about y angle 90") else: - cmds.append( f"brick x {w[0]} y {w[1]} z {w[2]}" ) + cmds.append( f"brick x {inner_world[0]} y {inner_world[1]} z {inner_world[2]}" ) wid = emit_get_last_id( ent_type , cmds) - cmds.append( f"subtract body {{ { ids } }} from body {{ { wid } }}" ) - cmds.append( f"Rotate body {{ {wid } }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( wid, translation[0,0], translation[1,0], translation[2,0], cmds) - return wid - cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) - return ids - if C_ == 0: - print( "Z", gq_type, A_, B_, C_, K_ ) - h = inner_world[2] if inner_world else w[2] - r1 = math.sqrt( abs( K_/A_ ) ) - r2 = math.sqrt( abs( K_/B_ ) ) - cmds.append( f"cylinder height {h} Major Radius {r1} Minor Radius {r2}") - ids = emit_get_last_id( ent_type , cmds) - if node.side != '-': - wid = 0 - if inner_world: - if hex: - cmds.append( f"create prism height {inner_world[2]} sides 6 radius { inner_world[0] / 2 } " ) - wid = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ {wid} }} about z angle 30" ) - cmds.append( f"rotate body {{ {wid} }} about y angle 90") - else: - cmds.append( f"brick x {inner_world[0]} y {inner_world[1]} z {inner_world[2]}" ) - wid = emit_get_last_id( ent_type , cmds) + else: + cmds.append( f"brick x {w[0]} y {w[1]} z {w[2]}" ) + wid = emit_get_last_id( ent_type , cmds) + cmds.append( f"subtract body {{ { ids } }} from body {{ { wid } }}" ) + cmds.append( f"Rotate body {{ {wid } }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( wid, translation[0,0], translation[1,0], translation[2,0], cmds) + return wid + cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) + return ids + if C_ == 0: + print( "Z", gq_type, A_, B_, C_, K_ ) + h = inner_world[2] if inner_world else w[2] + r1 = math.sqrt( abs( K_/A_ ) ) + r2 = math.sqrt( abs( K_/B_ ) ) + cmds.append( f"cylinder height {h} Major Radius {r1} Minor Radius {r2}") + ids = emit_get_last_id( ent_type , cmds) + if node.side != '-': + wid = 0 + if inner_world: + if hex: + cmds.append( f"create prism height {inner_world[2]} sides 6 radius { inner_world[0] / 2 } " ) + wid = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ {wid} }} about z angle 30" ) + cmds.append( f"rotate body {{ {wid} }} about y angle 90") else: - cmds.append( f"brick x {w[0]} y {w[1]} z {w[2]}" ) + cmds.append( f"brick x {inner_world[0]} y {inner_world[1]} z {inner_world[2]}" ) wid = emit_get_last_id( ent_type , cmds) - cmds.append( f"subtract body {{ { ids } }} from body {{ { wid } }}" ) - cmds.append( f"Rotate body {{ {wid } }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( wid, translation[0,0], translation[1,0], translation[2,0], cmds) - return wid - cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) - return ids - elif gq_type == ELLIPTIC_CONE : #3 - if A_ == 0: - h = inner_world[0] if inner_world else w[0] - minor = math.sqrt( abs( -A_/C_ ) ) - major = math.sqrt( abs( -A_/B_ ) ) - rot_angle = - 90 - rot_axis = 1 - cmds.append( f"create frustum height {h} Major Radius {major} Minor Radius {minor} top 0") - ids = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ { ids } }} about y angle -90") - cmds.append( f"copy body {{ { ids } }}") - mirror = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ { mirror } }} about 0 0 0 angle 180") - cmds.append( f"unit body {{ { ids } }} {{ { mirror } }}") - cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) - return ids - if B_ == 0: - h = inner_world[1] if inner_world else w[1] - minor = math.sqrt( abs( -B_/A_ ) ) - major = math.sqrt( abs( -B_/C_ ) ) - rot_angle = 90 - rot_axis = 0 - cmds.append( f"create frustum height {h} Major Radius {major} Minor Radius {minor} top 0") - ids = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ { ids } }} about x angle 90") - cmds.append( f"copy body {{ { ids } }}") - mirror = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ { mirror } }} about 0 0 0 angle 180") - cmds.append( f"unit body {{ { ids } }} {{ { mirror } }}") - cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) - return ids - if C_ == 0: - h = inner_world[2] if inner_world else w[2] - minor = math.sqrt( abs( -C_/A_ ) ) - major = math.sqrt( abs( -C_/B_ ) ) - rot_angle = 180 - rot_axis = 0 - cmds.append( f"create frustum height {h} Major Radius {major} Minor Radius {minor} top 0") - ids = emit_get_last_id( ent_type , cmds) - cmds.append( f"copy body {{ { ids } }}") - mirror = emit_get_last_id( ent_type , cmds) - cmds.append( f"rotate body {{ { mirror } }} about 0 0 0 angle 180") - cmds.append( f"unit body {{ { ids } }} {{ { mirror } }}") - cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") - move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) - return ids - else: - raise NotImplementedError(f"{surface.type} not implemented") + else: + cmds.append( f"brick x {w[0]} y {w[1]} z {w[2]}" ) + wid = emit_get_last_id( ent_type , cmds) + cmds.append( f"subtract body {{ { ids } }} from body {{ { wid } }}" ) + cmds.append( f"Rotate body {{ {wid } }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( wid, translation[0,0], translation[1,0], translation[2,0], cmds) + return wid + cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) + return ids + elif gq_type == ELLIPTIC_CONE : #3 + if A_ == 0: + h = inner_world[0] if inner_world else w[0] + minor = math.sqrt( abs( -A_/C_ ) ) + major = math.sqrt( abs( -A_/B_ ) ) + rot_angle = - 90 + rot_axis = 1 + cmds.append( f"create frustum height {h} Major Radius {major} Minor Radius {minor} top 0") + ids = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ { ids } }} about y angle -90") + cmds.append( f"copy body {{ { ids } }}") + mirror = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ { mirror } }} about 0 0 0 angle 180") + cmds.append( f"unit body {{ { ids } }} {{ { mirror } }}") + cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) + return ids + if B_ == 0: + h = inner_world[1] if inner_world else w[1] + minor = math.sqrt( abs( -B_/A_ ) ) + major = math.sqrt( abs( -B_/C_ ) ) + rot_angle = 90 + rot_axis = 0 + cmds.append( f"create frustum height {h} Major Radius {major} Minor Radius {minor} top 0") + ids = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ { ids } }} about x angle 90") + cmds.append( f"copy body {{ { ids } }}") + mirror = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ { mirror } }} about 0 0 0 angle 180") + cmds.append( f"unit body {{ { ids } }} {{ { mirror } }}") + cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) + return ids + if C_ == 0: + h = inner_world[2] if inner_world else w[2] + minor = math.sqrt( abs( -C_/A_ ) ) + major = math.sqrt( abs( -C_/B_ ) ) + rot_angle = 180 + rot_axis = 0 + cmds.append( f"create frustum height {h} Major Radius {major} Minor Radius {minor} top 0") + ids = emit_get_last_id( ent_type , cmds) + cmds.append( f"copy body {{ { ids } }}") + mirror = emit_get_last_id( ent_type , cmds) + cmds.append( f"rotate body {{ { mirror } }} about 0 0 0 angle 180") + cmds.append( f"unit body {{ { ids } }} {{ { mirror } }}") + cmds.append( f"Rotate body {{ {ids} }} about 0 0 0 direction {r_axis[0]} {r_axis[1]} {r_axis[2]} Angle {r_degs}") + move( ids, translation[0,0], translation[1,0], translation[2,0], cmds) + return ids else: raise NotImplementedError(f"{surface.type} not implemented") + else: + raise NotImplementedError(f"{surface.type} not implemented") elif isinstance(node, Complement): print( "Complement:" ) id = surface_to_cubit_journal(node.node, w, indent + 1, inner_world, ent_type = ent_type ) From d8e770d0985cb8be3573ad9d9d6408a27fb0d45b Mon Sep 17 00:00:00 2001 From: Stefano Segantin Date: Mon, 25 Nov 2024 17:03:38 -0500 Subject: [PATCH 2/2] fixed eigenvalue sign in gqs --- src/openmc_cad_adapter/gqs.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/openmc_cad_adapter/gqs.py b/src/openmc_cad_adapter/gqs.py index 380341a..5be38d2 100644 --- a/src/openmc_cad_adapter/gqs.py +++ b/src/openmc_cad_adapter/gqs.py @@ -42,15 +42,15 @@ 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 = np.sign(det_Ac) eigenvalues, eigenvectors = np.linalg.eig(Aa) signs = np.array([0, 0, 0]) for i in range(0, 3): if eigenvalues[i] > -1 * gq_tol: - signs[i] = 1 - else: signs[i] = -1 + else: + signs[i] = 1 S = 1 if np.abs( signs.sum() ) == 3 else -1