polygon example

[shimwell]

import openmc
import matplotlib.pyplot as plt

height = 20
width = 30
thickness = 2
center_x=100
center_z=0

surface_1 = openmc.model.Polygon(
    points= [
        (10,10),
        (7,13),
        (5,10),
        (7,5),
    ],
    basis='rz'
)

# offsetting is very useful when the geometry is not on the central axis
surface_2 = surface_1.offset(1)

region_1 = -surface_1
region_2 = +surface_1 & -surface_2

cell_1 = openmc.Cell(region=region_1)
cell_2 = openmc.Cell(region=region_2)

universe = openmc.Universe(cells=[cell_1, cell_2])
print(universe.bounding_box)
# set width and origin as bounding box can't be automatically found for polyline
universe.plot(width=(25,20),basis='xz', origin=(0,0,10))

plt.show()

[shimwell]

this is a better example from the test suit

import openmc
import numpy as np
import matplotlib.pyplot as plt
# define a 5 pointed star centered on 1, 1
star = np.array(
    [
        [1.0, 2.0],
        [0.70610737, 1.4045085],
        [0.04894348, 1.30901699],
        [0.52447174, 0.8454915],
        [0.41221475, 0.19098301],
        [1.0, 0.5],
        [1.58778525, 0.19098301],
        # [1.47552826, 0.8454915],
        [1.95105652, 1.30901699],
        [1.29389263, 1.4045085],
        [1.0, 2.0],
    ]
)

surf_poly = openmc.model.Polygon(star, basis="yz")
surf_poly2 = surf_poly.offset(0.5)

region_poly = -surf_poly
region_poly2 = -surf_poly2 & +surf_poly

cell_poly = openmc.Cell(region=region_poly)
# cell_poly.rotation=[45,0,0]
# cell_poly2 = openmc.Cell(region=region_poly2)

my_geometry = openmc.Geometry([cell_poly])
my_geometry.root_universe.plot(basis='yz')

plt.show()
# my_geometry.

[shimwell]

import openmc
import numpy as np
import matplotlib.pyplot as plt
# define a 5 pointed star centered on 1, 1
star = np.array(
    [
        [1.0, 2.0],
        [0.70610737, 1.4045085],
        [0.04894348, 1.30901699],
        [0.52447174, 0.8454915],
        [0.41221475, 0.19098301],
        [1.0, 0.5],
        [1.58778525, 0.19098301],
        [1.47552826, 0.8454915],
        [1.95105652, 1.30901699],
        [1.29389263, 1.4045085],
        [1.0, 2.0],
    ]
)

surf_poly = openmc.model.Polygon(star, basis="rz")
surf_poly.boundary_type='vacuum'
# surf_poly.rotate(5)

region_poly = surf_poly.region

cell_poly = openmc.Cell(region=region_poly)

my_geometry = openmc.Geometry([cell_poly])

my_geometry.root_universe.plot(basis='yz')

# plt.show()
my_geometry.root_universe.plot(basis='xy', origin=(0,0,1))

# plt.show()

material_1 = openmc.Material() 
material_1.add_element('Li', 1, percent_type='ao')
material_1.set_density('g/cm3', 0.5)

my_materials = openmc.Materials([material_1])

my_settings = openmc.Settings()
my_settings.batches = 100
my_settings.inactive = 0
my_settings.particles = 5000
my_settings.particle = "neutron"
my_settings.run_mode = 'fixed source'

my_source = openmc.Source()
my_source.space = openmc.stats.Point((1, 1, 1))
my_source.angle = openmc.stats.Isotropic()
my_source.energy = openmc.stats.Discrete([14e6], [1])
my_settings.source = my_source

mesh = openmc.RegularMesh()
mesh.lower_left = (-max(surf_poly.points[:,1]), -max(surf_poly.points[:,1]), min(surf_poly.points[:,0])) 
mesh.upper_right = (max(surf_poly.points[:,1]), max(surf_poly.points[:,1]), max(surf_poly.points[:,0])) 
mesh.dimension=[100, 100, 100] # only 1 cell in the Y dimension
print(mesh)

mesh_filter = openmc.MeshFilter(mesh)
mesh_tally = openmc.Tally(name='tallies_on_mesh')
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = ['flux']
my_tallies = openmc.Tallies([mesh_tally])

model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies)

output_filename = model.run()

