diff --git a/src/openmc_cad_adapter/surfaces.py b/src/openmc_cad_adapter/surfaces.py index 6d7f8ae..a7c6922 100644 --- a/src/openmc_cad_adapter/surfaces.py +++ b/src/openmc_cad_adapter/surfaces.py @@ -54,38 +54,38 @@ def to_cubit_surface_inner(self, ent_type, node, extents, inner_world=None, hex= cmds = [] n = np.array([self.coefficients[k] for k in ('a', 'b', 'c')]) - cd = self.coefficients['d'] - maxv = sys.float_info.min - for i, v in enumerate(n): - if abs(v) > maxv: - maxv = v - - ns = cd * n - - cmds.append( f"create surface rectangle width { 2*extents[0] } zplane") - sur = emit_get_last_id("surface", cmds) - surf = emit_get_last_id("body", cmds) - - n = n/np.linalg.norm(n) - ns = cd * n - zn = np.array([ 0, 0, 1 ]) - n3 = np.cross(n, zn) - dot = np.dot(n, zn) - cmds.append(f"# n3 {n3[0]} {n3[1]} {n3[2]}") - degs = math.degrees(math.acos(np.dot(n, zn))) - y = np.linalg.norm(n3) - x = dot - angle = - math.degrees(math.atan2( y, x )) - if n3[0] != 0 or n3[1] != 0 or n3[2] != 0: - cmds.append(f"Rotate body {{ {surf} }} about 0 0 0 direction {n3[0]} {n3[1]} {n3[2]} Angle {angle}") - cmds.append(f"body {{ { surf } }} move {ns[0]} {ns[1]} {ns[2]}") - cmds.append(f"brick x {extents[0]} y {extents[1]} z {extents[2]}" ) + distance = self.coefficients['d'] / np.linalg.norm(n) + + # Create cutter block larger than the world and rotate/translate it so + # the +z plane of the block is coincident with this general plane + max_extent = np.max(extents) + cmds.append(f"brick x {2*max_extent} y {2*max_extent} z {2*max_extent}" ) ids = emit_get_last_id( ent_type, cmds) - cmds.append(f"section body {{ {ids} }} with surface {{ {sur} }} {self.lreverse(node)}") - cmds.append(f"del surface {{ {sur} }}") + cmds.append(f"body {{ { ids } }} move 0.0 0.0 {-max_extent}") - cmds += self.boundary_condition(ids) - return ids, cmds + nhat = n / np.linalg.norm(n) + rhat = np.array([0.0, 0.0, 1.0]) + angle = math.degrees(math.acos(np.dot(nhat, rhat))) + + if not math.isclose(angle, 0.0, abs_tol=1e-6): + rot_axis = np.cross(rhat, nhat) + rot_axis /= np.linalg.norm(rot_axis) + axis = f"{rot_axis[0]} {rot_axis[1]} {rot_axis[2]}" + cmds.append(f"Rotate body {{ {ids} }} about 0 0 0 direction {axis} Angle {angle}") + + tvec = distance*nhat + cmds.append(f"body {{ { ids } }} move {tvec[0]} {tvec[1]} {tvec[2]}") + + cmds.append(f"brick x {extents[0]} y {extents[1]} z {extents[2]}" ) + wid = emit_get_last_id( ent_type, cmds) + # if positive half space we subtract the cutter block from the world + if node.side != '-': + cmds.append(f"subtract body {{ { ids } }} from body {{ { wid } }}") + # if negative half space we intersect the cutter block with the world + else: + cmds.append(f"intersect body {{ { ids } }} {{ { wid } }}") + + return wid, cmds @classmethod def from_openmc_surface_inner(cls, plane): diff --git a/test/gold/plane.jou b/test/gold/plane.jou index 5329c23..e188cec 100644 --- a/test/gold/plane.jou +++ b/test/gold/plane.jou @@ -5,88 +5,76 @@ graphics pause set journal off set default autosize off #CELL 1 -create surface rectangle width 1000 zplane -#{ id1 = Id("surface") } -#{ id2 = Id("body") } -# n3 0.7071067811865475 -0.7071067811865475 0.0 -Rotate body { id2 } about 0 0 0 direction 0.7071067811865475 -0.7071067811865475 0.0 Angle -90.0 -body { id2 } move -3.5355339059327373 -3.5355339059327373 -0.0 +brick x 1000 y 1000 z 1000 +#{ id1 = Id("body") } +body { id1 } move 0.0 0.0 -500 +Rotate body { id1 } about 0 0 0 direction -0.7071067811865476 0.7071067811865476 0.0 Angle 90.0 +body { id1 } move -2.4999999999999996 -2.4999999999999996 -0.0 brick x 500 y 500 z 500 +#{ id2 = Id("body") } +subtract body { id1 } from body { id2 } +brick x 1000 y 1000 z 1000 #{ id3 = Id("body") } -section body { id3 } with surface { id1 } reverse -del surface { id1 } -create surface rectangle width 1000 zplane -#{ id4 = Id("surface") } -#{ id5 = Id("body") } -# n3 0.7071067811865475 -0.7071067811865475 0.0 -Rotate body { id5 } about 0 0 0 direction 0.7071067811865475 -0.7071067811865475 0.0 Angle -90.0 -body { id5 } move 3.5355339059327373 3.5355339059327373 0.0 +body { id3 } move 0.0 0.0 -500 +Rotate body { id3 } about 0 0 0 direction -0.7071067811865476 0.7071067811865476 0.0 Angle 90.0 +body { id3 } move 2.4999999999999996 2.4999999999999996 0.0 brick x 500 y 500 z 500 +#{ id4 = Id("body") } +intersect body { id3 } { id4 } +#{ id5 = Id("body") } +intersect body { id2 } { id4 } #{ id6 = Id("body") } -section body { id6 } with surface { id4 } -del surface { id4 } -#{ id7 = Id("body") } -intersect body { id3 } { id6 } +#{id7 = ( id5 == id6 ) ? id4 : id6} +brick x 1000 y 1000 z 1000 #{ id8 = Id("body") } -#{id9 = ( id7 == id8 ) ? id6 : id8} -create surface rectangle width 1000 zplane -#{ id10 = Id("surface") } -#{ id11 = Id("body") } -# n3 0.7071067811865475 0.0 0.0 -Rotate body { id11 } about 0 0 0 direction 0.7071067811865475 0.0 0.0 Angle -45.0 -body { id11 } move -0.0 -3.5355339059327373 -3.5355339059327373 +body { id8 } move 0.0 0.0 -500 +Rotate body { id8 } about 0 0 0 direction -1.0 0.0 0.0 Angle 45.00000000000001 +body { id8 } move -0.0 -2.4999999999999996 -2.4999999999999996 brick x 500 y 500 z 500 -#{ id12 = Id("body") } -section body { id12 } with surface { id10 } reverse -del surface { id10 } +#{ id9 = Id("body") } +subtract body { id8 } from body { id9 } +#{ id10 = Id("body") } +intersect body { id7 } { id9 } +#{ id11 = Id("body") } +#{id12 = ( id10 == id11 ) ? id9 : id11} +brick x 1000 y 1000 z 1000 #{ id13 = Id("body") } -intersect body { id9 } { id12 } -#{ id14 = Id("body") } -#{id15 = ( id13 == id14 ) ? id12 : id14} -create surface rectangle width 1000 zplane -#{ id16 = Id("surface") } -#{ id17 = Id("body") } -# n3 0.7071067811865475 0.0 0.0 -Rotate body { id17 } about 0 0 0 direction 0.7071067811865475 0.0 0.0 Angle -45.0 -body { id17 } move 0.0 3.5355339059327373 3.5355339059327373 +body { id13 } move 0.0 0.0 -500 +Rotate body { id13 } about 0 0 0 direction -1.0 0.0 0.0 Angle 45.00000000000001 +body { id13 } move 0.0 2.4999999999999996 2.4999999999999996 brick x 500 y 500 z 500 +#{ id14 = Id("body") } +intersect body { id13 } { id14 } +#{ id15 = Id("body") } +intersect body { id12 } { id14 } +#{ id16 = Id("body") } +#{id17 = ( id15 == id16 ) ? id14 : id16} +brick x 1000 y 1000 z 1000 #{ id18 = Id("body") } -section body { id18 } with surface { id16 } -del surface { id16 } +body { id18 } move 0.0 0.0 -500 +Rotate body { id18 } about 0 0 0 direction 0.0 1.0 0.0 Angle 45.00000000000001 +body { id18 } move -2.4999999999999996 -0.0 -2.4999999999999996 +brick x 500 y 500 z 500 #{ id19 = Id("body") } -intersect body { id15 } { id18 } +subtract body { id18 } from body { id19 } #{ id20 = Id("body") } -#{id21 = ( id19 == id20 ) ? id18 : id20} -create surface rectangle width 1000 zplane -#{ id22 = Id("surface") } +intersect body { id17 } { id19 } +#{ id21 = Id("body") } +#{id22 = ( id20 == id21 ) ? id19 : id21} +brick x 1000 y 1000 z 1000 #{ id23 = Id("body") } -# n3 0.0 -0.7071067811865475 0.0 -Rotate body { id23 } about 0 0 0 direction 0.0 -0.7071067811865475 0.0 Angle -45.0 -body { id23 } move -3.5355339059327373 -0.0 -3.5355339059327373 +body { id23 } move 0.0 0.0 -500 +Rotate body { id23 } about 0 0 0 direction 0.0 1.0 0.0 Angle 45.00000000000001 +body { id23 } move 2.4999999999999996 0.0 2.4999999999999996 brick x 500 y 500 z 500 #{ id24 = Id("body") } -section body { id24 } with surface { id22 } reverse -del surface { id22 } +intersect body { id23 } { id24 } #{ id25 = Id("body") } -intersect body { id21 } { id24 } +intersect body { id22 } { id24 } #{ id26 = Id("body") } #{id27 = ( id25 == id26 ) ? id24 : id26} -create surface rectangle width 1000 zplane -#{ id28 = Id("surface") } -#{ id29 = Id("body") } -# n3 0.0 -0.7071067811865475 0.0 -Rotate body { id29 } about 0 0 0 direction 0.0 -0.7071067811865475 0.0 Angle -45.0 -body { id29 } move 3.5355339059327373 0.0 3.5355339059327373 -brick x 500 y 500 z 500 -#{ id30 = Id("body") } -section body { id30 } with surface { id28 } -del surface { id28 } -#{ id31 = Id("body") } -intersect body { id27 } { id30 } -#{ id32 = Id("body") } -#{id33 = ( id31 == id32 ) ? id30 : id32} -body { id33 } name "Cell_1" -group "mat:void" add body { id33 } +body { id27 } name "Cell_1" +group "mat:void" add body { id27 } graphics flush set default autosize on zoom reset