from abc import ABC, abstractmethod
from collections.abc import Iterable
from collections import OrderedDict
from math import pi
from numbers import Real, Integral
from pathlib import Path
import warnings
from xml.etree import ElementTree as ET
import numpy as np
import openmc.checkvalue as cv
import openmc
from ._xml import get_text
from .mixin import IDManagerMixin
from .surface import _BOUNDARY_TYPES
class MeshBase(IDManagerMixin, ABC):
"""A mesh that partitions geometry for tallying purposes.
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
"""
next_id = 1
used_ids = set()
def __init__(self, mesh_id=None, name=''):
# Initialize Mesh class attributes
self.id = mesh_id
self.name = name
@property
def name(self):
return self._name
@name.setter
def name(self, name):
if name is not None:
cv.check_type(f'name for mesh ID="{self._id}"', name, str)
self._name = name
else:
self._name = ''
def __repr__(self):
string = type(self).__name__ + '\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
return string
def _volume_dim_check(self):
if self.n_dimension != 3 or \
any([d == 0 for d in self.dimension]):
raise RuntimeError(f'Mesh {self.id} is not 3D. '
'Volumes cannot be provided.')
@classmethod
def from_hdf5(cls, group):
"""Create mesh from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
Returns
-------
openmc.MeshBase
Instance of a MeshBase subclass
"""
mesh_type = group['type'][()].decode()
if mesh_type == 'regular':
return RegularMesh.from_hdf5(group)
elif mesh_type == 'rectilinear':
return RectilinearMesh.from_hdf5(group)
elif mesh_type == 'cylindrical':
return CylindricalMesh.from_hdf5(group)
elif mesh_type == 'spherical':
return SphericalMesh.from_hdf5(group)
elif mesh_type == 'unstructured':
return UnstructuredMesh.from_hdf5(group)
else:
raise ValueError('Unrecognized mesh type: "' + mesh_type + '"')
@classmethod
def from_xml_element(cls, elem):
"""Generates a mesh from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.MeshBase
an openmc mesh object
"""
mesh_type = get_text(elem, 'type')
if mesh_type == 'regular' or mesh_type is None:
return RegularMesh.from_xml_element(elem)
elif mesh_type == 'rectilinear':
return RectilinearMesh.from_xml_element(elem)
elif mesh_type == 'cylindrical':
return CylindricalMesh.from_xml_element(elem)
elif mesh_type == 'spherical':
return SphericalMesh.from_xml_element(elem)
elif mesh_type == 'unstructured':
return UnstructuredMesh.from_xml_element(elem)
else:
raise ValueError(f'Unrecognized mesh type "{mesh_type}" found.')
class StructuredMesh(MeshBase):
"""A base class for structured mesh functionality
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
@property
@abstractmethod
def dimension(self):
pass
@property
@abstractmethod
def n_dimension(self):
pass
@property
@abstractmethod
def _grids(self):
pass
@property
def vertices(self):
"""Return coordinates of mesh vertices.
Returns
-------
vertices : numpy.ndarray
Returns a numpy.ndarray representing the coordinates of the mesh
vertices with a shape equal to (ndim, dim1 + 1, ..., dimn + 1). Can be
unpacked along the first dimension with xx, yy, zz = mesh.vertices.
"""
return np.stack(np.meshgrid(*self._grids, indexing='ij'), axis=0)
@property
def centroids(self):
"""Return coordinates of mesh element centroids.
Returns
-------
centroids : numpy.ndarray
Returns a numpy.ndarray representing the mesh element centroid
coordinates with a shape equal to (ndim, dim1, ..., dimn). Can be
unpacked along the first dimension with xx, yy, zz = mesh.centroids.
"""
ndim = self.n_dimension
vertices = self.vertices
s0 = (slice(None),) + (slice(0, -1),)*ndim
s1 = (slice(None),) + (slice(1, None),)*ndim
return (vertices[s0] + vertices[s1]) / 2
@property
def num_mesh_cells(self):
return np.prod(self.dimension)
def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True):
"""Creates a VTK object of the mesh
Parameters
----------
points : list or np.array
List of (X,Y,Z) 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.
Raises
------
ValueError
When the size of a dataset doesn't match the number of mesh cells
Returns
-------
vtk.vtkStructuredGrid
the VTK object
"""
import vtk
from vtk.util import numpy_support as nps
# check that the data sets are appropriately sized
for label, dataset in datasets.items():
errmsg = (
f"The size of the dataset '{label}' ({dataset.size}) should be"
f" equal to the number of mesh cells ({self.num_mesh_cells})"
)
if isinstance(dataset, np.ndarray):
if not dataset.size == self.num_mesh_cells:
raise ValueError(errmsg)
else:
if len(dataset) == self.num_mesh_cells:
raise ValueError(errmsg)
cv.check_type('label', label, str)
vtk_grid = vtk.vtkStructuredGrid()
vtk_grid.SetDimensions(*[dim + 1 for dim in 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
# maintain a list of the datasets as added
# to the VTK arrays to ensure they persist
# in memory until the file is written
datasets_out = []
for label, dataset in datasets.items():
dataset = np.asarray(dataset).flatten()
datasets_out.append(dataset)
if volume_normalization:
dataset /= self.volumes.T.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
[docs]class RegularMesh(StructuredMesh):
"""A regular Cartesian mesh in one, two, or three dimensions
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
dimension : Iterable of int
The number of mesh cells in each direction.
n_dimension : int
Number of mesh dimensions.
lower_left : Iterable of float
The lower-left corner of the structured mesh. If only two coordinate
are given, it is assumed that the mesh is an x-y mesh.
upper_right : Iterable of float
The upper-right corner of the structured mesh. If only two coordinate
are given, it is assumed that the mesh is an x-y mesh.
width : Iterable of float
The width of mesh cells in each direction.
indices : Iterable of tuple
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
(2, 1, 1), ...]
"""
def __init__(self, mesh_id=None, name=''):
super().__init__(mesh_id, name)
self._dimension = None
self._lower_left = None
self._upper_right = None
self._width = None
@property
def dimension(self):
return tuple(self._dimension)
@property
def n_dimension(self):
if self._dimension is not None:
return len(self._dimension)
else:
return None
@property
def lower_left(self):
return self._lower_left
@property
def upper_right(self):
if self._upper_right is not None:
return self._upper_right
elif self._width is not None:
if self._lower_left is not None and self._dimension is not None:
ls = self._lower_left
ws = self._width
dims = self._dimension
return [l + w * d for l, w, d in zip(ls, ws, dims)]
@property
def width(self):
if self._width is not None:
return self._width
elif self._upper_right is not None:
if self._lower_left is not None and self._dimension is not None:
us = self._upper_right
ls = self._lower_left
dims = self._dimension
return [(u - l) / d for u, l, d in zip(us, ls, dims)]
@property
def volumes(self):
"""Return Volumes for every mesh cell
Returns
-------
volumes : numpy.ndarray
Volumes
"""
self._volume_dim_check()
return np.full(self.dimension, np.prod(self.width))
@property
def total_volume(self):
return np.prod(self.dimension) * np.prod(self.width)
@property
def indices(self):
ndim = len(self._dimension)
if ndim == 3:
nx, ny, nz = self.dimension
return ((x, y, z)
for z in range(1, nz + 1)
for y in range(1, ny + 1)
for x in range(1, nx + 1))
elif ndim == 2:
nx, ny = self.dimension
return ((x, y)
for y in range(1, ny + 1)
for x in range(1, nx + 1))
else:
nx, = self.dimension
return ((x,) for x in range(1, nx + 1))
@property
def _grids(self):
ndim = len(self._dimension)
if ndim == 3:
x0, y0, z0 = self.lower_left
x1, y1, z1 = self.upper_right
nx, ny, nz = self.dimension
xarr = np.linspace(x0, x1, nx + 1)
yarr = np.linspace(y0, y1, ny + 1)
zarr = np.linspace(z0, z1, nz + 1)
return (xarr, yarr, zarr)
elif ndim == 2:
x0, y0 = self.lower_left
x1, y1 = self.upper_right
nx, ny = self.dimension
xarr = np.linspace(x0, x1, nx + 1)
yarr = np.linspace(y0, y1, ny + 1)
return (xarr, yarr)
else:
nx, = self.dimension
x0, = self.lower_left
x1, = self.upper_right
return (np.linspace(x0, x1, nx + 1),)
@dimension.setter
def dimension(self, dimension):
cv.check_type('mesh dimension', dimension, Iterable, Integral)
cv.check_length('mesh dimension', dimension, 1, 3)
self._dimension = dimension
@lower_left.setter
def lower_left(self, lower_left):
cv.check_type('mesh lower_left', lower_left, Iterable, Real)
cv.check_length('mesh lower_left', lower_left, 1, 3)
self._lower_left = lower_left
if self.upper_right is not None and any(np.isclose(self.upper_right, lower_left)):
raise ValueError("Mesh cannot have zero thickness in any dimension")
@upper_right.setter
def upper_right(self, upper_right):
cv.check_type('mesh upper_right', upper_right, Iterable, Real)
cv.check_length('mesh upper_right', upper_right, 1, 3)
self._upper_right = upper_right
if self._width is not None:
self._width = None
warnings.warn("Unsetting width attribute.")
if self.lower_left is not None and any(np.isclose(self.lower_left, upper_right)):
raise ValueError("Mesh cannot have zero thickness in any dimension")
@width.setter
def width(self, width):
cv.check_type('mesh width', width, Iterable, Real)
cv.check_length('mesh width', width, 1, 3)
self._width = width
if self._upper_right is not None:
self._upper_right = None
warnings.warn("Unsetting upper_right attribute.")
def __repr__(self):
string = super().__repr__()
string += '{0: <16}{1}{2}\n'.format('\tDimensions', '=\t', self.n_dimension)
string += '{0: <16}{1}{2}\n'.format('\tVoxels', '=\t', self._dimension)
string += '{0: <16}{1}{2}\n'.format('\tLower left', '=\t', self._lower_left)
string += '{0: <16}{1}{2}\n'.format('\tUpper Right', '=\t', self._upper_right)
string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._width)
return string
[docs] @classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.dimension = group['dimension'][()]
mesh.lower_left = group['lower_left'][()]
if 'width' in group:
mesh.width = group['width'][()]
elif 'upper_right' in group:
mesh.upper_right = group['upper_right'][()]
else:
raise IOError('Invalid mesh: must have one of "upper_right" or "width"')
return mesh
[docs] @classmethod
def from_rect_lattice(cls, lattice, division=1, mesh_id=None, name=''):
"""Create mesh from an existing rectangular lattice
Parameters
----------
lattice : openmc.RectLattice
Rectangular lattice used as a template for this mesh
division : int
Number of mesh cells per lattice cell.
If not specified, there will be 1 mesh cell per lattice cell.
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Returns
-------
openmc.RegularMesh
RegularMesh instance
"""
cv.check_type('rectangular lattice', lattice, openmc.RectLattice)
shape = np.array(lattice.shape)
width = lattice.pitch*shape
mesh = cls(mesh_id, name)
mesh.lower_left = lattice.lower_left
mesh.upper_right = lattice.lower_left + width
mesh.dimension = shape*division
return mesh
[docs] @classmethod
def from_domain(
cls,
domain,
dimension=(10, 10, 10),
mesh_id=None,
name=''
):
"""Create mesh from an existing openmc cell, region, universe or
geometry by making use of the objects bounding box property.
Parameters
----------
domain : {openmc.Cell, openmc.Region, openmc.Universe, openmc.Geometry}
The object passed in will be used as a template for this mesh. The
bounding box of the property of the object passed will be used to
set the lower_left and upper_right of the mesh instance
dimension : Iterable of int
The number of mesh cells in each direction.
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Returns
-------
openmc.RegularMesh
RegularMesh instance
"""
cv.check_type(
"domain",
domain,
(openmc.Cell, openmc.Region, openmc.Universe, openmc.Geometry),
)
mesh = cls(mesh_id, name)
mesh.lower_left = domain.bounding_box[0]
mesh.upper_right = domain.bounding_box[1]
mesh.dimension = dimension
return mesh
[docs] def to_xml_element(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
if self._dimension is not None:
subelement = ET.SubElement(element, "dimension")
subelement.text = ' '.join(map(str, self._dimension))
subelement = ET.SubElement(element, "lower_left")
subelement.text = ' '.join(map(str, self._lower_left))
if self._upper_right is not None:
subelement = ET.SubElement(element, "upper_right")
subelement.text = ' '.join(map(str, self._upper_right))
if self._width is not None:
subelement = ET.SubElement(element, "width")
subelement.text = ' '.join(map(str, self._width))
return element
[docs] @classmethod
def from_xml_element(cls, elem):
"""Generate mesh from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.Mesh
Mesh generated from XML element
"""
mesh_id = int(get_text(elem, 'id'))
mesh = cls(mesh_id)
mesh_type = get_text(elem, 'type')
if mesh_type is not None:
mesh.type = mesh_type
dimension = get_text(elem, 'dimension')
if dimension is not None:
mesh.dimension = [int(x) for x in dimension.split()]
lower_left = get_text(elem, 'lower_left')
if lower_left is not None:
mesh.lower_left = [float(x) for x in lower_left.split()]
upper_right = get_text(elem, 'upper_right')
if upper_right is not None:
mesh.upper_right = [float(x) for x in upper_right.split()]
width = get_text(elem, 'width')
if width is not None:
mesh.width = [float(x) for x in width.split()]
return mesh
[docs] def build_cells(self, bc=None):
"""Generates a lattice of universes with the same dimensionality
as the mesh object. The individual cells/universes produced
will not have material definitions applied and so downstream code
will have to apply that information.
Parameters
----------
bc : iterable of {'reflective', 'periodic', 'transmission', 'vacuum', or 'white'}
Boundary conditions for each of the four faces of a rectangle
(if applying to a 2D mesh) or six faces of a parallelepiped
(if applying to a 3D mesh) provided in the following order:
[x min, x max, y min, y max, z min, z max]. 2-D cells do not
contain the z min and z max entries. Defaults to 'reflective' for
all faces.
Returns
-------
root_cell : openmc.Cell
The cell containing the lattice representing the mesh geometry;
this cell is a single parallelepiped with boundaries matching
the outermost mesh boundary with the boundary conditions from bc
applied.
cells : iterable of openmc.Cell
The list of cells within each lattice position mimicking the mesh
geometry.
"""
if bc is None:
bc = ['reflective'] * 6
if len(bc) not in (4, 6):
raise ValueError('Boundary condition must be of length 4 or 6')
for entry in bc:
cv.check_value('bc', entry, _BOUNDARY_TYPES)
n_dim = self.n_dimension
# Build the cell which will contain the lattice
xplanes = [openmc.XPlane(self.lower_left[0], boundary_type=bc[0]),
openmc.XPlane(self.upper_right[0], boundary_type=bc[1])]
if n_dim == 1:
yplanes = [openmc.YPlane(-1e10, boundary_type='reflective'),
openmc.YPlane(1e10, boundary_type='reflective')]
else:
yplanes = [openmc.YPlane(self.lower_left[1], boundary_type=bc[2]),
openmc.YPlane(self.upper_right[1], boundary_type=bc[3])]
if n_dim <= 2:
# Would prefer to have the z ranges be the max supported float, but
# these values are apparently different between python and Fortran.
# Choosing a safe and sane default.
# Values of +/-1e10 are used here as there seems to be an
# inconsistency between what numpy uses as the max float and what
# Fortran expects for a real(8), so this avoids code complication
# and achieves the same goal.
zplanes = [openmc.ZPlane(-1e10, boundary_type='reflective'),
openmc.ZPlane(1e10, boundary_type='reflective')]
else:
zplanes = [openmc.ZPlane(self.lower_left[2], boundary_type=bc[4]),
openmc.ZPlane(self.upper_right[2], boundary_type=bc[5])]
root_cell = openmc.Cell()
root_cell.region = ((+xplanes[0] & -xplanes[1]) &
(+yplanes[0] & -yplanes[1]) &
(+zplanes[0] & -zplanes[1]))
# Build the universes which will be used for each of the (i,j,k)
# locations within the mesh.
# We will concurrently build cells to assign to these universes
cells = []
universes = []
for _ in self.indices:
cells.append(openmc.Cell())
universes.append(openmc.Universe())
universes[-1].add_cell(cells[-1])
lattice = openmc.RectLattice()
lattice.lower_left = self.lower_left
# Assign the universe and rotate to match the indexing expected for
# the lattice
if n_dim == 1:
universe_array = np.array([universes])
elif n_dim == 2:
universe_array = np.empty(self.dimension[::-1],
dtype=openmc.Universe)
i = 0
for y in range(self.dimension[1] - 1, -1, -1):
for x in range(self.dimension[0]):
universe_array[y][x] = universes[i]
i += 1
else:
universe_array = np.empty(self.dimension[::-1],
dtype=openmc.Universe)
i = 0
for z in range(self.dimension[2]):
for y in range(self.dimension[1] - 1, -1, -1):
for x in range(self.dimension[0]):
universe_array[z][y][x] = universes[i]
i += 1
lattice.universes = universe_array
if self.width is not None:
lattice.pitch = self.width
else:
dx = ((self.upper_right[0] - self.lower_left[0]) /
self.dimension[0])
if n_dim == 1:
lattice.pitch = [dx]
elif n_dim == 2:
dy = ((self.upper_right[1] - self.lower_left[1]) /
self.dimension[1])
lattice.pitch = [dx, dy]
else:
dy = ((self.upper_right[1] - self.lower_left[1]) /
self.dimension[1])
dz = ((self.upper_right[2] - self.lower_left[2]) /
self.dimension[2])
lattice.pitch = [dx, dy, dz]
# Fill Cell with the Lattice
root_cell.fill = lattice
return root_cell, cells
[docs] def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
"""Creates a VTK object of the mesh
Parameters
----------
filename : str or pathlib.Path
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
"""
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])
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 "
"OpenMC will not accept the name Mesh.")
return RegularMesh(*args, **kwargs)
[docs]class RectilinearMesh(StructuredMesh):
"""A 3D rectilinear Cartesian mesh
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
dimension : Iterable of int
The number of mesh cells in each direction.
n_dimension : int
Number of mesh dimensions (always 3 for a RectilinearMesh).
x_grid : numpy.ndarray
1-D array of mesh boundary points along the x-axis.
y_grid : numpy.ndarray
1-D array of mesh boundary points along the y-axis.
z_grid : numpy.ndarray
1-D array of mesh boundary points along the z-axis.
indices : Iterable of tuple
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
(2, 1, 1), ...]
"""
def __init__(self, mesh_id=None, name=''):
super().__init__(mesh_id, name)
self._x_grid = None
self._y_grid = None
self._z_grid = None
@property
def dimension(self):
return (len(self.x_grid) - 1,
len(self.y_grid) - 1,
len(self.z_grid) - 1)
@property
def n_dimension(self):
return 3
@property
def x_grid(self):
return self._x_grid
@property
def y_grid(self):
return self._y_grid
@property
def z_grid(self):
return self._z_grid
@property
def _grids(self):
return (self.x_grid, self.y_grid, self.z_grid)
@property
def volumes(self):
"""Return Volumes for every mesh cell
Returns
-------
volumes : numpy.ndarray
Volumes
"""
self._volume_dim_check()
V_x = np.diff(self.x_grid)
V_y = np.diff(self.y_grid)
V_z = np.diff(self.z_grid)
return np.multiply.outer(np.outer(V_x, V_y), V_z)
@property
def total_volume(self):
return np.sum(self.volumes)
@property
def indices(self):
nx = len(self.x_grid) - 1
ny = len(self.y_grid) - 1
nz = len(self.z_grid) - 1
return ((x, y, z)
for z in range(1, nz + 1)
for y in range(1, ny + 1)
for x in range(1, nx + 1))
@x_grid.setter
def x_grid(self, grid):
cv.check_type('mesh x_grid', grid, Iterable, Real)
self._x_grid = np.asarray(grid)
@y_grid.setter
def y_grid(self, grid):
cv.check_type('mesh y_grid', grid, Iterable, Real)
self._y_grid = np.asarray(grid)
@z_grid.setter
def z_grid(self, grid):
cv.check_type('mesh z_grid', grid, Iterable, Real)
self._z_grid = np.asarray(grid)
def __repr__(self):
fmt = '{0: <16}{1}{2}\n'
string = super().__repr__()
string += fmt.format('\tDimensions', '=\t', self.n_dimension)
x_grid_str = str(self._x_grid) if self._x_grid is None else len(self._x_grid)
string += fmt.format('\tN X pnts:', '=\t', x_grid_str)
if self._x_grid is not None:
string += fmt.format('\tX Min:', '=\t', self._x_grid[0])
string += fmt.format('\tX Max:', '=\t', self._x_grid[-1])
y_grid_str = str(self._y_grid) if self._y_grid is None else len(self._y_grid)
string += fmt.format('\tN Y pnts:', '=\t', y_grid_str)
if self._y_grid is not None:
string += fmt.format('\tY Min:', '=\t', self._y_grid[0])
string += fmt.format('\tY Max:', '=\t', self._y_grid[-1])
z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
if self._z_grid is not None:
string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
return string
[docs] @classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.x_grid = group['x_grid'][()]
mesh.y_grid = group['y_grid'][()]
mesh.z_grid = group['z_grid'][()]
return mesh
[docs] @classmethod
def from_xml_element(cls, elem):
"""Generate a rectilinear mesh from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.RectilinearMesh
Rectilinear mesh object
"""
id = int(get_text(elem, 'id'))
mesh = cls(id)
mesh.x_grid = [float(x) for x in get_text(elem, 'x_grid').split()]
mesh.y_grid = [float(y) for y in get_text(elem, 'y_grid').split()]
mesh.z_grid = [float(z) for z in get_text(elem, 'z_grid').split()]
return mesh
[docs] def to_xml_element(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
element.set("type", "rectilinear")
subelement = ET.SubElement(element, "x_grid")
subelement.text = ' '.join(map(str, self.x_grid))
subelement = ET.SubElement(element, "y_grid")
subelement.text = ' '.join(map(str, self.y_grid))
subelement = ET.SubElement(element, "z_grid")
subelement.text = ' '.join(map(str, self.z_grid))
return element
[docs] def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
"""Creates a VTK object of the mesh
Parameters
----------
filename : str or pathlib.Path
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
"""
# create points
pts_cartesian = self.vertices.T.reshape(-1, 3)
return super().write_data_to_vtk(
points=pts_cartesian,
filename=filename,
datasets=datasets,
volume_normalization=volume_normalization
)
[docs]class CylindricalMesh(StructuredMesh):
"""A 3D cylindrical mesh
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
dimension : Iterable of int
The number of mesh cells in each direction.
n_dimension : int
Number of mesh dimensions (always 3 for a CylindricalMesh).
r_grid : numpy.ndarray
1-D array of mesh boundary points along the r-axis.
Requirement is r >= 0.
phi_grid : numpy.ndarray
1-D array of mesh boundary points along the phi-axis in radians.
The default value is [0, 2Ï€], i.e. the full phi range.
z_grid : numpy.ndarray
1-D array of mesh boundary points along the z-axis.
origin : numpy.ndarray
1-D array of length 3 the (x,y,z) origin of the mesh in
cartesian coordinates
indices : Iterable of tuple
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
(2, 1, 1), ...]
"""
def __init__(self, mesh_id=None, name=''):
super().__init__(mesh_id, name)
self._r_grid = None
self._phi_grid = [0.0, 2*pi]
self._z_grid = None
self.origin = (0., 0., 0.)
@property
def dimension(self):
return (len(self.r_grid) - 1,
len(self.phi_grid) - 1,
len(self.z_grid) - 1)
@property
def n_dimension(self):
return 3
@property
def origin(self):
return self._origin
@property
def r_grid(self):
return self._r_grid
@property
def phi_grid(self):
return self._phi_grid
@property
def z_grid(self):
return self._z_grid
@property
def _grids(self):
return (self.r_grid, self.phi_grid, self.z_grid)
@property
def indices(self):
nr, np, nz = self.dimension
np = len(self.phi_grid) - 1
nz = len(self.z_grid) - 1
return ((r, p, z)
for z in range(1, nz + 1)
for p in range(1, np + 1)
for r in range(1, nr + 1))
@origin.setter
def origin(self, coords):
cv.check_type('mesh origin', coords, Iterable, Real)
cv.check_length("mesh origin", coords, 3)
self._origin = np.asarray(coords)
@r_grid.setter
def r_grid(self, grid):
cv.check_type('mesh r_grid', grid, Iterable, Real)
self._r_grid = np.asarray(grid)
@phi_grid.setter
def phi_grid(self, grid):
cv.check_type('mesh phi_grid', grid, Iterable, Real)
self._phi_grid = np.asarray(grid)
@z_grid.setter
def z_grid(self, grid):
cv.check_type('mesh z_grid', grid, Iterable, Real)
self._z_grid = np.asarray(grid)
def __repr__(self):
fmt = '{0: <16}{1}{2}\n'
string = super().__repr__()
string += fmt.format('\tDimensions', '=\t', self.n_dimension)
string += fmt.format('\tOrigin', '=\t', self.origin)
r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
if self._r_grid is not None:
string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
if self._phi_grid is not None:
string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
if self._z_grid is not None:
string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
return string
[docs] @classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.r_grid = group['r_grid'][()]
mesh.phi_grid = group['phi_grid'][()]
mesh.z_grid = group['z_grid'][()]
if 'origin' in group:
mesh.origin = group['origin'][()]
return mesh
[docs] @classmethod
def from_domain(
cls,
domain,
dimension=(10, 10, 10),
mesh_id=None,
phi_grid_bounds=(0.0, 2*pi),
name=''
):
"""Creates a regular CylindricalMesh from an existing openmc domain.
Parameters
----------
domain : openmc.Cell or openmc.Region or openmc.Universe or openmc.Geometry
The object passed in will be used as a template for this mesh. The
bounding box of the property of the object passed will be used to
set the r_grid, z_grid ranges.
dimension : Iterable of int
The number of equally spaced mesh cells in each direction (r_grid,
phi_grid, z_grid)
mesh_id : int
Unique identifier for the mesh
phi_grid_bounds : numpy.ndarray
Mesh bounds points along the phi-axis in radians. The default value
is (0, 2Ï€), i.e., the full phi range.
name : str
Name of the mesh
Returns
-------
openmc.RegularMesh
RegularMesh instance
"""
cv.check_type(
"domain",
domain,
(openmc.Cell, openmc.Region, openmc.Universe, openmc.Geometry),
)
mesh = cls(mesh_id, name)
# loaded once to avoid reading h5m file repeatedly
cached_bb = domain.bounding_box
max_bounding_box_radius = max(
[
cached_bb[0][0],
cached_bb[0][1],
cached_bb[1][0],
cached_bb[1][1],
]
)
mesh.r_grid = np.linspace(
0,
max_bounding_box_radius,
num=dimension[0]+1
)
mesh.phi_grid = np.linspace(
phi_grid_bounds[0],
phi_grid_bounds[1],
num=dimension[1]+1
)
mesh.z_grid = np.linspace(
cached_bb[0][2],
cached_bb[1][2],
num=dimension[2]+1
)
return mesh
[docs] def to_xml_element(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
element.set("type", "cylindrical")
subelement = ET.SubElement(element, "r_grid")
subelement.text = ' '.join(map(str, self.r_grid))
subelement = ET.SubElement(element, "phi_grid")
subelement.text = ' '.join(map(str, self.phi_grid))
subelement = ET.SubElement(element, "z_grid")
subelement.text = ' '.join(map(str, self.z_grid))
subelement = ET.SubElement(element, "origin")
subelement.text = ' '.join(map(str, self.origin))
return element
[docs] @classmethod
def from_xml_element(cls, elem):
"""Generate a cylindrical mesh from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.CylindricalMesh
Cylindrical mesh object
"""
mesh_id = int(get_text(elem, 'id'))
mesh = cls(mesh_id)
mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()]
mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()]
return mesh
@property
def volumes(self):
"""Return Volumes for every mesh cell
Returns
-------
volumes : Iterable of float
Volumes
"""
self._volume_dim_check()
V_r = np.diff(np.asarray(self.r_grid)**2 / 2)
V_p = np.diff(self.phi_grid)
V_z = np.diff(self.z_grid)
return np.multiply.outer(np.outer(V_r, V_p), V_z)
[docs] def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
"""Creates a VTK object of the mesh
Parameters
----------
filename : str or pathlib.Path
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
"""
# create points
pts_cylindrical = self.vertices.T.reshape(-1, 3)
pts_cartesian = np.copy(pts_cylindrical)
r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1]
pts_cartesian[:, 0] = r * np.cos(phi) + self.origin[0]
pts_cartesian[:, 1] = r * np.sin(phi) + self.origin[1]
pts_cartesian[:, 2] += self.origin[2]
return super().write_data_to_vtk(
points=pts_cartesian,
filename=filename,
datasets=datasets,
volume_normalization=volume_normalization
)
[docs]class SphericalMesh(StructuredMesh):
"""A 3D spherical mesh
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
dimension : Iterable of int
The number of mesh cells in each direction.
n_dimension : int
Number of mesh dimensions (always 3 for a SphericalMesh).
r_grid : numpy.ndarray
1-D array of mesh boundary points along the r-axis.
Requirement is r >= 0.
theta_grid : numpy.ndarray
1-D array of mesh boundary points along the theta-axis in radians.
The default value is [0, π], i.e. the full theta range.
phi_grid : numpy.ndarray
1-D array of mesh boundary points along the phi-axis in radians.
The default value is [0, 2Ï€], i.e. the full phi range.
origin : numpy.ndarray
1-D array of length 3 the (x,y,z) origin of the mesh in
cartesian coordinates
indices : Iterable of tuple
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
(2, 1, 1), ...]
"""
def __init__(self, mesh_id=None, name=''):
super().__init__(mesh_id, name)
self._r_grid = None
self._theta_grid = [0, pi]
self._phi_grid = [0, 2*pi]
self.origin = (0., 0., 0.)
@property
def dimension(self):
return (len(self.r_grid) - 1,
len(self.theta_grid) - 1,
len(self.phi_grid) - 1)
@property
def n_dimension(self):
return 3
@property
def origin(self):
return self._origin
@property
def r_grid(self):
return self._r_grid
@property
def theta_grid(self):
return self._theta_grid
@property
def phi_grid(self):
return self._phi_grid
@property
def _grids(self):
return (self.r_grid, self.theta_grid, self.phi_grid)
@property
def indices(self):
nr, nt, np = self.dimension
nt = len(self.theta_grid) - 1
np = len(self.phi_grid) - 1
return ((r, t, p)
for p in range(1, np + 1)
for t in range(1, nt + 1)
for r in range(1, nr + 1))
@origin.setter
def origin(self, coords):
cv.check_type('mesh origin', coords, Iterable, Real)
cv.check_length("mesh origin", coords, 3)
self._origin = np.asarray(coords)
@r_grid.setter
def r_grid(self, grid):
cv.check_type('mesh r_grid', grid, Iterable, Real)
self._r_grid = np.asarray(grid)
@theta_grid.setter
def theta_grid(self, grid):
cv.check_type('mesh theta_grid', grid, Iterable, Real)
self._theta_grid = np.asarray(grid)
@phi_grid.setter
def phi_grid(self, grid):
cv.check_type('mesh phi_grid', grid, Iterable, Real)
self._phi_grid = np.asarray(grid)
def __repr__(self):
fmt = '{0: <16}{1}{2}\n'
string = super().__repr__()
string += fmt.format('\tDimensions', '=\t', self.n_dimension)
string += fmt.format('\tOrigin', '=\t', self.origin)
r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
if self._r_grid is not None:
string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
theta_grid_str = str(self._theta_grid) if self._theta_grid is None else len(self._theta_grid)
string += fmt.format('\tN Theta pnts:', '=\t', theta_grid_str)
if self._theta_grid is not None:
string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0])
string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1])
phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
if self._phi_grid is not None:
string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
return string
[docs] @classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.r_grid = group['r_grid'][()]
mesh.theta_grid = group['theta_grid'][()]
mesh.phi_grid = group['phi_grid'][()]
if 'origin' in group:
mesh.origin = group['origin'][()]
return mesh
[docs] def to_xml_element(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
element.set("type", "spherical")
subelement = ET.SubElement(element, "r_grid")
subelement.text = ' '.join(map(str, self.r_grid))
subelement = ET.SubElement(element, "theta_grid")
subelement.text = ' '.join(map(str, self.theta_grid))
subelement = ET.SubElement(element, "phi_grid")
subelement.text = ' '.join(map(str, self.phi_grid))
subelement = ET.SubElement(element, "origin")
subelement.text = ' '.join(map(str, self.origin))
return element
[docs] @classmethod
def from_xml_element(cls, elem):
"""Generate a spherical mesh from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.SphericalMesh
Spherical mesh object
"""
mesh_id = int(get_text(elem, 'id'))
mesh = cls(mesh_id)
mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()]
mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()]
return mesh
@property
def volumes(self):
"""Return Volumes for every mesh cell
Returns
-------
volumes : Iterable of float
Volumes
"""
self._volume_dim_check()
V_r = np.diff(np.asarray(self.r_grid)**3 / 3)
V_t = np.diff(-np.cos(self.theta_grid))
V_p = np.diff(self.phi_grid)
return np.multiply.outer(np.outer(V_r, V_t), V_p)
[docs] def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
"""Creates a VTK object of the mesh
Parameters
----------
filename : str or pathlib.Path
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
"""
# create points
pts_spherical = self.vertices.T.reshape(-1, 3)
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) + self.origin[0]
pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) + self.origin[1]
pts_cartesian[:, 2] = r * np.cos(phi) + self.origin[2]
return super().write_data_to_vtk(
points=pts_cartesian,
filename=filename,
datasets=datasets,
volume_normalization=volume_normalization
)
[docs]class UnstructuredMesh(MeshBase):
"""A 3D unstructured mesh
.. versionadded:: 0.12
.. versionchanged:: 0.12.2
Support for libMesh unstructured meshes was added.
Parameters
----------
filename : str or pathlib.Path
Location of the unstructured mesh file
library : {'moab', 'libmesh'}
Mesh library used for the unstructured mesh tally
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
length_multiplier: float
Constant multiplier to apply to mesh coordinates
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
filename : str
Name of the file containing the unstructured mesh
length_multiplier: float
Multiplicative factor to apply to mesh coordinates
library : {'moab', 'libmesh'}
Mesh library used for the unstructured mesh tally
output : bool
Indicates whether or not automatic tally output should
be generated for this mesh
volumes : Iterable of float
Volumes of the unstructured mesh elements
centroids : numpy.ndarray
Centroids of the mesh elements with array shape (n_elements, 3)
vertices : numpy.ndarray
Coordinates of the mesh vertices with array shape (n_elements, 3)
.. versionadded:: 0.13.1
connectivity : numpy.ndarray
Connectivity of the elements with array shape (n_elements, 8)
.. versionadded:: 0.13.1
element_types : Iterable of integers
Mesh element types
.. versionadded:: 0.13.1
total_volume : float
Volume of the unstructured mesh in total
"""
_UNSUPPORTED_ELEM = -1
_LINEAR_TET = 0
_LINEAR_HEX = 1
def __init__(self, filename, library, mesh_id=None, name='',
length_multiplier=1.0):
super().__init__(mesh_id, name)
self.filename = filename
self._volumes = None
self._n_elements = None
self._conectivity = None
self._vertices = None
self.library = library
self._output = False
self.length_multiplier = length_multiplier
@property
def filename(self):
return self._filename
@filename.setter
def filename(self, filename):
cv.check_type('Unstructured Mesh filename', filename, (str, Path))
self._filename = filename
@property
def library(self):
return self._library
@library.setter
def library(self, lib):
cv.check_value('Unstructured mesh library', lib, ('moab', 'libmesh'))
self._library = lib
@property
def size(self):
return self._size
@size.setter
def size(self, size):
cv.check_type("Unstructured mesh size", size, Integral)
self._size = size
@property
def output(self):
return self._output
@output.setter
def output(self, val):
cv.check_type("Unstructured mesh output value", val, bool)
self._output = val
@property
def volumes(self):
"""Return Volumes for every mesh cell if
populated by a StatePoint file
Returns
-------
volumes : numpy.ndarray
Volumes
"""
return self._volumes
@volumes.setter
def volumes(self, volumes):
cv.check_type("Unstructured mesh volumes", volumes, Iterable, Real)
self._volumes = volumes
@property
def total_volume(self):
return np.sum(self.volumes)
@property
def vertices(self):
return self._vertices
@property
def connectivity(self):
return self._connectivity
@property
def element_types(self):
return self._element_types
@property
def centroids(self):
return np.array([self.centroid(i) for i in range(self.n_elements)])
@property
def n_elements(self):
if self._n_elements is None:
raise RuntimeError("No information about this mesh has "
"been loaded from a statepoint file.")
return self._n_elements
@n_elements.setter
def n_elements(self, val):
cv.check_type('Number of elements', val, Integral)
self._n_elements = val
@property
def length_multiplier(self):
return self._length_multiplier
@length_multiplier.setter
def length_multiplier(self, length_multiplier):
cv.check_type("Unstructured mesh length multiplier",
length_multiplier,
Real)
self._length_multiplier = length_multiplier
@property
def dimension(self):
return (self.n_elements,)
@property
def n_dimension(self):
return 3
def __repr__(self):
string = super().__repr__()
string += '{: <16}=\t{}\n'.format('\tFilename', self.filename)
string += '{: <16}=\t{}\n'.format('\tMesh Library', self.library)
if self.length_multiplier != 1.0:
string += '{: <16}=\t{}\n'.format('\tLength multiplier',
self.length_multiplier)
return string
[docs] def centroid(self, bin):
"""Return the vertex averaged centroid of an element
Parameters
----------
bin : int
Bin ID for the returned centroid
Returns
-------
numpy.ndarray
x, y, z values of the element centroid
"""
conn = self.connectivity[bin]
# remove invalid connectivity values
conn = conn[conn >= 0]
coords = self.vertices[conn]
return coords.mean(axis=0)
[docs] def write_vtk_mesh(self, **kwargs):
"""Map data to unstructured VTK mesh elements.
.. deprecated:: 0.13
Use :func:`UnstructuredMesh.write_data_to_vtk` instead.
Parameters
----------
filename : str or pathlib.Path
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
Whether or not to normalize the data by the
volume of the mesh elements
"""
warnings.warn(
"The 'UnstructuredMesh.write_vtk_mesh' method has been renamed "
"to 'write_data_to_vtk' and will be removed in a future version "
" of OpenMC.", FutureWarning
)
self.write_data_to_vtk(**kwargs)
[docs] def write_data_to_vtk(self, filename=None, datasets=None, volume_normalization=True):
"""Map data to unstructured VTK mesh elements.
Parameters
----------
filename : str or pathlib.Path
Name of the VTK file to write
datasets : dict
Dictionary whose keys are the data labels
and values are numpy appropriately sized arrays
of the data
volume_normalization : bool
Whether or not to normalize the data by the
volume of the mesh elements
"""
import vtk
from vtk.util import numpy_support as nps
if self.connectivity is None or self.vertices is None:
raise RuntimeError('This mesh has not been '
'loaded from a statepoint file.')
if filename is None:
filename = f'mesh_{self.id}.vtk'
writer = vtk.vtkUnstructuredGridWriter()
writer.SetFileName(str(filename))
grid = vtk.vtkUnstructuredGrid()
vtk_pnts = vtk.vtkPoints()
vtk_pnts.SetData(nps.numpy_to_vtk(self.vertices))
grid.SetPoints(vtk_pnts)
n_skipped = 0
elems = []
for elem_type, conn in zip(self.element_types, self.connectivity):
if elem_type == self._LINEAR_TET:
elem = vtk.vtkTetra()
elif elem_type == self._LINEAR_HEX:
elem = vtk.vtkHexahedron()
elif elem_type == self._UNSUPPORTED_ELEM:
n_skipped += 1
else:
raise RuntimeError(f'Invalid element type {elem_type} found')
for i, c in enumerate(conn):
if c == -1:
break
elem.GetPointIds().SetId(i, c)
elems.append(elem)
if n_skipped > 0:
warnings.warn(f'{n_skipped} elements were not written because '
'they are not of type linear tet/hex')
for elem in elems:
grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())
# check that datasets are the correct size
datasets_out = []
if datasets is not None:
for name, data in datasets.items():
if data.shape != self.dimension:
raise ValueError(f'Cannot apply dataset "{name}" with '
f'shape {data.shape} to mesh {self.id} '
f'with dimensions {self.dimension}')
if volume_normalization:
for name, data in datasets.items():
if np.issubdtype(data.dtype, np.integer):
warnings.warn(f'Integer data set "{name}" will '
'not be volume-normalized.')
continue
data /= self.volumes
# add data to the mesh
for name, data in datasets.items():
datasets_out.append(data)
arr = vtk.vtkDoubleArray()
arr.SetName(name)
arr.SetNumberOfTuples(data.size)
for i in range(data.size):
arr.SetTuple1(i, data.flat[i])
grid.GetCellData().AddArray(arr)
writer.SetInputData(grid)
writer.Write()
[docs] @classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
filename = group['filename'][()].decode()
library = group['library'][()].decode()
mesh = cls(filename, library, mesh_id=mesh_id)
vol_data = group['volumes'][()]
mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
mesh.n_elements = mesh.volumes.size
vertices = group['vertices'][()]
mesh._vertices = vertices.reshape((-1, 3))
connectivity = group['connectivity'][()]
mesh._connectivity = connectivity.reshape((-1, 8))
mesh._element_types = group['element_types'][()]
if 'length_multiplier' in group:
mesh.length_multiplier = group['length_multiplier'][()]
return mesh
[docs] def to_xml_element(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
element.set("type", "unstructured")
element.set("library", self._library)
subelement = ET.SubElement(element, "filename")
subelement.text = str(self.filename)
if self._length_multiplier != 1.0:
element.set("length_multiplier", str(self.length_multiplier))
return element
[docs] @classmethod
def from_xml_element(cls, elem):
"""Generate unstructured mesh object from XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.UnstructuredMesh
UnstructuredMesh generated from an XML element
"""
mesh_id = int(get_text(elem, 'id'))
filename = get_text(elem, 'filename')
library = get_text(elem, 'library')
length_multiplier = float(get_text(elem, 'length_multiplier', 1.0))
return cls(filename, library, mesh_id, '', length_multiplier)
def _read_meshes(elem):
"""Generate dictionary of meshes from a given XML node
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
dict
A dictionary with mesh IDs as keys and openmc.MeshBase
instanaces as values
"""
out = dict()
for mesh_elem in elem.findall('mesh'):
mesh = MeshBase.from_xml_element(mesh_elem)
out[mesh.id] = mesh
return out