results = openmc.StatePoint(output_filename)

my_tally = results.get_tally(scores=['flux'])

plot = my_tally.plot_mesh_tally()
plot.show()

# my_slice = my_tally.get_slice(scores=['flux'])

# my_slice.mean.shape = (mesh.dimension[0], mesh.dimension[1], mesh.dimension[2])

# plt.imshow(my_slice.mean[:, :, 50])
# plt.show()

[shimwell]

example with transport

import openmc
import matplotlib.pyplot as plt


surface_1 = openmc.model.Polygon(
    points= [
        (10,10), # left
        (7,13), # top
        (5,10), # right
        (7,5), # bottom
    ],
    basis='rz'
)
surface_2 = openmc.model.Polygon(
    points= [
        (10,10), # right
        (7,13), # top
        (15,10),
        (7,5),
    ],
    basis='rz'
)

material_1 = openmc.Material()
material_1.add_nuclide('Li7',1)
material_1.set_density('g/cm3',1)

materials = openmc.Materials([material_1])

surface_vac = openmc.Sphere(r=30, boundary_type='vacuum')

region_1 = -surface_1
region_2 = -surface_2
region_void = +surface_1 & +surface_2 & -surface_vac

cell_1 = openmc.Cell(region=region_1)
cell_2 = openmc.Cell(region=region_2)
cell_3 = openmc.Cell(region=region_void)

geometry = openmc.Geometry([cell_1, cell_2, cell_3])
print(geometry.bounding_box)
# set width and origin as bounding box can't be automatically found for polyline
geometry.plot(width=(35,20),basis='xz', origin=(0,0,10), color_by='cell')

# plt.show()

source = openmc.IndependentSource()

settings = openmc.Settings()
settings.batches = 5
settings.particles = 1000
settings.source = source
settings.run_mode = 'fixed source'

model = openmc.Model(geometry, materials, settings)
model.run()

[shimwell]

varible offset plasma

line_coords = [
    [600.0, 0.0],
    [537.0570808585081, 176.33557568774194],
    [418.8892073777862, 285.31695488854604],
    [338.58025161128063, 285.3169548885461],
    [306.924607710337, 176.33557568774197],
    [300.0, 3.6739403974420595e-14],
    [306.924607710337, -176.3355756877419],
    [338.5802516112806, -285.31695488854604],
    [418.8892073777862, -285.3169548885461],
    [537.0570808585081, -176.335575687742],
    [600.0, 0.0]
]

# line_coords = [(0.1, 0.1), (1, 1), (0.1, 2), (2, 2), (3, 1), (1, 0.1)]

# Function to iterate through the list with a sliding window
def sliding_window(lst, window_size, step_size):
    for i in range(0, len(lst) - window_size + 1, step_size):
        yield lst[i:i + window_size]

def count_sliding_windows(lst, window_size, step_size):
    if window_size > len(lst):
        return 0
    return (len(lst) - window_size) // step_size + 1

# Define window size and step size
window_size = 3
step_size = 1

num_windows = count_sliding_windows(line_coords, window_size, step_size)
# offsets = [0.5, 0.1, 0.2,0.3]
import numpy as np
offsets = np.array([5,7,9,13,14,13,9,7,5])

all_left_hand_sides = []
# Iterate and print each window
for coords, offset in zip(sliding_window(line_coords, window_size, step_size), offsets):
    print(coords, offset)
    line = shapely.LineString(coords)
    left_hand_side = line.buffer(offset, single_sided=True, quad_segs=4, join_style='round')
    all_coordinates = shapely.geometry.mapping(left_hand_side)['coordinates'][0]

    all_left_hand_sides.append(left_hand_side)

    x_values = [coord[0] for coord in all_coordinates]
    y_values = [coord[1] for coord in all_coordinates]
    plt.plot(x_values, y_values, marker='o')

for another_one in all_left_hand_sides[1:]:
    all_left_hand_sides[0] = all_left_hand_sides[0].union(another_one)
combined = all_left_hand_sides[0]

plt.show()

coordinates=shapely.geometry.mapping(combined)['coordinates'][0]

import matplotlib.pyplot as plt
x_values = [coord[0] for coord in line_coords]
y_values = [coord[1] for coord in line_coords]
plt.plot(x_values, y_values, marker='o')
plt.grid(True)
plt.show()