From d71da6bb6bc6fcdd23cc13dbc64440f9df77d3cd Mon Sep 17 00:00:00 2001 From: Jonathan Date: Wed, 20 Oct 2021 01:09:51 +0100 Subject: [PATCH 1/7] adding tally accepting functions --- examples/plot_mesh_with_slice.py | 4 +- regular_mesh_plotter/__init__.py | 5 +- regular_mesh_plotter/core.py | 108 +++++++++++++++++++++++++++---- tests/test_plot_mesh.py | 8 +-- tests/test_slice_stl.py | 10 +-- 5 files changed, 108 insertions(+), 27 deletions(-) diff --git a/examples/plot_mesh_with_slice.py b/examples/plot_mesh_with_slice.py index e40d246..3d302bb 100644 --- a/examples/plot_mesh_with_slice.py +++ b/examples/plot_mesh_with_slice.py @@ -1,7 +1,7 @@ -from reg +from regular_mesh_plotter import plot_geometry_slice -plot_stl_slice("example.stl", plane_normal=[0, 0, 1]) +plot_geometry_slice("example.stl", plane_normal=[0, 0, 1]) diff --git a/regular_mesh_plotter/__init__.py b/regular_mesh_plotter/__init__.py index fe09240..12e3d54 100644 --- a/regular_mesh_plotter/__init__.py +++ b/regular_mesh_plotter/__init__.py @@ -1,4 +1,5 @@ -from .core import plot_mesh +from .core import plot_geometry_mesh +from .core import plot_regular_mesh_values from .core import get_tally_extent -from .core import plot_stl_slice +from .core import plot_regular_mesh_values_with_geometry diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index c233532..2742d0a 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -5,8 +5,8 @@ import numpy as np from pathlib import Path -def plot_stl_slice( - stl_or_mesh, +def plot_geometry_mesh( + mesh_file_or_trimesh_object, plane_origin: List[float] = None, plane_normal: List[float] = [0, 0, 1], rotate_plot: float = 0, @@ -28,12 +28,12 @@ def plot_stl_slice( A matplotlib.pyplot object """ - if isinstance(stl_or_mesh, str): - if not Path(stl_or_mesh).is_file(): - raise FileNotFoundError(f'file {stl_or_mesh} not found.') - mesh = trimesh.load_mesh(stl_or_mesh, process=False) + if isinstance(mesh_file_or_trimesh_object, str): + if not Path(mesh_file_or_trimesh_object).is_file(): + raise FileNotFoundError(f'file {mesh_file_or_trimesh_object} not found.') + mesh = trimesh.load_mesh(mesh_file_or_trimesh_object, process=False) else: - mesh = stl_or_mesh + mesh = mesh_file_or_trimesh_object if plane_origin is None: plane_origin = mesh.centroid @@ -65,7 +65,7 @@ def plot_stl_slice( return plt -def plot_mesh( +def plot_regular_mesh_values( values: np.ndarray, filename: Optional[str] = None, scale=None, # LogNorm(), @@ -100,6 +100,83 @@ def plot_mesh( # plt.close() return plt +def plot_regular_mesh_values_with_geometry( + values: np.ndarray, + mesh_file_or_trimesh_object, + filename: Optional[str] = None, + scale=None, # LogNorm(), + vmin=None, + label="", + extent=None, + x_label='X [cm]', + y_label='Y [cm]', + plane_origin: List[float] = None, + plane_normal: List[float] = [0, 0, 1], + rotate_plot: float = 0, +): + slice = plot_geometry_mesh( + mesh_file_or_trimesh_object=mesh_file_or_trimesh_object, + plane_origin = plane_origin, + plane_normal = plane_normal, + rotate_plot = rotate_plot + ) + + both = plot_regular_mesh_values( + values = values, + filename = filename, + scale=scale, # LogNorm(), + vmin=vmin, + label=label, + base_plt=slice, + extent=extent, + x_label=x_label, + y_label=y_label, + ) + + return both + +def plot_regular_mesh_tally( + tally: np.ndarray, + filename: Optional[str] = None, + scale=None, # LogNorm(), + vmin=None, + label="", + base_plt=None, + x_label='X [cm]', + y_label='Y [cm]', +): + + if base_plt: + plt = base_plt + else: + import matplotlib.pyplot as plt + plt.plot() + + import openmc_post_processor as opp + values = statepoint.process_tally( + tally=my_tally_1, + ) + + extent = get_tally_extent(tally) + + image_map = plt.imshow( + values, + norm=scale, + vmin=vmin, + extent=extent + ) + + plt.xlabel(x_label) + plt.ylabel(y_label) + + # image_map = fig.imshow(values, norm=scale, vmin=vmin) + plt.colorbar(image_map, label=label) + if filename: + plt.savefig(filename, dpi=300) + # fig.clear() + # plt.close() + return plt + def get_tally_extent( tally, ): @@ -107,17 +184,20 @@ def get_tally_extent( for filter in tally.filters: if isinstance(filter, openmc.MeshFilter): mesh_filter = filter - print(mesh_filter) - print(mesh_filter.mesh.lower_left) - print(mesh_filter.mesh.upper_right) - + # print(mesh_filter) + # print(mesh_filter.mesh.lower_left) + # print(mesh_filter.mesh.upper_right) + # print(mesh_filter.mesh.width) + # print(mesh_filter.mesh.__dict__) + # print(mesh_filter.mesh.dimension) + extent_x = (min(mesh_filter.mesh.lower_left[0],mesh_filter.mesh.upper_right[0]), max(mesh_filter.mesh.lower_left[0],mesh_filter.mesh.upper_right[0])) extent_y = (min(mesh_filter.mesh.lower_left[1],mesh_filter.mesh.upper_right[1]), max(mesh_filter.mesh.lower_left[1],mesh_filter.mesh.upper_right[1])) extent_z = (min(mesh_filter.mesh.lower_left[2],mesh_filter.mesh.upper_right[2]), max(mesh_filter.mesh.lower_left[2],mesh_filter.mesh.upper_right[2])) - if 1 in mesh_filter.mesh.width: + if 1 in mesh_filter.mesh.dimension.tolist(): print('2d mesh tally') - index_of_1d = mesh_filter.mesh.width.tolist().index(1) + index_of_1d = mesh_filter.mesh.dimension.tolist().index(1) print('index', index_of_1d) if index_of_1d == 0: return extent_y + extent_z diff --git a/tests/test_plot_mesh.py b/tests/test_plot_mesh.py index 6584e6d..d3ac81a 100644 --- a/tests/test_plot_mesh.py +++ b/tests/test_plot_mesh.py @@ -1,10 +1,10 @@ import unittest -from regular_mesh_plotter import plot_mesh +from regular_mesh_plotter import plot_regular_mesh_tally import numpy as np import matplotlib -class TestPlotSpectra(unittest.TestCase): +class TestPlotRegularMeshTally(unittest.TestCase): def setUp(self): self.values = np.array( @@ -72,8 +72,8 @@ def setUp(self): ] ) - def test_plot_mesh(self): + def test_plot_regular_mesh_tally(self): - test_plot = plot_mesh(values=self.values, filename='test.png') + test_plot = plot_regular_mesh_tally(values=self.values, filename='test.png') assert isinstance(test_plot, type(matplotlib.pyplot)) diff --git a/tests/test_slice_stl.py b/tests/test_slice_stl.py index e6480d9..43ad022 100644 --- a/tests/test_slice_stl.py +++ b/tests/test_slice_stl.py @@ -1,7 +1,7 @@ import unittest import trimesh -from regular_mesh_plotter import plot_stl_slice +from regular_mesh_plotter import plot_geometry_mesh class TestShape(unittest.TestCase): @@ -11,16 +11,16 @@ def setUp(self): def test_z_axis_slice(self): - plot_stl_slice(self.mesh, plane_normal=[0, 0, 1]) + plot_geometry_mesh(self.mesh, plane_normal=[0, 0, 1]) def test_offset_z_axis_slice(self): - plot_stl_slice(self.mesh, plane_origin=[0, 0, 10], plane_normal=[0, 0, 1]) + plot_geometry_mesh(self.mesh, plane_origin=[0, 0, 10], plane_normal=[0, 0, 1]) def test_x_axis_slice(self): - plot_stl_slice(self.mesh, plane_normal=[1, 0, 0]) + plot_geometry_mesh(self.mesh, plane_normal=[1, 0, 0]) def test_offset_x_axis_slice(self): - plot_stl_slice(self.mesh, plane_origin=[10, 0, 0], plane_normal=[1, 0, 0]) + plot_geometry_mesh(self.mesh, plane_origin=[10, 0, 0], plane_normal=[1, 0, 0]) From 31055b8775137f9287daff8c9a241afda80cc46c Mon Sep 17 00:00:00 2001 From: Jonathan Date: Wed, 20 Oct 2021 23:48:25 +0100 Subject: [PATCH 2/7] started new test strategy --- examples/plot_mesh_with_geometry_slice.py | 18 +++++++ examples/plot_mesh_with_slice.py | 7 --- regular_mesh_plotter/__init__.py | 5 +- regular_mesh_plotter/core.py | 49 ++++++++++++++++--- ...lice_stl.py => test_plot_geometry_mesh.py} | 0 tests/test_plot_regular_mesh_tally.py | 2 + ...t_plot_regular_mesh_tally_with_geometry.py | 2 + ...sh.py => test_plot_regular_mesh_values.py} | 4 +- ...plot_regular_mesh_values_with_geometry.py} | 2 +- 9 files changed, 72 insertions(+), 17 deletions(-) create mode 100644 examples/plot_mesh_with_geometry_slice.py delete mode 100644 examples/plot_mesh_with_slice.py rename tests/{test_slice_stl.py => test_plot_geometry_mesh.py} (100%) create mode 100644 tests/test_plot_regular_mesh_tally.py create mode 100644 tests/test_plot_regular_mesh_tally_with_geometry.py rename tests/{test_plot_mesh.py => test_plot_regular_mesh_values.py} (93%) rename tests/{test_plot_mesh_with_slice.py => test_plot_regular_mesh_values_with_geometry.py} (92%) diff --git a/examples/plot_mesh_with_geometry_slice.py b/examples/plot_mesh_with_geometry_slice.py new file mode 100644 index 0000000..eca6fc9 --- /dev/null +++ b/examples/plot_mesh_with_geometry_slice.py @@ -0,0 +1,18 @@ + + +from regular_mesh_plotter import plot_regular_mesh_tally_with_geometry +import openmc + +# loads in the statepoint file containing tallies +statepoint = openmc.StatePoint(filepath="statepoint.3.h5") + +# gets one tally from the available tallies +my_tally = statepoint.get_tally(name="neutron_effective_dose_on_2D_mesh_xy") + +# creates the matplotlib mesh plot with geometry +plot = plot_regular_mesh_tally_with_geometry( + tally=my_tally, + mesh_file_or_trimesh_object="dagmc.h5m", +) + +plot.show() \ No newline at end of file diff --git a/examples/plot_mesh_with_slice.py b/examples/plot_mesh_with_slice.py deleted file mode 100644 index 3d302bb..0000000 --- a/examples/plot_mesh_with_slice.py +++ /dev/null @@ -1,7 +0,0 @@ - - -from regular_mesh_plotter import plot_geometry_slice - - -plot_geometry_slice("example.stl", plane_normal=[0, 0, 1]) - diff --git a/regular_mesh_plotter/__init__.py b/regular_mesh_plotter/__init__.py index 12e3d54..fbb3b8d 100644 --- a/regular_mesh_plotter/__init__.py +++ b/regular_mesh_plotter/__init__.py @@ -1,5 +1,8 @@ from .core import plot_geometry_mesh from .core import plot_regular_mesh_values -from .core import get_tally_extent from .core import plot_regular_mesh_values_with_geometry +from .core import plot_regular_mesh_tally +from .core import plot_regular_mesh_tally_with_geometry +from .core import get_tally_extent + diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index 2742d0a..ad7470e 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -5,6 +5,9 @@ import numpy as np from pathlib import Path +import openmc_post_processor as opp +import openmc + def plot_geometry_mesh( mesh_file_or_trimesh_object, plane_origin: List[float] = None, @@ -64,7 +67,6 @@ def plot_geometry_mesh( plt.savefig(filename, dpi=300) return plt - def plot_regular_mesh_values( values: np.ndarray, filename: Optional[str] = None, @@ -135,8 +137,42 @@ def plot_regular_mesh_values_with_geometry( return both +def plot_regular_mesh_tally_with_geometry( + tally, + mesh_file_or_trimesh_object, + filename: Optional[str] = None, + scale=None, # LogNorm(), + vmin=None, + label="", + x_label='X [cm]', + y_label='Y [cm]', + plane_origin: List[float] = None, + plane_normal: List[float] = [0, 0, 1], + rotate_plot: float = 0, +): + + slice = plot_geometry_mesh( + mesh_file_or_trimesh_object=mesh_file_or_trimesh_object, + plane_origin = plane_origin, + plane_normal = plane_normal, + rotate_plot = rotate_plot + ) + + both = plot_regular_mesh_tally( + tally=tally, + filename=filename, + scale=scale, # LogNorm(), + vmin=vmin, + label=label, + base_plt=slice, + x_label=x_label, + y_label=y_label, + ) + + return both + def plot_regular_mesh_tally( - tally: np.ndarray, + tally, filename: Optional[str] = None, scale=None, # LogNorm(), vmin=None, @@ -152,13 +188,14 @@ def plot_regular_mesh_tally( import matplotlib.pyplot as plt plt.plot() - import openmc_post_processor as opp - values = statepoint.process_tally( - tally=my_tally_1, + values = opp.process_dose_tally( + tally=tally, ) extent = get_tally_extent(tally) + print('extent', extent) + image_map = plt.imshow( values, norm=scale, @@ -180,7 +217,7 @@ def plot_regular_mesh_tally( def get_tally_extent( tally, ): - import openmc + for filter in tally.filters: if isinstance(filter, openmc.MeshFilter): mesh_filter = filter diff --git a/tests/test_slice_stl.py b/tests/test_plot_geometry_mesh.py similarity index 100% rename from tests/test_slice_stl.py rename to tests/test_plot_geometry_mesh.py diff --git a/tests/test_plot_regular_mesh_tally.py b/tests/test_plot_regular_mesh_tally.py new file mode 100644 index 0000000..23d44aa --- /dev/null +++ b/tests/test_plot_regular_mesh_tally.py @@ -0,0 +1,2 @@ + +g \ No newline at end of file diff --git a/tests/test_plot_regular_mesh_tally_with_geometry.py b/tests/test_plot_regular_mesh_tally_with_geometry.py new file mode 100644 index 0000000..23d44aa --- /dev/null +++ b/tests/test_plot_regular_mesh_tally_with_geometry.py @@ -0,0 +1,2 @@ + +g \ No newline at end of file diff --git a/tests/test_plot_mesh.py b/tests/test_plot_regular_mesh_values.py similarity index 93% rename from tests/test_plot_mesh.py rename to tests/test_plot_regular_mesh_values.py index d3ac81a..e02bb29 100644 --- a/tests/test_plot_mesh.py +++ b/tests/test_plot_regular_mesh_values.py @@ -1,5 +1,5 @@ import unittest -from regular_mesh_plotter import plot_regular_mesh_tally +from regular_mesh_plotter import plot_regular_mesh_values import numpy as np import matplotlib @@ -74,6 +74,6 @@ def setUp(self): def test_plot_regular_mesh_tally(self): - test_plot = plot_regular_mesh_tally(values=self.values, filename='test.png') + test_plot = plot_regular_mesh_values(values=self.values, filename='test.png') assert isinstance(test_plot, type(matplotlib.pyplot)) diff --git a/tests/test_plot_mesh_with_slice.py b/tests/test_plot_regular_mesh_values_with_geometry.py similarity index 92% rename from tests/test_plot_mesh_with_slice.py rename to tests/test_plot_regular_mesh_values_with_geometry.py index 6d310db..8108b4b 100644 --- a/tests/test_plot_mesh_with_slice.py +++ b/tests/test_plot_regular_mesh_values_with_geometry.py @@ -1,4 +1,4 @@ - +# plot_regular_mesh_values_with_geometry # todo test the combination of mesh and stl slice # stl_slice = plot_stl_slice( # stl_or_mesh='steel.stl', From 89f3b60a62b92b3d2de886bf033ed1d627de8b1e Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 21 Oct 2021 17:52:35 +0100 Subject: [PATCH 3/7] making use of dagmc_geometry_slice_plotter --- examples/examples_from_readme.py | 77 +++++++++ regular_mesh_plotter/__init__.py | 1 - regular_mesh_plotter/core.py | 147 ++++++------------ setup.py | 3 +- tests/test_plot_geometry_mesh.py | 26 ---- tests/test_plot_regular_mesh_tally.py | 2 - ...t_plot_regular_mesh_tally_with_geometry.py | 2 - tests/test_plot_regular_mesh_values.py | 21 ++- ..._plot_regular_mesh_values_with_geometry.py | 118 +++++++++++--- 9 files changed, 239 insertions(+), 158 deletions(-) create mode 100644 examples/examples_from_readme.py delete mode 100644 tests/test_plot_geometry_mesh.py diff --git a/examples/examples_from_readme.py b/examples/examples_from_readme.py new file mode 100644 index 0000000..3c23db6 --- /dev/null +++ b/examples/examples_from_readme.py @@ -0,0 +1,77 @@ +# this examples requires additional python packages + +from dagmc_geometry_slice_plotter import plot_slice_of_dagmc_geometry +from stl_to_h5m import stl_to_h5m + +import paramak + +my_reactor = paramak.SubmersionTokamak( + inner_bore_radial_thickness=30, + inboard_tf_leg_radial_thickness=30, + center_column_shield_radial_thickness=30, + divertor_radial_thickness=80, + inner_plasma_gap_radial_thickness=50, + plasma_radial_thickness=200, + outer_plasma_gap_radial_thickness=50, + firstwall_radial_thickness=30, + blanket_rear_wall_radial_thickness=30, + number_of_tf_coils=16, + rotation_angle=360, + support_radial_thickness=90, + inboard_blanket_radial_thickness=30, + outboard_blanket_radial_thickness=30, + elongation=2.00, + triangularity=0.50, + pf_coil_case_thicknesses=[10, 10, 10, 10], + pf_coil_radial_thicknesses=[20, 50, 50, 20], + pf_coil_vertical_thicknesses=[20, 50, 50, 20], + pf_coil_radial_position=[500, 550, 550, 500], + pf_coil_vertical_position=[270, 100, -100, -270], + rear_blanket_to_tf_gap=50, + outboard_tf_coil_radial_thickness=30, + outboard_tf_coil_poloidal_thickness=30, +) + + +stl_filenames = my_reactor.export_stl() + + +stl_to_h5m( + files_with_tags=[ + (stl_filename, name) + for name, stl_filename in zip(my_reactor.name, stl_filenames) + ], + h5m_filename="dagmc.h5m", +) + + +plot = plot_slice_of_dagmc_geometry( + dagmc_file_or_trimesh_object="dagmc.h5m", + plane_normal=[0, 0, 1], + output_filename="my_plot1.png", +) + + +plot = plot_slice_of_dagmc_geometry( + dagmc_file_or_trimesh_object="dagmc.h5m", + plane_origin=[0, 0, 300], + plane_normal=[0, 0, 1], + output_filename="my_plot2.png", +) + + +plot = plot_slice_of_dagmc_geometry( + dagmc_file_or_trimesh_object="dagmc.h5m", + plane_normal=[0, 1, 0], + rotate_plot=45, + output_filename="my_plot3.png", +) + +plot.savefig("big.png", dpi=600) + +plot_slice_of_dagmc_geometry( + dagmc_file_or_trimesh_object="dagmc.h5m", + plane_normal=[1, 0, 0], + rotate_plot=270, + output_filename="my_plot4.png", +) diff --git a/regular_mesh_plotter/__init__.py b/regular_mesh_plotter/__init__.py index fbb3b8d..49becfc 100644 --- a/regular_mesh_plotter/__init__.py +++ b/regular_mesh_plotter/__init__.py @@ -1,5 +1,4 @@ -from .core import plot_geometry_mesh from .core import plot_regular_mesh_values from .core import plot_regular_mesh_values_with_geometry from .core import plot_regular_mesh_tally diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index ad7470e..2e510b2 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -4,68 +4,11 @@ from typing import List, Optional import numpy as np from pathlib import Path +import dagmc_geometry_slice_plotter as dgsp import openmc_post_processor as opp import openmc -def plot_geometry_mesh( - mesh_file_or_trimesh_object, - plane_origin: List[float] = None, - plane_normal: List[float] = [0, 0, 1], - rotate_plot: float = 0, - filename: Optional[str] = None, -) -> plt: - """Slices through a 3D geometry in STL file format and extracts a slice of - the geometry at the provided plane and origin - - Args: - plane_origin: the origin of the plain, if None then the mesh.centroid - will be used. - plane_normal: the plane to slice the geometry on. Defaults to slicing - along the Z plane which is input as [0, 0, 1]. - rotate_plot: the angle in degrees to rotate the plot by. Useful when - used in conjunction with changing plane_normal to orientate the - plot correctly. - - Return: - A matplotlib.pyplot object - """ - - if isinstance(mesh_file_or_trimesh_object, str): - if not Path(mesh_file_or_trimesh_object).is_file(): - raise FileNotFoundError(f'file {mesh_file_or_trimesh_object} not found.') - mesh = trimesh.load_mesh(mesh_file_or_trimesh_object, process=False) - else: - mesh = mesh_file_or_trimesh_object - - if plane_origin is None: - plane_origin = mesh.centroid - slice = mesh.section( - plane_origin=plane_origin, - plane_normal=plane_normal, - ) - - slice_2D, to_3D = slice.to_planar() - - # keep plot axis scaled the same - plt.axes().set_aspect("equal", "datalim") - - if rotate_plot != 0: - base = plt.gca().transData - rot = transforms.Affine2D().rotate_deg(rotate_plot) - - for entity in slice_2D.entities: - - discrete = entity.discrete(slice_2D.vertices) - - if rotate_plot != 0: - plt.plot(*discrete.T, color="black", linewidth=1, transform=rot + base) - else: - plt.plot(*discrete.T, color="black", linewidth=1) - - if filename: - plt.savefig(filename, dpi=300) - return plt def plot_regular_mesh_values( values: np.ndarray, @@ -75,21 +18,17 @@ def plot_regular_mesh_values( label="", base_plt=None, extent=None, - x_label='X [cm]', - y_label='Y [cm]', + x_label="X [cm]", + y_label="Y [cm]", ): if base_plt: plt = base_plt else: import matplotlib.pyplot as plt + plt.plot() - image_map = plt.imshow( - values, - norm=scale, - vmin=vmin, - extent=extent - ) + image_map = plt.imshow(values, norm=scale, vmin=vmin, extent=extent) plt.xlabel(x_label) plt.ylabel(y_label) @@ -102,30 +41,32 @@ def plot_regular_mesh_values( # plt.close() return plt + def plot_regular_mesh_values_with_geometry( values: np.ndarray, - mesh_file_or_trimesh_object, + dagmc_file_or_trimesh_object, filename: Optional[str] = None, scale=None, # LogNorm(), vmin=None, label="", extent=None, - x_label='X [cm]', - y_label='Y [cm]', + x_label="X [cm]", + y_label="Y [cm]", plane_origin: List[float] = None, plane_normal: List[float] = [0, 0, 1], rotate_plot: float = 0, ): - slice = plot_geometry_mesh( - mesh_file_or_trimesh_object=mesh_file_or_trimesh_object, - plane_origin = plane_origin, - plane_normal = plane_normal, - rotate_plot = rotate_plot + + slice = dgsp.plot_slice_of_dagmc_geometry( + dagmc_file_or_trimesh_object=dagmc_file_or_trimesh_object, + plane_origin=plane_origin, + plane_normal=plane_normal, + rotate_plot=rotate_plot, ) both = plot_regular_mesh_values( - values = values, - filename = filename, + values=values, + filename=filename, scale=scale, # LogNorm(), vmin=vmin, label=label, @@ -137,25 +78,26 @@ def plot_regular_mesh_values_with_geometry( return both + def plot_regular_mesh_tally_with_geometry( tally, - mesh_file_or_trimesh_object, + dagmc_file_or_trimesh_object, filename: Optional[str] = None, scale=None, # LogNorm(), vmin=None, label="", - x_label='X [cm]', - y_label='Y [cm]', + x_label="X [cm]", + y_label="Y [cm]", plane_origin: List[float] = None, plane_normal: List[float] = [0, 0, 1], rotate_plot: float = 0, ): - slice = plot_geometry_mesh( - mesh_file_or_trimesh_object=mesh_file_or_trimesh_object, - plane_origin = plane_origin, - plane_normal = plane_normal, - rotate_plot = rotate_plot + slice = dgsp.plot_slice_of_dagmc_geometry( + dagmc_file_or_trimesh_object=dagmc_file_or_trimesh_object, + plane_origin=plane_origin, + plane_normal=plane_normal, + rotate_plot=rotate_plot, ) both = plot_regular_mesh_tally( @@ -171,6 +113,7 @@ def plot_regular_mesh_tally_with_geometry( return both + def plot_regular_mesh_tally( tally, filename: Optional[str] = None, @@ -178,14 +121,15 @@ def plot_regular_mesh_tally( vmin=None, label="", base_plt=None, - x_label='X [cm]', - y_label='Y [cm]', + x_label="X [cm]", + y_label="Y [cm]", ): if base_plt: plt = base_plt else: import matplotlib.pyplot as plt + plt.plot() values = opp.process_dose_tally( @@ -194,14 +138,9 @@ def plot_regular_mesh_tally( extent = get_tally_extent(tally) - print('extent', extent) + print("extent", extent) - image_map = plt.imshow( - values, - norm=scale, - vmin=vmin, - extent=extent - ) + image_map = plt.imshow(values, norm=scale, vmin=vmin, extent=extent) plt.xlabel(x_label) plt.ylabel(y_label) @@ -214,6 +153,7 @@ def plot_regular_mesh_tally( # plt.close() return plt + def get_tally_extent( tally, ): @@ -228,14 +168,23 @@ def get_tally_extent( # print(mesh_filter.mesh.__dict__) # print(mesh_filter.mesh.dimension) - extent_x = (min(mesh_filter.mesh.lower_left[0],mesh_filter.mesh.upper_right[0]), max(mesh_filter.mesh.lower_left[0],mesh_filter.mesh.upper_right[0])) - extent_y = (min(mesh_filter.mesh.lower_left[1],mesh_filter.mesh.upper_right[1]), max(mesh_filter.mesh.lower_left[1],mesh_filter.mesh.upper_right[1])) - extent_z = (min(mesh_filter.mesh.lower_left[2],mesh_filter.mesh.upper_right[2]), max(mesh_filter.mesh.lower_left[2],mesh_filter.mesh.upper_right[2])) - + extent_x = ( + min(mesh_filter.mesh.lower_left[0], mesh_filter.mesh.upper_right[0]), + max(mesh_filter.mesh.lower_left[0], mesh_filter.mesh.upper_right[0]), + ) + extent_y = ( + min(mesh_filter.mesh.lower_left[1], mesh_filter.mesh.upper_right[1]), + max(mesh_filter.mesh.lower_left[1], mesh_filter.mesh.upper_right[1]), + ) + extent_z = ( + min(mesh_filter.mesh.lower_left[2], mesh_filter.mesh.upper_right[2]), + max(mesh_filter.mesh.lower_left[2], mesh_filter.mesh.upper_right[2]), + ) + if 1 in mesh_filter.mesh.dimension.tolist(): - print('2d mesh tally') + print("2d mesh tally") index_of_1d = mesh_filter.mesh.dimension.tolist().index(1) - print('index', index_of_1d) + print("index", index_of_1d) if index_of_1d == 0: return extent_y + extent_z if index_of_1d == 1: diff --git a/setup.py b/setup.py index bcd4020..42cd22d 100644 --- a/setup.py +++ b/setup.py @@ -37,6 +37,7 @@ "matplotlib>=3.4.2", "trimesh", "shapely", - "scipy" + "scipy", + "dagmc_geometry_slice_plotter" ], ) diff --git a/tests/test_plot_geometry_mesh.py b/tests/test_plot_geometry_mesh.py deleted file mode 100644 index 43ad022..0000000 --- a/tests/test_plot_geometry_mesh.py +++ /dev/null @@ -1,26 +0,0 @@ -import unittest - -import trimesh -from regular_mesh_plotter import plot_geometry_mesh - - -class TestShape(unittest.TestCase): - def setUp(self): - - self.mesh = trimesh.load_mesh("tests/example.stl", process=False) - - def test_z_axis_slice(self): - - plot_geometry_mesh(self.mesh, plane_normal=[0, 0, 1]) - - def test_offset_z_axis_slice(self): - - plot_geometry_mesh(self.mesh, plane_origin=[0, 0, 10], plane_normal=[0, 0, 1]) - - def test_x_axis_slice(self): - - plot_geometry_mesh(self.mesh, plane_normal=[1, 0, 0]) - - def test_offset_x_axis_slice(self): - - plot_geometry_mesh(self.mesh, plane_origin=[10, 0, 0], plane_normal=[1, 0, 0]) diff --git a/tests/test_plot_regular_mesh_tally.py b/tests/test_plot_regular_mesh_tally.py index 23d44aa..e69de29 100644 --- a/tests/test_plot_regular_mesh_tally.py +++ b/tests/test_plot_regular_mesh_tally.py @@ -1,2 +0,0 @@ - -g \ No newline at end of file diff --git a/tests/test_plot_regular_mesh_tally_with_geometry.py b/tests/test_plot_regular_mesh_tally_with_geometry.py index 23d44aa..e69de29 100644 --- a/tests/test_plot_regular_mesh_tally_with_geometry.py +++ b/tests/test_plot_regular_mesh_tally_with_geometry.py @@ -1,2 +0,0 @@ - -g \ No newline at end of file diff --git a/tests/test_plot_regular_mesh_values.py b/tests/test_plot_regular_mesh_values.py index e02bb29..834d98b 100644 --- a/tests/test_plot_regular_mesh_values.py +++ b/tests/test_plot_regular_mesh_values.py @@ -1,10 +1,13 @@ +import os import unittest -from regular_mesh_plotter import plot_regular_mesh_values -import numpy as np +from pathlib import Path + import matplotlib +import numpy as np +from regular_mesh_plotter import plot_regular_mesh_values -class TestPlotRegularMeshTally(unittest.TestCase): +class TestPlotRegularMeshValues(unittest.TestCase): def setUp(self): self.values = np.array( @@ -72,8 +75,16 @@ def setUp(self): ] ) - def test_plot_regular_mesh_tally(self): + def test_plot_regular_mesh_values(self): - test_plot = plot_regular_mesh_values(values=self.values, filename='test.png') + test_plot = plot_regular_mesh_values(values=self.values) assert isinstance(test_plot, type(matplotlib.pyplot)) + + def test_plot_regular_mesh_values_with_output(self): + + os.system('rm test.png') + + plot_regular_mesh_values(values=self.values, filename='test.png') + + assert Path('test.png').is_file() diff --git a/tests/test_plot_regular_mesh_values_with_geometry.py b/tests/test_plot_regular_mesh_values_with_geometry.py index 8108b4b..a8c934b 100644 --- a/tests/test_plot_regular_mesh_values_with_geometry.py +++ b/tests/test_plot_regular_mesh_values_with_geometry.py @@ -1,22 +1,96 @@ -# plot_regular_mesh_values_with_geometry -# todo test the combination of mesh and stl slice -# stl_slice = plot_stl_slice( -# stl_or_mesh='steel.stl', -# plane_origin = None, -# plane_normal = [0, 0, 1], -# rotate_plot = 0, -# filename='slice.png' -# ) - -# extent = get_tally_extent(my_tally) -# print(extent) -# plot_mesh( -# extent=extent, -# values= result, -# scale=None, # LogNorm(), -# vmin=None, -# label="picosievert / cm / pulse", -# base_plt=stl_slice, -# filename= 'test.png', -# # vmin=1e6, -# ) \ No newline at end of file +import os +import unittest +from pathlib import Path + +import matplotlib +import numpy as np +from regular_mesh_plotter import plot_regular_mesh_values_with_geometry + + +class TestPlotRegularMeshValuesWithGeometry(unittest.TestCase): + def setUp(self): + + self.values = np.array( + [ + [ + 1.64565869e09, + 2.59505642e09, + 3.06732422e09, + 2.64141597e09, + 1.76520128e09, + 2.14909584e09, + 3.70513916e09, + 4.99737560e09, + 3.86795536e09, + 2.21272467e09, + ], + [ + 2.92406594e09, + 5.91396360e09, + 9.39595883e09, + 5.94102629e09, + 2.78174231e09, + 3.47563407e09, + 1.02570496e10, + 2.50416310e10, + 1.01248003e10, + 3.34937674e09, + ], + [ + 3.99684099e09, + 1.72147289e10, + 3.68431465e11, + 1.64968908e10, + 3.74168705e09, + 3.90640820e09, + 1.65505774e10, + 3.66837062e11, + 1.65568070e10, + 3.65646431e09, + ], + [ + 3.27796129e09, + 1.01456714e10, + 2.45757058e10, + 1.00180127e10, + 3.45105436e09, + 2.63214911e09, + 5.95924816e09, + 9.43589769e09, + 5.71740137e09, + 2.75163850e09, + ], + [ + 2.37324680e09, + 3.77339226e09, + 5.01889988e09, + 3.58500172e09, + 2.16754228e09, + 1.81599509e09, + 2.57229036e09, + 3.09622197e09, + 2.50136006e09, + 1.72280196e09, + ], + ] + ) + + def test_plot_regular_mesh_values_with_geometry(self): + + test_plot = plot_regular_mesh_values_with_geometry( + values=self.values, + dagmc_file_or_trimesh_object='example.stl' # this could be a h5m file + ) + + assert isinstance(test_plot, type(matplotlib.pyplot)) + + def test_plot_regular_mesh_values_with_geometry(self): + + test_plot = plot_regular_mesh_values_with_geometry( + values=self.values, + dagmc_file_or_trimesh_object='example.stl', # this could be a h5m file + filename='test.png' + ) + + assert isinstance(test_plot, type(matplotlib.pyplot)) + assert Path('test.png').is_file() From aa4370684a4826d4c150bb37174676d20c423db8 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 21 Oct 2021 18:14:04 +0100 Subject: [PATCH 4/7] added openmc_post_processor --- setup.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 42cd22d..cb86a3f 100644 --- a/setup.py +++ b/setup.py @@ -38,6 +38,7 @@ "trimesh", "shapely", "scipy", - "dagmc_geometry_slice_plotter" + "dagmc_geometry_slice_plotter", + "openmc_post_processor" ], ) From fefd700cc11b2ded20e52f1a03fc4d68b1c0ca88 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 21 Oct 2021 18:18:40 +0100 Subject: [PATCH 5/7] added openmc to ci --- .github/workflows/ci_with_install.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/ci_with_install.yml b/.github/workflows/ci_with_install.yml index fb7ba6c..fbc3ce0 100644 --- a/.github/workflows/ci_with_install.yml +++ b/.github/workflows/ci_with_install.yml @@ -29,6 +29,10 @@ jobs: run: | python setup.py install + - name: install openmc + run: | + conda install -c conda-forge openmc + - name: install packages for tests run: | pip install -r requirements-test.txt From 003f6fed1607863ecb9a8b5f8b55152bc1e9a3ee Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 22 Oct 2021 17:38:42 +0100 Subject: [PATCH 6/7] added ability to rotate mesh plot --- regular_mesh_plotter/core.py | 63 ++++++++++++++++++++---------------- 1 file changed, 36 insertions(+), 27 deletions(-) diff --git a/regular_mesh_plotter/core.py b/regular_mesh_plotter/core.py index 2e510b2..2ff1f32 100644 --- a/regular_mesh_plotter/core.py +++ b/regular_mesh_plotter/core.py @@ -5,6 +5,7 @@ import numpy as np from pathlib import Path import dagmc_geometry_slice_plotter as dgsp +from matplotlib import transforms import openmc_post_processor as opp import openmc @@ -20,6 +21,7 @@ def plot_regular_mesh_values( extent=None, x_label="X [cm]", y_label="Y [cm]", + rotate_plot: float = 0, ): if base_plt: @@ -28,7 +30,18 @@ def plot_regular_mesh_values( import matplotlib.pyplot as plt plt.plot() - image_map = plt.imshow(values, norm=scale, vmin=vmin, extent=extent) + + if rotate_plot != 0: + x_center = sum(extent[:2]) / 2 + y_center = sum(extent[2:]) / 2 + base = plt.gca().transData + rot = transforms.Affine2D().rotate_deg_around(x_center, y_center, rotate_plot) + + image_map = plt.imshow( + values, norm=scale, vmin=vmin, extent=extent, transform=rot + base + ) + else: + image_map = plt.imshow(values, norm=scale, vmin=vmin, extent=extent) plt.xlabel(x_label) plt.ylabel(y_label) @@ -37,8 +50,6 @@ def plot_regular_mesh_values( plt.colorbar(image_map, label=label) if filename: plt.savefig(filename, dpi=300) - # fig.clear() - # plt.close() return plt @@ -54,14 +65,15 @@ def plot_regular_mesh_values_with_geometry( y_label="Y [cm]", plane_origin: List[float] = None, plane_normal: List[float] = [0, 0, 1], - rotate_plot: float = 0, + rotate_mesh: float = 0, + rotate_geometry: float = 0, ): slice = dgsp.plot_slice_of_dagmc_geometry( dagmc_file_or_trimesh_object=dagmc_file_or_trimesh_object, plane_origin=plane_origin, plane_normal=plane_normal, - rotate_plot=rotate_plot, + rotate_plot=rotate_geometry, ) both = plot_regular_mesh_values( @@ -74,6 +86,7 @@ def plot_regular_mesh_values_with_geometry( extent=extent, x_label=x_label, y_label=y_label, + rotate_plot=rotate_mesh, ) return both @@ -90,14 +103,15 @@ def plot_regular_mesh_tally_with_geometry( y_label="Y [cm]", plane_origin: List[float] = None, plane_normal: List[float] = [0, 0, 1], - rotate_plot: float = 0, + rotate_mesh: float = 0, + rotate_geometry: float = 0, ): slice = dgsp.plot_slice_of_dagmc_geometry( dagmc_file_or_trimesh_object=dagmc_file_or_trimesh_object, plane_origin=plane_origin, plane_normal=plane_normal, - rotate_plot=rotate_plot, + rotate_plot=rotate_geometry, ) both = plot_regular_mesh_tally( @@ -109,6 +123,7 @@ def plot_regular_mesh_tally_with_geometry( base_plt=slice, x_label=x_label, y_label=y_label, + rotate_plot=rotate_mesh, ) return both @@ -123,35 +138,29 @@ def plot_regular_mesh_tally( base_plt=None, x_label="X [cm]", y_label="Y [cm]", + rotate_plot: float = 0, ): - if base_plt: - plt = base_plt - else: - import matplotlib.pyplot as plt - - plt.plot() - values = opp.process_dose_tally( tally=tally, ) extent = get_tally_extent(tally) - print("extent", extent) - - image_map = plt.imshow(values, norm=scale, vmin=vmin, extent=extent) - - plt.xlabel(x_label) - plt.ylabel(y_label) + plot = plot_regular_mesh_values( + values=values, + filename=filename, + scale=scale, + vmin=vmin, + label=label, + base_plt=base_plt, + extent=extent, + x_label=x_label, + y_label=y_label, + rotate_plot=rotate_plot, + ) - # image_map = fig.imshow(values, norm=scale, vmin=vmin) - plt.colorbar(image_map, label=label) - if filename: - plt.savefig(filename, dpi=300) - # fig.clear() - # plt.close() - return plt + return plot def get_tally_extent( From 2228beaffeb71db47317ea2bcf6f2d6102d4978f Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 22 Oct 2021 18:02:22 +0100 Subject: [PATCH 7/7] added rotate mesh option and started readme --- README.md | 66 ++++++++++++++++++++++- examples/plot_mesh_with_geometry_slice.py | 3 +- 2 files changed, 65 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index ba746a1..59d7134 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,64 @@ -# dose_map_plotter -A Python package for plotting neutronics dose maps +[![N|Python](https://www.python.org/static/community_logos/python-powered-w-100x40.png)](https://www.python.org) + +[![CI with install](https://github.com/fusion-energy/regular_mesh_plotter/actions/workflows/ci_with_install.yml/badge.svg?branch=develop)](https://github.com/fusion-energy/regular_mesh_plotter/actions/workflows/ci_with_install.yml) + +[![PyPI](https://img.shields.io/pypi/v/regular-mesh-plotter?color=brightgreen&label=pypi&logo=grebrightgreenen&logoColor=green)](https://pypi.org/project/regular-mesh-plotter/) + +[![codecov](https://codecov.io/gh/fusion-energy/regular_mesh_plotter/branch/main/graph/badge.svg)](https://codecov.io/gh/fusion-energy/regular_mesh_plotter) + +## A minimal Python package that plots 2D mesh tally results with the underlying DAGMC geometry + +# Installation + +```bash +pip install regular_mesh_plotter +``` + +Mesh results in the form of Numpy arrays or OpenMC.tally objects can be plotted +with a single API call. + +A Matplotlib.pyplot object is returned by all functions so one can make changes +to the legend, axis, colour map etc. However some key options are accessable +in the function call directly. + +There are additional options that allow + +- rotation of the mesh tally results +- rotation of the DAGMC geometry slice +- saving the plot as an image file +- specifying contour lines TODO +- changing axis and colour bar labels +- changing colour scale applied +- truncation of values +- The plane_normal of the DAGMC geometry + +The resulting plots can be used to show dose maps, activation, reaction rate +and other mesh tally results. + + +Example 1 shows a Numpy array plotted +```python +TODO +``` + +Example 2 shows a Numpy array plotted with an underlying DAGMC geometry +```python +TODO +``` + +Example 3 shows a OpenMC tally plotted with an underlying DAGMC geometry +```python +TODO +``` + +Example 4 shows how to rotate the underlying DAGMC geometry and mesh tally. +This is sometimes necessary as the slice and mesh can get out of alignment +when changing the +```python +TODO +``` + +# Related packages + +If you want to plot the DAGMC geometry without a mesh tally then take a look at +the [dagmc_geometry_slice_plotter](https://github.com/fusion-energy/dagmc_geometry_slice_plotter) package diff --git a/examples/plot_mesh_with_geometry_slice.py b/examples/plot_mesh_with_geometry_slice.py index eca6fc9..fb2622c 100644 --- a/examples/plot_mesh_with_geometry_slice.py +++ b/examples/plot_mesh_with_geometry_slice.py @@ -1,5 +1,4 @@ - from regular_mesh_plotter import plot_regular_mesh_tally_with_geometry import openmc @@ -15,4 +14,4 @@ mesh_file_or_trimesh_object="dagmc.h5m", ) -plot.show() \ No newline at end of file +plot.show()