From fc615451d04c358d0f32bb72837fc52992fb39b0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 15:54:49 +0200 Subject: [PATCH 01/53] added method to Tally --- openmc/tallies.py | 86 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index d0355f14ee..f559c392cc 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -8,6 +8,8 @@ from pathlib import Path from xml.etree import ElementTree as ET import h5py +import vtk +import vtk.util.numpy_support as nps import numpy as np import pandas as pd import scipy.sparse as sps @@ -457,6 +459,54 @@ class Tally(IDManagerMixin): self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) self._sparse = False + def write_to_vtk(self, filename): + mesh = None + for f in self.filters: + if isinstance(f, openmc.MeshFilter): + if isinstance(f.mesh, (openmc.RegularMesh, openmc.CylindricalMesh)): + mesh = f.mesh + break + + if not mesh: + raise ValueError( + "write_to_vtk only works with openmc.RegularMesh, openmc.CylindricalMesh" + ) + + if isinstance(mesh, openmc.RegularMesh): + vtk_grid = voxels_to_vtk( + x_vals=np.linspace( + mesh.lower_left[0], + mesh.upper_right[0], + num=mesh.dimension[0] + 1, + ), + y_vals=np.linspace( + mesh.lower_left[1], + mesh.upper_right[1], + num=mesh.dimension[1] + 1, + ), + z_vals=np.linspace( + mesh.lower_left[2], + mesh.upper_right[2], + num=mesh.dimension[2] + 1, + ), + mean=self.mean, + std_dev=self.std_dev, + cylindrical=False, + ) + elif isinstance(mesh, openmc.CylindricalMesh): + vtk_grid = voxels_to_vtk( + x_vals=mesh.r_grid, + y_vals=mesh.phi_grid, + z_vals=mesh.z_grid, + mean=self.mean, + std_dev=self.std_dev, + cylindrical=True, + ) + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + def remove_score(self, score): """Remove a score from the tally @@ -3229,3 +3279,39 @@ class Tallies(cv.CheckedList): tallies.append(tally) return cls(tallies) + + +def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): + vtk_box = vtk.vtkStructuredGrid() + + vtk_box.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + + # points + points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + if cylindrical: + points_cartesian = np.copy(points) + r = points[:, 0] + phi = points[:, 1] + z = points[:, 2] + points_cartesian[:, 0] = r * np.cos(phi) + points_cartesian[:, 1] = r * np.sin(phi) + points_cartesian[:, 2] = z + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) + + vtk_box.SetPoints(vtkPts) + + # # data + mean_array = vtk.vtkDoubleArray() + mean_array.SetName("mean") + mean_array.SetArray(mean, mean.size, True) + + std_dev_array = vtk.vtkDoubleArray() + std_dev_array.SetName("std_dev") + std_dev_array.SetArray(std_dev, std_dev.size, True) + + vtk_box.GetCellData().AddArray(mean_array) + vtk_box.GetCellData().AddArray(std_dev_array) + + return vtk_box From e870ab2040e6d2a9894eb99e413a5174f83ae546 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:00:54 +0200 Subject: [PATCH 02/53] added RectilinearMesh --- openmc/tallies.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index f559c392cc..20e6e39964 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -463,13 +463,13 @@ class Tally(IDManagerMixin): mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): - if isinstance(f.mesh, (openmc.RegularMesh, openmc.CylindricalMesh)): + if isinstance(f.mesh, (openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh)): mesh = f.mesh break if not mesh: raise ValueError( - "write_to_vtk only works with openmc.RegularMesh, openmc.CylindricalMesh" + "write_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" ) if isinstance(mesh, openmc.RegularMesh): @@ -493,6 +493,15 @@ class Tally(IDManagerMixin): std_dev=self.std_dev, cylindrical=False, ) + if isinstance(mesh, openmc.RectilinearMesh): + vtk_grid = voxels_to_vtk( + x_vals=mesh.x_grid, + y_vals=mesh.y_grid, + z_vals=mesh.z_grid, + mean=self.mean, + std_dev=self.std_dev, + cylindrical=False, + ) elif isinstance(mesh, openmc.CylindricalMesh): vtk_grid = voxels_to_vtk( x_vals=mesh.r_grid, From 06c5d3ec584c9ef9a62b08c561945b34e2e6b4e9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:04:33 +0200 Subject: [PATCH 03/53] added comments --- openmc/tallies.py | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 20e6e39964..6797253e1c 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -460,6 +460,7 @@ class Tally(IDManagerMixin): self._sparse = False def write_to_vtk(self, filename): + # check that tally has a MeshFilter mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): @@ -511,6 +512,8 @@ class Tally(IDManagerMixin): std_dev=self.std_dev, cylindrical=True, ) + + # write the .vtk file writer = vtk.vtkStructuredGridWriter() writer.SetFileName(filename) writer.SetInputData(vtk_grid) @@ -3291,36 +3294,32 @@ class Tallies(cv.CheckedList): def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): - vtk_box = vtk.vtkStructuredGrid() + vtk_grid = vtk.vtkStructuredGrid() - vtk_box.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) - # points + # create points points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - if cylindrical: + if cylindrical: # transform points to cartesian coordinates points_cartesian = np.copy(points) - r = points[:, 0] - phi = points[:, 1] - z = points[:, 2] + r, phi, z = points[:, 0], points[:, 1], points[:, 2] points_cartesian[:, 0] = r * np.cos(phi) points_cartesian[:, 1] = r * np.sin(phi) points_cartesian[:, 2] = z vtkPts = vtk.vtkPoints() vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) + vtk_grid.SetPoints(vtkPts) - vtk_box.SetPoints(vtkPts) - - # # data + # add mean and std dev data mean_array = vtk.vtkDoubleArray() mean_array.SetName("mean") mean_array.SetArray(mean, mean.size, True) + vtk_grid.GetCellData().AddArray(mean_array) std_dev_array = vtk.vtkDoubleArray() std_dev_array.SetName("std_dev") std_dev_array.SetArray(std_dev, std_dev.size, True) + vtk_grid.GetCellData().AddArray(std_dev_array) - vtk_box.GetCellData().AddArray(mean_array) - vtk_box.GetCellData().AddArray(std_dev_array) - - return vtk_box + return vtk_grid From 30b06fa2910a51ac04a249170bff27d54cb3f26a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:08:36 +0200 Subject: [PATCH 04/53] docstrings --- openmc/tallies.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 6797253e1c..77f0b8cd39 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -460,6 +460,15 @@ class Tally(IDManagerMixin): self._sparse = False def write_to_vtk(self, filename): + """Writes the tally to a vtk file + + Args: + filename (str): the filename (must end with .vtk) + + Raises: + ValueError: if no MeshFilter with appropriate mesh was found + (SphericalMesh not supported) + """ # check that tally has a MeshFilter mesh = None for f in self.filters: @@ -512,7 +521,7 @@ class Tally(IDManagerMixin): std_dev=self.std_dev, cylindrical=True, ) - + # write the .vtk file writer = vtk.vtkStructuredGridWriter() writer.SetFileName(filename) @@ -3294,6 +3303,20 @@ class Tallies(cv.CheckedList): def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): + """Creates a vtk object from a list of X, Y, Z values and mean/std_dev data. + + Args: + x_vals (list): X values. + y_vals (list): Y values. + z_vals (list): Z values. + mean (np.array): the tally mean. + std_dev (np.array): the tally standard deviation. + cylindrical (bool, optional): If set to True, cylindrical coordinates + (r, phi, z) are assumed. Defaults to True. + + Returns: + vtkStructuredGrid: a vtk object containing tally data on the appropriate grid + """ vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) From 06c6bcd7d269aa3b3317c8d64f9cf4318170468c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:09:32 +0200 Subject: [PATCH 05/53] added TODO --- openmc/tallies.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index 77f0b8cd39..c0f410524f 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3317,6 +3317,9 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): Returns: vtkStructuredGrid: a vtk object containing tally data on the appropriate grid """ + + # TODO: should this be a method of Tally? + vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) From 877a62c0af072e14c3b62991f21a6648352dc8bb Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:43:42 +0200 Subject: [PATCH 06/53] fixed cylindrical --- openmc/tallies.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index c0f410524f..2e31e8ee26 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3332,6 +3332,7 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): points_cartesian[:, 0] = r * np.cos(phi) points_cartesian[:, 1] = r * np.sin(phi) points_cartesian[:, 2] = z + points = points_cartesian vtkPts = vtk.vtkPoints() vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) From 1eb8670f2b93cbbb42912dd7903c1d1ae18b1d33 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:46:10 +0200 Subject: [PATCH 07/53] mean and std_dev are np.zeros if None (for testing without running a sim) --- openmc/tallies.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index 2e31e8ee26..75c7ea8f3b 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3341,11 +3341,23 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): # add mean and std dev data mean_array = vtk.vtkDoubleArray() mean_array.SetName("mean") + if mean is None: + mean = np.zeros( + (len(x_vals) - 1) + * (len(y_vals) - 1) + * (len(z_vals) - 1) + ) mean_array.SetArray(mean, mean.size, True) vtk_grid.GetCellData().AddArray(mean_array) std_dev_array = vtk.vtkDoubleArray() std_dev_array.SetName("std_dev") + if std_dev is None: + std_dev = np.zeros( + (len(x_vals) - 1) + * (len(y_vals) - 1) + * (len(z_vals) - 1) + ) std_dev_array.SetArray(std_dev, std_dev.size, True) vtk_grid.GetCellData().AddArray(std_dev_array) From 782d987fc836494166ffb66ab9e9f81888666582 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Tue, 28 Jun 2022 20:18:52 +0200 Subject: [PATCH 08/53] Optional vtk import Co-authored-by: Jonathan Shimwell --- openmc/tallies.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 75c7ea8f3b..b795fc83a2 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -8,8 +8,6 @@ from pathlib import Path from xml.etree import ElementTree as ET import h5py -import vtk -import vtk.util.numpy_support as nps import numpy as np import pandas as pd import scipy.sparse as sps @@ -470,6 +468,12 @@ class Tally(IDManagerMixin): (SphericalMesh not supported) """ # check that tally has a MeshFilter + try: + import vtk + import vtk.util.numpy_support as nps + except ImportError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ImportError(msg) mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): From 6eaab920eeaddac452ad07e25015357c3658f139 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:18:10 +0200 Subject: [PATCH 09/53] better error catching --- openmc/tallies.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index b795fc83a2..faa7c3efe0 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -471,9 +471,12 @@ class Tally(IDManagerMixin): try: import vtk import vtk.util.numpy_support as nps - except ImportError: + except ModuleNotFoundError: msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ImportError(msg) + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err + mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): From fb40ba894addb3ab87226e80e631aaa6288dafa6 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:18:22 +0200 Subject: [PATCH 10/53] moved comment statement --- openmc/tallies.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index faa7c3efe0..bec0ddc6f7 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -467,7 +467,6 @@ class Tally(IDManagerMixin): ValueError: if no MeshFilter with appropriate mesh was found (SphericalMesh not supported) """ - # check that tally has a MeshFilter try: import vtk import vtk.util.numpy_support as nps @@ -477,6 +476,7 @@ class Tally(IDManagerMixin): except ImportError as err: raise err + # check that tally has a MeshFilter mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): From dacd3b75b855d8f28da8f7eedea4e6508c20fd80 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:20:25 +0200 Subject: [PATCH 11/53] moved import to voxels_to_vtk --- openmc/tallies.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index bec0ddc6f7..fbf2962b9b 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -467,15 +467,6 @@ class Tally(IDManagerMixin): ValueError: if no MeshFilter with appropriate mesh was found (SphericalMesh not supported) """ - try: - import vtk - import vtk.util.numpy_support as nps - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - # check that tally has a MeshFilter mesh = None for f in self.filters: @@ -3324,7 +3315,14 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): Returns: vtkStructuredGrid: a vtk object containing tally data on the appropriate grid """ - + try: + import vtk + import vtk.util.numpy_support as nps + except ModuleNotFoundError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err # TODO: should this be a method of Tally? vtk_grid = vtk.vtkStructuredGrid() From 9ec0a0a3df1a90783e1b0e60701f972275ffd17e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:42:32 +0200 Subject: [PATCH 12/53] refactoring --- openmc/tallies.py | 105 +++++++++++++++++++++++----------------------- 1 file changed, 53 insertions(+), 52 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index fbf2962b9b..a25d535e5f 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -464,61 +464,29 @@ class Tally(IDManagerMixin): filename (str): the filename (must end with .vtk) Raises: - ValueError: if no MeshFilter with appropriate mesh was found - (SphericalMesh not supported) + ValueError: if no MeshFilter was found """ + try: + import vtk + except ModuleNotFoundError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err + # check that tally has a MeshFilter mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): - if isinstance(f.mesh, (openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh)): - mesh = f.mesh - break + mesh = f.mesh + break if not mesh: raise ValueError( - "write_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" + "write_to_vtk requires a MeshFilter in the tally filters" ) - if isinstance(mesh, openmc.RegularMesh): - vtk_grid = voxels_to_vtk( - x_vals=np.linspace( - mesh.lower_left[0], - mesh.upper_right[0], - num=mesh.dimension[0] + 1, - ), - y_vals=np.linspace( - mesh.lower_left[1], - mesh.upper_right[1], - num=mesh.dimension[1] + 1, - ), - z_vals=np.linspace( - mesh.lower_left[2], - mesh.upper_right[2], - num=mesh.dimension[2] + 1, - ), - mean=self.mean, - std_dev=self.std_dev, - cylindrical=False, - ) - if isinstance(mesh, openmc.RectilinearMesh): - vtk_grid = voxels_to_vtk( - x_vals=mesh.x_grid, - y_vals=mesh.y_grid, - z_vals=mesh.z_grid, - mean=self.mean, - std_dev=self.std_dev, - cylindrical=False, - ) - elif isinstance(mesh, openmc.CylindricalMesh): - vtk_grid = voxels_to_vtk( - x_vals=mesh.r_grid, - y_vals=mesh.phi_grid, - z_vals=mesh.z_grid, - mean=self.mean, - std_dev=self.std_dev, - cylindrical=True, - ) + vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) # write the .vtk file writer = vtk.vtkStructuredGridWriter() @@ -3300,20 +3268,20 @@ class Tallies(cv.CheckedList): return cls(tallies) -def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): +def voxels_to_vtk(mesh, mean, std_dev): """Creates a vtk object from a list of X, Y, Z values and mean/std_dev data. Args: - x_vals (list): X values. - y_vals (list): Y values. - z_vals (list): Z values. + mesh (openmc.StructuredMesh): The tallied mesh (SphericalMesh is not yet supported). mean (np.array): the tally mean. std_dev (np.array): the tally standard deviation. - cylindrical (bool, optional): If set to True, cylindrical coordinates - (r, phi, z) are assumed. Defaults to True. Returns: vtkStructuredGrid: a vtk object containing tally data on the appropriate grid + + Raises: + ValueError: if mesh is not openmc.RegularMesh, openmc.RectilinearMesh or + openmc.CylindricalMesh """ try: import vtk @@ -3325,13 +3293,46 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): raise err # TODO: should this be a method of Tally? + system_of_coordinates = "cartesian" + + if isinstance(mesh, openmc.RegularMesh): + print('coucou') + x_vals = np.linspace( + mesh.lower_left[0], + mesh.upper_right[0], + num=mesh.dimension[0] + 1, + ) + y_vals = np.linspace( + mesh.lower_left[1], + mesh.upper_right[1], + num=mesh.dimension[1] + 1, + ) + z_vals = np.linspace( + mesh.lower_left[2], + mesh.upper_right[2], + num=mesh.dimension[2] + 1, + ) + elif isinstance(mesh, openmc.RectilinearMesh): + x_vals = mesh.x_grid + y_vals = mesh.y_grid + z_vals = mesh.z_grid + elif isinstance(mesh, openmc.CylindricalMesh): + x_vals = mesh.r_grid + y_vals = mesh.phi_grid + z_vals = mesh.z_grid + system_of_coordinates = "cylindrical" + else: + print(type(mesh)) + raise ValueError( + "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" + ) vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) # create points points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - if cylindrical: # transform points to cartesian coordinates + if system_of_coordinates == "cylindrical": # transform points to cartesian coordinates points_cartesian = np.copy(points) r, phi, z = points[:, 0], points[:, 1], points[:, 2] points_cartesian[:, 0] = r * np.cos(phi) From 11cee75e46ebc7e9756d81dbfae3452ff824cc10 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:12:28 +0200 Subject: [PATCH 13/53] put write_to_vtk method at the end --- openmc/tallies.py | 74 +++++++++++++++++++++++------------------------ 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index a25d535e5f..198f5581b2 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -457,43 +457,6 @@ class Tally(IDManagerMixin): self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) self._sparse = False - def write_to_vtk(self, filename): - """Writes the tally to a vtk file - - Args: - filename (str): the filename (must end with .vtk) - - Raises: - ValueError: if no MeshFilter was found - """ - try: - import vtk - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - - # check that tally has a MeshFilter - mesh = None - for f in self.filters: - if isinstance(f, openmc.MeshFilter): - mesh = f.mesh - break - - if not mesh: - raise ValueError( - "write_to_vtk requires a MeshFilter in the tally filters" - ) - - vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - def remove_score(self, score): """Remove a score from the tally @@ -3050,6 +3013,43 @@ class Tally(IDManagerMixin): new_tally.sparse = self.sparse return new_tally + def write_to_vtk(self, filename): + """Writes the tally to a vtk file + + Args: + filename (str): the filename (must end with .vtk) + + Raises: + ValueError: if no MeshFilter was found + """ + try: + import vtk + except ModuleNotFoundError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err + + # check that tally has a MeshFilter + mesh = None + for f in self.filters: + if isinstance(f, openmc.MeshFilter): + mesh = f.mesh + break + + if not mesh: + raise ValueError( + "write_to_vtk requires a MeshFilter in the tally filters" + ) + + vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + class Tallies(cv.CheckedList): """Collection of Tallies used for an OpenMC simulation. From bff8bc4b459030448749e587fb0fde2cf53d48f9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:16:40 +0200 Subject: [PATCH 14/53] added tests --- tests/unit_tests/test_tallies.py | 62 ++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index aeeb0612d5..b6fc6cc478 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -1,4 +1,8 @@ import numpy as np +import pytest +import vtk +from os.path import exists + import openmc @@ -38,3 +42,61 @@ def test_xml_roundtrip(run_in_tmpdir): assert new_tally.triggers[0].trigger_type == tally.triggers[0].trigger_type assert new_tally.triggers[0].threshold == tally.triggers[0].threshold assert new_tally.triggers[0].scores == tally.triggers[0].scores + + +cylinder_mesh = openmc.CylindricalMesh() +cylinder_mesh.r_grid = np.linspace(1, 2, num=30) +cylinder_mesh.phi_grid = np.linspace(0, np.pi / 2, num=50) +cylinder_mesh.z_grid = np.linspace(0, 1, num=30) + +regular_mesh = openmc.RegularMesh() +regular_mesh.lower_left = [0, 0, 0] +regular_mesh.upper_right = [1, 1, 1] +regular_mesh.dimension = [10, 5, 6] + +rectilinear_mesh = openmc.RectilinearMesh() +rectilinear_mesh.x_grid = np.linspace(0, 1) +rectilinear_mesh.y_grid = np.linspace(0, 1) +rectilinear_mesh.z_grid = np.linspace(0, 1) + + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) +def test_voxels_to_vtk(mesh): + vtk_grid = openmc.voxels_to_vtk(mesh, mean=None, std_dev=None) + assert isinstance(vtk_grid, vtk.vtkStructuredGrid) + + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) +def test_write_to_vtk(mesh, tmpdir): + # build + tally = openmc.Tally() + tally.filters = [openmc.MeshFilter(mesh)] + filename = tmpdir / "out.vtk" + # run + tally.write_to_vtk(filename) + # test + assert exists(filename) + + +def test_write_to_vtk_raises_error_when_no_meshfilter(): + # build + tally = openmc.Tally() + + # test + expected_err_msg = "write_to_vtk requires a MeshFilter in the tally filters" + with pytest.raises(ValueError, match=expected_err_msg): + tally.write_to_vtk("out.vtk") + + +def test_voxel_to_vtk_raises_error_with_wrong_mesh(): + # build + tally = openmc.Tally() + spherical_mesh = openmc.SphericalMesh() + spherical_mesh.r_grid = np.linspace(1, 2) + spherical_mesh.phi_grid = np.linspace(1, 2) + spherical_mesh.theta_grid = np.linspace(1, 2) + tally.filters = [openmc.MeshFilter(spherical_mesh)] + # test + expected_err_msg = "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" + with pytest.raises(ValueError, match=expected_err_msg): + tally.write_to_vtk("out.vtk") \ No newline at end of file From b61ec91d781caf8e7ca593a98c76d916a8efa94e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:50:28 +0200 Subject: [PATCH 15/53] removed print statement --- openmc/tallies.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 198f5581b2..88f747249d 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3296,7 +3296,6 @@ def voxels_to_vtk(mesh, mean, std_dev): system_of_coordinates = "cartesian" if isinstance(mesh, openmc.RegularMesh): - print('coucou') x_vals = np.linspace( mesh.lower_left[0], mesh.upper_right[0], From cac84d6e1f76b69d28da8674fcf81011170cf446 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:50:42 +0200 Subject: [PATCH 16/53] removd print statement --- openmc/tallies.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 88f747249d..020d4a3d9e 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3321,7 +3321,6 @@ def voxels_to_vtk(mesh, mean, std_dev): z_vals = mesh.z_grid system_of_coordinates = "cylindrical" else: - print(type(mesh)) raise ValueError( "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" ) From 7be28ef156a96fb4c525130d9c52e6966cf43539 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:08:58 +0200 Subject: [PATCH 17/53] added vtk_grid to meshes --- openmc/mesh.py | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 7f7b078795..dbdfb10b02 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -621,7 +621,45 @@ class RegularMesh(StructuredMesh): root_cell.fill = lattice return root_cell, cells + + def vtk_grid(self, filename=None): + import vtk + from vtk.util import numpy_support as nps + x_vals = np.linspace( + self.lower_left[0], + self.upper_right[0], + num=self.dimension[0] + 1, + ) + y_vals = np.linspace( + self.lower_left[1], + self.upper_right[1], + num=self.dimension[1] + 1, + ) + z_vals = np.linspace( + self.lower_left[2], + self.upper_right[2], + num=self.dimension[2] + 1, + ) + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + + # create points + pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid def Mesh(*args, **kwargs): warnings.warn("Mesh has been renamed RegularMesh. Future versions of " @@ -821,6 +859,33 @@ class RectilinearMesh(StructuredMesh): return element + def vtk_grid(self, filename=None): + import vtk + from vtk.util import numpy_support as nps + + x_vals = self.x_grid + y_vals = self.y_grid + z_vals = self.z_grid + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + + # create points + pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid + class CylindricalMesh(StructuredMesh): """A 3D cylindrical mesh @@ -1012,6 +1077,37 @@ class CylindricalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_p), V_z) + def vtk_grid(self, filename=None): + import vtk + from vtk.util import numpy_support as nps + + r_vals = self.r_grid + phi_vals = self.phi_grid + z_vals = self.z_grid + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(r_vals), len(phi_vals), len(z_vals)) + + # create points + pts_cylindrical = np.array([[r, phi, z] for z in z_vals for phi in phi_vals for r in r_vals]) + pts_cartesian = np.copy(pts_cylindrical) + r, phi, z = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] + pts_cartesian[:, 0] = r * np.cos(phi) + pts_cartesian[:, 1] = r * np.sin(phi) + pts_cartesian[:, 2] = z + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid class SphericalMesh(StructuredMesh): """A 3D spherical mesh @@ -1204,6 +1300,9 @@ class SphericalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_t), V_p) + def vtk_grid(self, filename=None): + raise NotImplementedError("vtk_grid not implemented for SphericalMesh") + class UnstructuredMesh(MeshBase): """A 3D unstructured mesh From 4ab6617bc4fe681bff68bce6dd4b22bf57a821d9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:09:18 +0200 Subject: [PATCH 18/53] write_to_vtk makes use of mesh.vtk_grid() --- openmc/tallies.py | 113 +++++----------------------------------------- 1 file changed, 12 insertions(+), 101 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 020d4a3d9e..e67365df66 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3042,7 +3042,18 @@ class Tally(IDManagerMixin): "write_to_vtk requires a MeshFilter in the tally filters" ) - vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) + vtk_grid = mesh.vtk_grid() + + # add mean and std dev data + mean_array = vtk.vtkDoubleArray() + mean_array.SetName("mean") + mean_array.SetArray(self.mean, self.mean.size, True) + vtk_grid.GetCellData().AddArray(mean_array) + + std_dev_array = vtk.vtkDoubleArray() + std_dev_array.SetName("std_dev") + std_dev_array.SetArray(self.std_dev, self.std_dev.size, True) + vtk_grid.GetCellData().AddArray(std_dev_array) # write the .vtk file writer = vtk.vtkStructuredGridWriter() @@ -3266,103 +3277,3 @@ class Tallies(cv.CheckedList): tallies.append(tally) return cls(tallies) - - -def voxels_to_vtk(mesh, mean, std_dev): - """Creates a vtk object from a list of X, Y, Z values and mean/std_dev data. - - Args: - mesh (openmc.StructuredMesh): The tallied mesh (SphericalMesh is not yet supported). - mean (np.array): the tally mean. - std_dev (np.array): the tally standard deviation. - - Returns: - vtkStructuredGrid: a vtk object containing tally data on the appropriate grid - - Raises: - ValueError: if mesh is not openmc.RegularMesh, openmc.RectilinearMesh or - openmc.CylindricalMesh - """ - try: - import vtk - import vtk.util.numpy_support as nps - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - # TODO: should this be a method of Tally? - - system_of_coordinates = "cartesian" - - if isinstance(mesh, openmc.RegularMesh): - x_vals = np.linspace( - mesh.lower_left[0], - mesh.upper_right[0], - num=mesh.dimension[0] + 1, - ) - y_vals = np.linspace( - mesh.lower_left[1], - mesh.upper_right[1], - num=mesh.dimension[1] + 1, - ) - z_vals = np.linspace( - mesh.lower_left[2], - mesh.upper_right[2], - num=mesh.dimension[2] + 1, - ) - elif isinstance(mesh, openmc.RectilinearMesh): - x_vals = mesh.x_grid - y_vals = mesh.y_grid - z_vals = mesh.z_grid - elif isinstance(mesh, openmc.CylindricalMesh): - x_vals = mesh.r_grid - y_vals = mesh.phi_grid - z_vals = mesh.z_grid - system_of_coordinates = "cylindrical" - else: - raise ValueError( - "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" - ) - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) - - # create points - points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - if system_of_coordinates == "cylindrical": # transform points to cartesian coordinates - points_cartesian = np.copy(points) - r, phi, z = points[:, 0], points[:, 1], points[:, 2] - points_cartesian[:, 0] = r * np.cos(phi) - points_cartesian[:, 1] = r * np.sin(phi) - points_cartesian[:, 2] = z - points = points_cartesian - - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # add mean and std dev data - mean_array = vtk.vtkDoubleArray() - mean_array.SetName("mean") - if mean is None: - mean = np.zeros( - (len(x_vals) - 1) - * (len(y_vals) - 1) - * (len(z_vals) - 1) - ) - mean_array.SetArray(mean, mean.size, True) - vtk_grid.GetCellData().AddArray(mean_array) - - std_dev_array = vtk.vtkDoubleArray() - std_dev_array.SetName("std_dev") - if std_dev is None: - std_dev = np.zeros( - (len(x_vals) - 1) - * (len(y_vals) - 1) - * (len(z_vals) - 1) - ) - std_dev_array.SetArray(std_dev, std_dev.size, True) - vtk_grid.GetCellData().AddArray(std_dev_array) - - return vtk_grid From 8958b74da73e47757147ab80b22b325926291142 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:09:32 +0200 Subject: [PATCH 19/53] adapted tests --- tests/unit_tests/test_tallies.py | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index b6fc6cc478..e49596f6da 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -43,6 +43,23 @@ def test_xml_roundtrip(run_in_tmpdir): assert new_tally.triggers[0].threshold == tally.triggers[0].threshold assert new_tally.triggers[0].scores == tally.triggers[0].scores +def run_dummy_sim(tally): + mat = openmc.Material() + mat.add_nuclide('Zr90', 1.0) + mat.set_density('g/cm3', 1.0) + + model = openmc.Model() + sph = openmc.Sphere(r=25.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sph) + model.geometry = openmc.Geometry([cell]) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 2 + model.settings.particles = 50 + + model.tallies = openmc.Tallies([tally]) + + model.run() cylinder_mesh = openmc.CylindricalMesh() cylinder_mesh.r_grid = np.linspace(1, 2, num=30) @@ -59,19 +76,13 @@ rectilinear_mesh.x_grid = np.linspace(0, 1) rectilinear_mesh.y_grid = np.linspace(0, 1) rectilinear_mesh.z_grid = np.linspace(0, 1) - -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) -def test_voxels_to_vtk(mesh): - vtk_grid = openmc.voxels_to_vtk(mesh, mean=None, std_dev=None) - assert isinstance(vtk_grid, vtk.vtkStructuredGrid) - - @pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) def test_write_to_vtk(mesh, tmpdir): # build tally = openmc.Tally() tally.filters = [openmc.MeshFilter(mesh)] filename = tmpdir / "out.vtk" + run_dummy_sim(tally) # run tally.write_to_vtk(filename) # test @@ -97,6 +108,6 @@ def test_voxel_to_vtk_raises_error_with_wrong_mesh(): spherical_mesh.theta_grid = np.linspace(1, 2) tally.filters = [openmc.MeshFilter(spherical_mesh)] # test - expected_err_msg = "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" - with pytest.raises(ValueError, match=expected_err_msg): + expected_err_msg = "vtk_grid not implemented for SphericalMesh" + with pytest.raises(NotImplementedError, match=expected_err_msg): tally.write_to_vtk("out.vtk") \ No newline at end of file From 2a78c06cf512e579eaeef007a63a2c90da9cb4ed Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:21:45 +0200 Subject: [PATCH 20/53] added docstrings --- openmc/mesh.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index dbdfb10b02..f6e3953282 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -623,6 +623,15 @@ class RegularMesh(StructuredMesh): return root_cell, cells def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ import vtk from vtk.util import numpy_support as nps @@ -860,6 +869,15 @@ class RectilinearMesh(StructuredMesh): return element def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ import vtk from vtk.util import numpy_support as nps @@ -1078,6 +1096,15 @@ class CylindricalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_p), V_z) def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ import vtk from vtk.util import numpy_support as nps @@ -1301,6 +1328,16 @@ class SphericalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_t), V_p) def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ + # FIXME raise NotImplementedError("vtk_grid not implemented for SphericalMesh") From 913bc99564b97490d90152aa68ed9a65838bdc3e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:33:17 +0200 Subject: [PATCH 21/53] minor refactore --- openmc/mesh.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index f6e3953282..d89a980ac6 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1108,15 +1108,12 @@ class CylindricalMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps - r_vals = self.r_grid - phi_vals = self.phi_grid - z_vals = self.z_grid vtk_grid = vtk.vtkStructuredGrid() - vtk_grid.SetDimensions(len(r_vals), len(phi_vals), len(z_vals)) + vtk_grid.SetDimensions(len(self.r_grid), len(self.phi_grid), len(self.z_grid)) # create points - pts_cylindrical = np.array([[r, phi, z] for z in z_vals for phi in phi_vals for r in r_vals]) + pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) pts_cartesian = np.copy(pts_cylindrical) r, phi, z = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] pts_cartesian[:, 0] = r * np.cos(phi) From aeb66509098b597a5cd53fca88d5e21096c691a0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:54:52 +0200 Subject: [PATCH 22/53] implemented spherical grid --- openmc/mesh.py | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index d89a980ac6..260ed5f43c 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1334,8 +1334,33 @@ class SphericalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object """ - # FIXME - raise NotImplementedError("vtk_grid not implemented for SphericalMesh") + import vtk + from vtk.util import numpy_support as nps + + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) + + # create points + pts_cylindrical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) + pts_cartesian = np.copy(pts_cylindrical) + r, theta, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] + pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) + pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) + pts_cartesian[:, 2] = r * np.cos(phi) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid class UnstructuredMesh(MeshBase): From 7de718a85c1ceaa0a6c485184829357d0f5c6372 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:56:32 +0200 Subject: [PATCH 23/53] adapted tests --- tests/unit_tests/test_tallies.py | 21 ++++++--------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index e49596f6da..e3f9a486da 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -76,7 +76,12 @@ rectilinear_mesh.x_grid = np.linspace(0, 1) rectilinear_mesh.y_grid = np.linspace(0, 1) rectilinear_mesh.z_grid = np.linspace(0, 1) -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) +spherical_mesh = openmc.SphericalMesh() +spherical_mesh.r_grid = np.linspace(1, 2) +spherical_mesh.phi_grid = np.linspace(1, 2) +spherical_mesh.theta_grid = np.linspace(1, 2) + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) def test_write_to_vtk(mesh, tmpdir): # build tally = openmc.Tally() @@ -97,17 +102,3 @@ def test_write_to_vtk_raises_error_when_no_meshfilter(): expected_err_msg = "write_to_vtk requires a MeshFilter in the tally filters" with pytest.raises(ValueError, match=expected_err_msg): tally.write_to_vtk("out.vtk") - - -def test_voxel_to_vtk_raises_error_with_wrong_mesh(): - # build - tally = openmc.Tally() - spherical_mesh = openmc.SphericalMesh() - spherical_mesh.r_grid = np.linspace(1, 2) - spherical_mesh.phi_grid = np.linspace(1, 2) - spherical_mesh.theta_grid = np.linspace(1, 2) - tally.filters = [openmc.MeshFilter(spherical_mesh)] - # test - expected_err_msg = "vtk_grid not implemented for SphericalMesh" - with pytest.raises(NotImplementedError, match=expected_err_msg): - tally.write_to_vtk("out.vtk") \ No newline at end of file From 0f23bf2109505f1ebb2de692552f096b25a5e93b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:21:59 +0200 Subject: [PATCH 24/53] write_data_to_vtk SphericalMesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 260ed5f43c..c780827ea6 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1324,7 +1324,7 @@ class SphericalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_t), V_p) - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -1337,6 +1337,19 @@ class SphericalMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + cv.check_type('label', label, str) + vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) @@ -1353,12 +1366,26 @@ class SphericalMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 15cc2a0804d3fa146bf6e337f4987ae634490513 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:23:55 +0200 Subject: [PATCH 25/53] fixed variable name --- openmc/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index c780827ea6..3dbc170d60 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1355,9 +1355,9 @@ class SphericalMesh(StructuredMesh): vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) # create points - pts_cylindrical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) - pts_cartesian = np.copy(pts_cylindrical) - r, theta, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] + pts_spherical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) + pts_cartesian = np.copy(pts_spherical) + r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2] pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) pts_cartesian[:, 2] = r * np.cos(phi) From 3b1d1756e4716e5e70f58a3142dc7efaed2d7913 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:26:29 +0200 Subject: [PATCH 26/53] write_data_to_vtk CylindricalMesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 3dbc170d60..029d714f0b 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1095,7 +1095,7 @@ class CylindricalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_p), V_z) - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -1108,6 +1108,19 @@ class CylindricalMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + cv.check_type('label', label, str) + vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(self.r_grid), len(self.phi_grid), len(self.z_grid)) @@ -1124,12 +1137,26 @@ class CylindricalMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 722b572dc314ba3dc5802da57a43e0d124323556 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:31:35 +0200 Subject: [PATCH 27/53] write_data_to_vtk RegularMesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 029d714f0b..179e9a1713 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -622,7 +622,7 @@ class RegularMesh(StructuredMesh): return root_cell, cells - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -635,6 +635,19 @@ class RegularMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1] * self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1] * self.dimension[2] + cv.check_type('label', label, str) + x_vals = np.linspace( self.lower_left[0], self.upper_right[0], @@ -661,12 +674,26 @@ class RegularMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 5ecac76b6960f635f375fd48f3cb567733073d3f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:34:05 +0200 Subject: [PATCH 28/53] write_data_to_vtk rectilinear mesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 179e9a1713..db2e1593c5 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -895,7 +895,7 @@ class RectilinearMesh(StructuredMesh): return element - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -908,6 +908,19 @@ class RectilinearMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + cv.check_type('label', label, str) + x_vals = self.x_grid y_vals = self.y_grid z_vals = self.z_grid @@ -922,12 +935,26 @@ class RectilinearMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 3233ffc759708b6df7a12c8234404897fa329e51 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:34:58 +0200 Subject: [PATCH 29/53] removed Tally.write_to_vtk --- openmc/tallies.py | 48 ----------------------------------------------- 1 file changed, 48 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index e67365df66..d0355f14ee 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3013,54 +3013,6 @@ class Tally(IDManagerMixin): new_tally.sparse = self.sparse return new_tally - def write_to_vtk(self, filename): - """Writes the tally to a vtk file - - Args: - filename (str): the filename (must end with .vtk) - - Raises: - ValueError: if no MeshFilter was found - """ - try: - import vtk - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - - # check that tally has a MeshFilter - mesh = None - for f in self.filters: - if isinstance(f, openmc.MeshFilter): - mesh = f.mesh - break - - if not mesh: - raise ValueError( - "write_to_vtk requires a MeshFilter in the tally filters" - ) - - vtk_grid = mesh.vtk_grid() - - # add mean and std dev data - mean_array = vtk.vtkDoubleArray() - mean_array.SetName("mean") - mean_array.SetArray(self.mean, self.mean.size, True) - vtk_grid.GetCellData().AddArray(mean_array) - - std_dev_array = vtk.vtkDoubleArray() - std_dev_array.SetName("std_dev") - std_dev_array.SetArray(self.std_dev, self.std_dev.size, True) - vtk_grid.GetCellData().AddArray(std_dev_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - class Tallies(cv.CheckedList): """Collection of Tallies used for an OpenMC simulation. From c58465214b96f7245979f792f4bf98201551e560 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:01:47 +0200 Subject: [PATCH 30/53] removed tally tests --- tests/unit_tests/test_tallies.py | 63 -------------------------------- 1 file changed, 63 deletions(-) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index e3f9a486da..bfff741a97 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -1,7 +1,4 @@ import numpy as np -import pytest -import vtk -from os.path import exists import openmc @@ -42,63 +39,3 @@ def test_xml_roundtrip(run_in_tmpdir): assert new_tally.triggers[0].trigger_type == tally.triggers[0].trigger_type assert new_tally.triggers[0].threshold == tally.triggers[0].threshold assert new_tally.triggers[0].scores == tally.triggers[0].scores - -def run_dummy_sim(tally): - mat = openmc.Material() - mat.add_nuclide('Zr90', 1.0) - mat.set_density('g/cm3', 1.0) - - model = openmc.Model() - sph = openmc.Sphere(r=25.0, boundary_type='vacuum') - cell = openmc.Cell(fill=mat, region=-sph) - model.geometry = openmc.Geometry([cell]) - - model.settings.run_mode = 'fixed source' - model.settings.batches = 2 - model.settings.particles = 50 - - model.tallies = openmc.Tallies([tally]) - - model.run() - -cylinder_mesh = openmc.CylindricalMesh() -cylinder_mesh.r_grid = np.linspace(1, 2, num=30) -cylinder_mesh.phi_grid = np.linspace(0, np.pi / 2, num=50) -cylinder_mesh.z_grid = np.linspace(0, 1, num=30) - -regular_mesh = openmc.RegularMesh() -regular_mesh.lower_left = [0, 0, 0] -regular_mesh.upper_right = [1, 1, 1] -regular_mesh.dimension = [10, 5, 6] - -rectilinear_mesh = openmc.RectilinearMesh() -rectilinear_mesh.x_grid = np.linspace(0, 1) -rectilinear_mesh.y_grid = np.linspace(0, 1) -rectilinear_mesh.z_grid = np.linspace(0, 1) - -spherical_mesh = openmc.SphericalMesh() -spherical_mesh.r_grid = np.linspace(1, 2) -spherical_mesh.phi_grid = np.linspace(1, 2) -spherical_mesh.theta_grid = np.linspace(1, 2) - -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) -def test_write_to_vtk(mesh, tmpdir): - # build - tally = openmc.Tally() - tally.filters = [openmc.MeshFilter(mesh)] - filename = tmpdir / "out.vtk" - run_dummy_sim(tally) - # run - tally.write_to_vtk(filename) - # test - assert exists(filename) - - -def test_write_to_vtk_raises_error_when_no_meshfilter(): - # build - tally = openmc.Tally() - - # test - expected_err_msg = "write_to_vtk requires a MeshFilter in the tally filters" - with pytest.raises(ValueError, match=expected_err_msg): - tally.write_to_vtk("out.vtk") From 15b98ea541fd27f6e02b5939bd482155cf44b56b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:02:12 +0200 Subject: [PATCH 31/53] added a test to meshes --- tests/unit_tests/test_mesh_to_vtk.py | 35 ++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 tests/unit_tests/test_mesh_to_vtk.py diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py new file mode 100644 index 0000000000..b5c8ca60a9 --- /dev/null +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -0,0 +1,35 @@ +import numpy as np +from os.path import exists +import pytest + +import openmc + + +regular_mesh = openmc.RegularMesh() +regular_mesh.lower_left = (0, 0, 0) +regular_mesh.upper_right = (1, 1, 1) +regular_mesh.dimension = [30, 20, 10] + +rectilinear_mesh = openmc.RectilinearMesh() +rectilinear_mesh.x_grid = np.linspace(1, 2, num=30) +rectilinear_mesh.y_grid = np.linspace(1, 2, num=30) +rectilinear_mesh.z_grid = np.linspace(1, 2, num=30) + +cylinder_mesh = openmc.CylindricalMesh() +cylinder_mesh.r_grid = np.linspace(1, 2, num=30) +cylinder_mesh.phi_grid = np.linspace(0, np.pi, num=50) +cylinder_mesh.z_grid = np.linspace(0, 1, num=30) + +spherical_mesh = openmc.SphericalMesh() +spherical_mesh.r_grid = np.linspace(1, 2, num=30) +spherical_mesh.phi_grid = np.linspace(0, np.pi, num=50) +spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30) + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) +def test_write_data_to_vtk(mesh, tmpdir): + filename = tmpdir / "out.vtk" + + data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) + + mesh.write_data_to_vtk(filename=filename, datasets={"label": data}) + assert exists(filename) \ No newline at end of file From 5174145c1ec9daacac8cfcd33e4b9212d9889fa1 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:21:36 +0200 Subject: [PATCH 32/53] added checks to test --- tests/unit_tests/test_mesh_to_vtk.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index b5c8ca60a9..e487b3668b 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,6 +1,8 @@ import numpy as np from os.path import exists import pytest +import vtk +from vtk.util import numpy_support as nps import openmc @@ -27,9 +29,30 @@ spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30) @pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) def test_write_data_to_vtk(mesh, tmpdir): + # BUILD filename = tmpdir / "out.vtk" data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) - mesh.write_data_to_vtk(filename=filename, datasets={"label": data}) - assert exists(filename) \ No newline at end of file + # RUN + mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) + + # TEST + assert exists(filename) + + # read file + reader = vtk.vtkStructuredGridReader() + reader.SetFileName(filename) + reader.Update() + + # check name of datasets + vtk_grid = reader.GetOutput() + array1 = vtk_grid.GetCellData().GetArray(0) + array2 = vtk_grid.GetCellData().GetArray(1) + + assert array1.GetName() == "label1" + assert array2.GetName() == "label2" + + # check size of datasets + assert nps.vtk_to_numpy(array1).size == data.size + assert nps.vtk_to_numpy(array2).size == data.size \ No newline at end of file From 5603e328b8650e874b2cf0235bba6fe792850129 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:24:44 +0200 Subject: [PATCH 33/53] docstrings --- openmc/mesh.py | 48 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 8 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index db2e1593c5..1372c5196e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -626,8 +626,16 @@ class RegularMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object @@ -899,8 +907,16 @@ class RectilinearMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object @@ -1153,8 +1169,16 @@ class CylindricalMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object @@ -1409,8 +1433,16 @@ class SphericalMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object From c77eb6097d6584820256ad497520273669299f9a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:18:31 +0200 Subject: [PATCH 34/53] removed check volume for structured meshes --- openmc/mesh.py | 36 ------------------------------------ 1 file changed, 36 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 1372c5196e..a63834941f 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -633,21 +633,12 @@ class RegularMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): @@ -914,21 +905,12 @@ class RectilinearMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): @@ -1176,21 +1158,12 @@ class CylindricalMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): @@ -1440,21 +1413,12 @@ class SphericalMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): From 1a7c0791bf92e46475343323d7ac0b939bcb9a5a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:19:12 +0200 Subject: [PATCH 35/53] z is unchanged --- openmc/mesh.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index a63834941f..1ccb6e86e0 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1179,10 +1179,9 @@ class CylindricalMesh(StructuredMesh): # create points pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) pts_cartesian = np.copy(pts_cylindrical) - r, phi, z = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] + r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1] pts_cartesian[:, 0] = r * np.cos(phi) pts_cartesian[:, 1] = r * np.sin(phi) - pts_cartesian[:, 2] = z vtkPts = vtk.vtkPoints() vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) From 9867a4340abbbd2fa093acb1fd04ed1745ab6513 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:19:56 +0200 Subject: [PATCH 36/53] [skip ci] use .num_mesh_cells Co-authored-by: Paul Romano --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index e487b3668b..b57235edda 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -32,7 +32,7 @@ def test_write_data_to_vtk(mesh, tmpdir): # BUILD filename = tmpdir / "out.vtk" - data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) + data = np.random.random(mesh.num_mesh_cells) # RUN mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) From 7e074fe354fc60887c8c4d6540c28e1a6a338ebc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:20:17 +0200 Subject: [PATCH 37/53] [skip ci] str(filename) Co-authored-by: Paul Romano --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 1372c5196e..e024710b48 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -699,7 +699,7 @@ class RegularMesh(StructuredMesh): # write the .vtk file writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) + writer.SetFileName(str(filename)) writer.SetInputData(vtk_grid) writer.Write() From 1155ffb739d8a80e5520c1013596864874a57cf5 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:21:15 +0200 Subject: [PATCH 38/53] added pathlib import --- tests/unit_tests/test_mesh_to_vtk.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index e487b3668b..cfb0407ca0 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,5 +1,6 @@ import numpy as np from os.path import exists +from pathlib import Path import pytest import vtk from vtk.util import numpy_support as nps From aeeb078344023ebb47ec7238b763c89a541a05fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:22:01 +0200 Subject: [PATCH 39/53] [skip ci] pathlib instead of os Co-authored-by: Paul Romano --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index b57235edda..3b7ccf9fa6 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -30,7 +30,7 @@ spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30) @pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) def test_write_data_to_vtk(mesh, tmpdir): # BUILD - filename = tmpdir / "out.vtk" + filename = Path(tmpdir) / "out.vtk" data = np.random.random(mesh.num_mesh_cells) From f99e93b5839f4166972fa8e83d93de6eccb5b619 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:22:26 +0200 Subject: [PATCH 40/53] [skip ci] is_file() instead of exists() Co-authored-by: Paul Romano --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 3b7ccf9fa6..31fbd4c882 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -38,7 +38,7 @@ def test_write_data_to_vtk(mesh, tmpdir): mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) # TEST - assert exists(filename) + assert filename.is_file() # read file reader = vtk.vtkStructuredGridReader() From 1266e433c8bd3f525e6b1ea9bafef174238be09c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:23:06 +0200 Subject: [PATCH 41/53] removed unused import --- tests/unit_tests/test_mesh_to_vtk.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 572834f84f..5da62cdef4 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,5 +1,4 @@ import numpy as np -from os.path import exists from pathlib import Path import pytest import vtk From 58ccb82cf50a5ba76e622f96969ba4db6752a885 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:26:27 +0200 Subject: [PATCH 42/53] back to .dimension --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 5da62cdef4..451c869c2e 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -32,7 +32,7 @@ def test_write_data_to_vtk(mesh, tmpdir): # BUILD filename = Path(tmpdir) / "out.vtk" - data = np.random.random(mesh.num_mesh_cells) + data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) # RUN mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) From 890c3ad0d58fa59532dc143e56c0100e4540fd0d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:26:46 +0200 Subject: [PATCH 43/53] pytest.importorskip --- tests/unit_tests/test_mesh_to_vtk.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 451c869c2e..f651f5acb5 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,7 +1,8 @@ import numpy as np from pathlib import Path import pytest -import vtk + +vtk = pytest.importorskip("vtk") from vtk.util import numpy_support as nps import openmc From 48a0d15b1aa6bbbfa2f3b3c85252a9d05c923028 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:36:54 +0200 Subject: [PATCH 44/53] added error raise + new test --- openmc/mesh.py | 53 ++++++++++++++++++++++------ tests/unit_tests/test_mesh_to_vtk.py | 17 ++++++++- 2 files changed, 59 insertions(+), 11 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index b883a4a096..ada69887ad 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -635,16 +635,22 @@ class RegularMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1] * self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1] * self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) x_vals = np.linspace( @@ -907,16 +913,22 @@ class RectilinearMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) x_vals = self.x_grid @@ -1160,16 +1172,22 @@ class CylindricalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) vtk_grid = vtk.vtkStructuredGrid() @@ -1414,16 +1432,23 @@ class SphericalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) + cv.check_type('label', label, str) vtk_grid = vtk.vtkStructuredGrid() @@ -1642,6 +1667,11 @@ class UnstructuredMesh(MeshBase): volume_normalization : bool Whether or not to normalize the data by the volume of the mesh elements + + Raises + ------ + RuntimeError + when the size of a dataset doesn't match the number of cells """ import vtk @@ -1658,11 +1688,14 @@ class UnstructuredMesh(MeshBase): " mesh information from a statepoint file.") # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.n_elements + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.n_elements + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) # create data arrays for the cells/points diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index f651f5acb5..9022feecd0 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -56,4 +56,19 @@ def test_write_data_to_vtk(mesh, tmpdir): # check size of datasets assert nps.vtk_to_numpy(array1).size == data.size - assert nps.vtk_to_numpy(array2).size == data.size \ No newline at end of file + assert nps.vtk_to_numpy(array2).size == data.size + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) +def test_write_data_to_vtk_size_mismatch(mesh): + """Checks that an error is raised when the size of the dataset + doesn't match the mesh number of cells + + Args: + mesh (openmc.StructuredMesh): the mesh to test + """ + right_size = mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2] + data = np.random.random(right_size + 1) + + expected_error_msg = "The size of the dataset label should be equal to the number of cells" + with pytest.raises(RuntimeError, match=expected_error_msg): + mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data}) \ No newline at end of file From 8857b769324f299ad410926b2542c7b7a24a895f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:50:28 +0200 Subject: [PATCH 45/53] refactored by adding StructuredMesh.write_data_to_vtk() --- openmc/mesh.py | 260 +++++++++++++------------------------------------ 1 file changed, 70 insertions(+), 190 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index ada69887ad..5083b215cf 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -193,6 +193,51 @@ class StructuredMesh(MeshBase): s1 = (slice(1, None),)*ndim + (slice(None),) return (vertices[s0] + vertices[s1]) / 2 + def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True): + import vtk + from vtk.util import numpy_support as nps + + # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) + else: + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) + cv.check_type('label', label, str) + + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(*self.dimension) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) + vtk_grid.SetPoints(vtkPts) + + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(str(filename)) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid class RegularMesh(StructuredMesh): """A regular Cartesian mesh in one, two, or three dimensions @@ -635,23 +680,7 @@ class RegularMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) x_vals = np.linspace( self.lower_left[0], @@ -668,39 +697,16 @@ class RegularMesh(StructuredMesh): self.upper_right[2], num=self.dimension[2] + 1, ) - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) # create points pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(str(filename)) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) def Mesh(*args, **kwargs): warnings.warn("Mesh has been renamed RegularMesh. Future versions of " @@ -913,60 +919,16 @@ class RectilinearMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) - - x_vals = self.x_grid - y_vals = self.y_grid - z_vals = self.z_grid - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) - # create points - pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + pts_cartesian = np.array([[x, y, z] for z in self.z_grid for y in self.y_grid for x in self.x_grid]) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) class CylindricalMesh(StructuredMesh): @@ -1172,28 +1134,7 @@ class CylindricalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) - - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(self.r_grid), len(self.phi_grid), len(self.z_grid)) - # create points pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) pts_cartesian = np.copy(pts_cylindrical) @@ -1201,32 +1142,12 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 0] = r * np.cos(phi) pts_cartesian[:, 1] = r * np.sin(phi) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) class SphericalMesh(StructuredMesh): """A 3D spherical mesh @@ -1432,28 +1353,7 @@ class SphericalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - - cv.check_type('label', label, str) - - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) # create points pts_spherical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) @@ -1463,32 +1363,12 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) pts_cartesian[:, 2] = r * np.cos(phi) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) class UnstructuredMesh(MeshBase): From 0c4de41f5076acbd93efaad6cdb970e26c890964 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:56:27 +0200 Subject: [PATCH 46/53] docstrings --- openmc/mesh.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 5083b215cf..449edb8d6a 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -194,6 +194,24 @@ class StructuredMesh(MeshBase): return (vertices[s0] + vertices[s1]) / 2 def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True): + """Creates a VTK object of the mesh + + Args: + points (list or np.array): List of (X,Y,Y) tuples. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ + import vtk from vtk.util import numpy_support as nps From ad73f03a0d712311b0d5ef80f0456d2c35b61dbe Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 15:02:20 +0200 Subject: [PATCH 47/53] converted docstrings to numpy style --- openmc/mesh.py | 129 +++++++++++++++++++++++++++++++------------------ 1 file changed, 81 insertions(+), 48 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 449edb8d6a..aa338f97e1 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -196,20 +196,29 @@ class StructuredMesh(MeshBase): def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - points (list or np.array): List of (X,Y,Y) tuples. - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + points : list or np.array + List of (X,Y,Y) tuples. + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells + Raises + ------ + RuntimeError + When the size of a dataset doesn't match the number of cells - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ import vtk @@ -688,16 +697,22 @@ class RegularMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ x_vals = np.linspace( @@ -927,16 +942,22 @@ class RectilinearMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ # create points pts_cartesian = np.array([[x, y, z] for z in self.z_grid for y in self.y_grid for x in self.x_grid]) @@ -1142,16 +1163,22 @@ class CylindricalMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ # create points pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) @@ -1361,16 +1388,22 @@ class SphericalMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ # create points From a2299063f5b4061a2eb2cd762004439fadbab57c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 15:02:52 +0200 Subject: [PATCH 48/53] converted docstrings to numpy style --- tests/unit_tests/test_mesh_to_vtk.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 9022feecd0..9247540a09 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -63,8 +63,10 @@ def test_write_data_to_vtk_size_mismatch(mesh): """Checks that an error is raised when the size of the dataset doesn't match the mesh number of cells - Args: - mesh (openmc.StructuredMesh): the mesh to test + Parameters + ---------- + mesh : openmc.StructuredMesh + The mesh to test """ right_size = mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2] data = np.random.random(right_size + 1) From 7cba2ad03012658481c6edbc04f070892e90dd3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 18:46:10 +0200 Subject: [PATCH 49/53] [skip ci] X Y Y typo Co-authored-by: Paul Romano --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index aa338f97e1..90f4c55a9a 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -199,7 +199,7 @@ class StructuredMesh(MeshBase): Parameters ---------- points : list or np.array - List of (X,Y,Y) tuples. + List of (X,Y,Z) tuples. filename : str Name of the VTK file to write. datasets : dict From 5180ab3f8a6031cfbcd708c04e7fd23f8936085a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 18:46:43 +0200 Subject: [PATCH 50/53] [skip ci] removed default value Co-authored-by: Paul Romano --- openmc/mesh.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 90f4c55a9a..73a49bfc19 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -208,7 +208,6 @@ class StructuredMesh(MeshBase): volume_normalization : bool, optional Whether or not to normalize the data by the volume of the mesh elements. - Defaults to True. Raises ------ From 94ef0c40abe6f4a208aaaccad36c17ff486dd9b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 18:47:24 +0200 Subject: [PATCH 51/53] [skip ci] refactore upperright and lowerleft Co-authored-by: Paul Romano --- openmc/mesh.py | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 73a49bfc19..7155b0ba04 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -714,21 +714,10 @@ class RegularMesh(StructuredMesh): the VTK object """ - x_vals = np.linspace( - self.lower_left[0], - self.upper_right[0], - num=self.dimension[0] + 1, - ) - y_vals = np.linspace( - self.lower_left[1], - self.upper_right[1], - num=self.dimension[1] + 1, - ) - z_vals = np.linspace( - self.lower_left[2], - self.upper_right[2], - num=self.dimension[2] + 1, - ) + ll, ur = self.lower_left, self.upper_right + x_vals = np.linspace(ll[0], ur[0], num=self.dimension[0] + 1) + y_vals = np.linspace(ll[1], ur[1], num=self.dimension[1] + 1) + z_vals = np.linspace(ll[2], ur[2], num=self.dimension[2] + 1) # create points pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) From 96c89f138daa97f78cb996449a32a88ca58754ff Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 19:01:39 +0200 Subject: [PATCH 52/53] multiple lines --- openmc/mesh.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 7155b0ba04..209002e42a 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1169,7 +1169,14 @@ class CylindricalMesh(StructuredMesh): the VTK object """ # create points - pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) + pts_cylindrical = np.array( + [ + [r, phi, z] + for z in self.z_grid + for phi in self.phi_grid + for r in self.r_grid + ] + ) pts_cartesian = np.copy(pts_cylindrical) r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1] pts_cartesian[:, 0] = r * np.cos(phi) @@ -1395,7 +1402,14 @@ class SphericalMesh(StructuredMesh): """ # create points - pts_spherical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) + pts_spherical = np.array( + [ + [r, theta, phi] + for phi in self.phi_grid + for theta in self.theta_grid + for r in self.r_grid + ] + ) pts_cartesian = np.copy(pts_spherical) r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2] pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) From fabe3e607c92891ec740aefa632eba010428eef2 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 19:02:25 +0200 Subject: [PATCH 53/53] pathlib --- openmc/mesh.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 209002e42a..2a53c8a91e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -698,7 +698,7 @@ class RegularMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -932,7 +932,7 @@ class RectilinearMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -1153,7 +1153,7 @@ class CylindricalMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -1385,7 +1385,7 @@ class SphericalMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -1592,7 +1592,7 @@ class UnstructuredMesh(MeshBase): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels