diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 326c16182b..cc09709d79 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -783,6 +783,11 @@ public: int set_grid(); + enum class Orientation { + y, //!< Flat side of lattice parallel to y-axis + x //!< Flat side of lattice parallel to x-axis + }; + // Data members int num_rings_; double pitch_; @@ -791,15 +796,9 @@ public: Direction r_; Direction q_dual_; Direction r_dual_; - -private: - enum class Orientation { - y, //!< Flat side of lattice parallel to y-axis - x //!< Flat side of lattice parallel to x-axis - }; - Orientation orientation_ {Orientation::y}; +private: StructuredMesh::MeshDistance find_z_crossing( const Position& r, const Direction& u, double l, int shell) const; diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 19e6f74d7a..f9223e1484 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -18,7 +18,7 @@ from ..mesh import MeshMaterialVolumes __all__ = [ 'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh', - 'SphericalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes' + 'SphericalMesh', 'HexagonalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes' ] @@ -108,6 +108,11 @@ _dll.openmc_spherical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double), _dll.openmc_spherical_mesh_set_grid.restype = c_int _dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler +_dll.openmc_hexagonal_mesh_get_grid.argtypes = [c_int32, + POINTER(POINTER(c_double)), POINTER(c_int), POINTER(c_double), + POINTER(POINTER(c_double)), POINTER(c_double), POINTER(c_char_p)] +_dll.openmc_hexagonal_mesh_get_grid.restype = c_int +_dll.openmc_hexagonal_mesh_get_grid.errcheck = _error_handler class Mesh(_FortranObjectWithID): """Base class to represent mesh objects @@ -727,6 +732,91 @@ class SphericalMesh(Mesh): ntheta, phi_grid, nphi) +class HexagonalMesh(Mesh): + """HexagonalMesh stored internally. + + This class exposes a mesh that is stored internally in the OpenMC + library. To obtain a view of a mesh with a given ID, use the + :data:`openmc.lib.meshes` mapping. + + Parameters + ---------- + index : int + Index in the `meshes` array. + + Attributes + ---------- + id : int + ID of the mesh + z_grid : numpy.ndarray + 1-D array of mesh boundary points along the z-axis. + pitch : float + Radial pitch of the hexagonal mesh in cm. + num_rings : int + Number of radial ring positions in the xy-plane + orientation : {'x', 'y'} + The orientation of the lattice. The 'x' orientation means that each + lattice element has two faces that are perpendicular to the x-axis, + while the 'y' orientation means that each lattice element has two faces + that are perpendicular to the y-axis. By default, the orientation is + 'y'. + origin : numpy.ndarray + 1-D array of length 3 the (x,y,z) origin of the mesh in + cartesian coordinates + n_elements : int + Total number of mesh elements. + volumes : numpy.ndarray + Volume of each mesh element in [cm^3] + bounding_box : openmc.BoundingBox + Axis-aligned bounding box of the mesh + + """ + mesh_type = 'hexagonal' + + def __init__(self, uid=None, new=True, index=None): + super().__init__(uid, new, index) + + @property + def n_elements(self): + z_grid, nr, *_ = self._get_parameters() + return (z_grid.size-1)*(3*nr*(nr-1)+1) + + @property + def num_rings(self): + return self._get_parameters()[1] + + @property + def pitch(self): + return self._get_parameters()[2] + + @property + def orientation(self): + return self._get_parameters()[4] + + @property + def origin(self): + return self._get_parameters()[3] + + @property + def z_grid(self): + return self._get_parameters()[0] + + def _get_parameters(self): + gz = POINTER(c_double)() + nz = c_int() + nr = c_int() + orig = POINTER(c_double)() + orient = c_char() + p = c_double() + # Call C API to get grid parameters + _dll.openmc_hexagonal_mesh_get_grid(self._index, gz, nz, nr, orig, p, orient) + + # Convert grid parameters to Numpy arrays + grid_z = as_array(gz, (nz.value,)) + origin = as_array(orig, (3,)) + + return (grid_z, nr, pitch, origin, orientation.decode()) + class UnstructuredMesh(Mesh): pass @@ -736,6 +826,7 @@ _MESH_TYPE_MAP = { 'rectilinear': RectilinearMesh, 'cylindrical': CylindricalMesh, 'spherical': SphericalMesh, + 'hexagonal': HexagonalMesh, 'unstructured': UnstructuredMesh } diff --git a/openmc/mesh.py b/openmc/mesh.py index 7feac1d195..5899a54ab4 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2638,10 +2638,18 @@ class HexagonalMesh(StructuredMesh): cv.check_type('mesh radial pitch', pitch, Real) cv.check_greater_than('mesh radial pitch', pitch, 0.0) self._pitch = pitch + + @property + def n_elements(self): + return (self.z_grid.size - 1)*(3*self.num_rings*(self.num_rings - 1) + 1) + + @property + def dimension(self): + return (self.n_elements,) @property def n_dimension(self): - return 4 + return 4 @property def lower_left(self): @@ -2695,7 +2703,7 @@ class HexagonalMesh(StructuredMesh): idx = [] for j in range(self.z_grid.size-1): idx.append((0, 0, 0, j)) - for rad in range(1, self.num_rings+1): + for rad in range(1, self.num_rings): for i in range(rad): idx.append((i, rad-i, -rad, j)) idx.append((rad, -i, i-rad, j)) diff --git a/src/mesh.cpp b/src/mesh.cpp index 5808d7c21e..724393b755 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2490,7 +2490,7 @@ std::string HexagonalMesh::surface_label(int surface) const int HexagonalMesh::n_bins() const { - return (1 + 3 * (num_rings_ + 1) * num_rings_) * (grid_.size() - 1); + return (1 + 3 * num_rings_ * (num_rings_ - 1)) * (grid_.size() - 1); } bool HexagonalMesh::valid_index(const MeshIndex& ijk, int k) const @@ -2507,7 +2507,7 @@ int HexagonalMesh::get_bin_from_indices(const MeshIndex& ijk) const int q = ijk[0]; int r = ijk[1]; int k = ijk[2]; - int hexes = 3 * num_rings_ * (num_rings_ + 1) + 1; + int hexes = 3 * num_rings_ * (num_rings_ - 1) + 1; int bin = (k - 1) * (grid_.size() - 1) * hexes; int rad = std::max({std::abs(r), std::abs(q), std::abs(r + q)}); if (rad == 0) @@ -2576,7 +2576,7 @@ StructuredMesh::MeshIndex HexagonalMesh::get_indices( StructuredMesh::MeshIndex HexagonalMesh::get_indices_from_bin(int bin) const { MeshIndex ijk = {0, 0, 0}; - int hexes = 3 * num_rings_ * (num_rings_ + 1) + 1; + int hexes = 3 * num_rings_ * (num_rings_ - 1) + 1; ijk[2] = static_cast(std::floor(bin / hexes)) + 1; int sp_idx = bin % hexes; if (sp_idx == 0) { @@ -3277,6 +3277,28 @@ extern "C" int openmc_spherical_mesh_set_grid(int32_t index, index, grid_x, nx, grid_y, ny, grid_z, nz); } +//! Get the hexagonal mesh grid +extern "C" int openmc_hexagonal_mesh_get_grid(int32_t index, double** grid_z, + int* nz, int* nr, double** origin, double* pitch, const char** orient) +{ + if (int err = check_mesh_type(index)) + return err; + HexagonalMesh* m = dynamic_cast(model::meshes[index].get()); + + if (m->grid_.empty()) { + set_errmsg("Mesh parameters have not been set."); + return OPENMC_E_ALLOCATE; + } + + *grid_z = m->grid_.data(); + *nz = m->grid_.size(); + *nr = m->num_rings_; + *origin = &m->origin_.x; + *pitch = m->pitch_; + *orient = (m->orientation_ == HexagonalMesh::Orientation::y) ? "y" : "x"; + return 0; +} + #ifdef OPENMC_DAGMC_ENABLED const std::string MOABMesh::mesh_lib_type = "moab";