Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

General Plane fix #17

Merged
merged 3 commits into from
Nov 22, 2024
Merged

General Plane fix #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4be1a51
fix for general planes using volume intersection/subtration instead
eepeterson Nov 21, 2024
759a811
update plane gold file
eepeterson Nov 21, 2024
4bea4db
Removing commented line
pshriwise Nov 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 30 additions & 30 deletions src/openmc_cad_adapter/surfaces.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
118 changes: 53 additions & 65 deletions test/gold/plane.jou
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading