bugfix

Author: GuySten

include/openmc/mesh.h

@@ -783,6 +783,11 @@ class HexagonalMesh : public PeriodicStructuredMesh {
 
   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 @@ class HexagonalMesh : public PeriodicStructuredMesh {
   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;
 

openmc/lib/mesh.py

@@ -18,7 +18,7 @@
 
 __all__ = [
     'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh',
-    'SphericalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes'
+    'SphericalMesh', 'HexagonalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes'
 ]
 
 
@@ -108,6 +108,11 @@
 _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 @@ def set_grid(self, r_grid, theta_grid, phi_grid):
                                               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 @@ class UnstructuredMesh(Mesh):
     'rectilinear': RectilinearMesh,
     'cylindrical': CylindricalMesh,
     'spherical': SphericalMesh,
+    'hexagonal': HexagonalMesh,
     'unstructured': UnstructuredMesh
 }
 

openmc/mesh.py

@@ -2638,10 +2638,18 @@ def pitch(self, pitch):
         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 @@ def indices(self):
         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))

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<int>(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<HexagonalMesh>(index))
+    return err;
+  HexagonalMesh* m = dynamic_cast<HexagonalMesh*>(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";