Minimal implementation of hexagonal mesh

docs/source/io_formats/tallies.rst

@@ -362,7 +362,7 @@ attributes/sub-elements:
 
   :type:
     The type of mesh. This can be either "regular", "rectilinear",
-    "cylindrical", "spherical", or "unstructured".
+    "cylindrical", "spherical", "hexagonal", or "unstructured".
 
   :dimension:
     The number of mesh cells in each direction. (For regular mesh only.)
@@ -379,6 +379,9 @@ attributes/sub-elements:
     The upper-right corner of the structured mesh. If only two coordinates are
     given, it is assumed that the mesh is an x-y mesh. (For regular mesh only.)
 
+  :pitch:
+    The mesh radial pitch. (For hexagonal mesh only.)
+
   :width:
     The width of mesh cells in each direction. (For regular mesh only.)
 
@@ -389,7 +392,7 @@ attributes/sub-elements:
     The mesh divisions along the y-axis. (For rectilinear mesh only.)
 
   :z_grid:
-    The mesh divisions along the z-axis. (For rectilinear and cylindrical meshes only.)
+    The mesh divisions along the z-axis. (For rectilinear, cylindrical and hexagonal meshes only.)
 
   :r_grid:
     The mesh divisions along the r-axis. (For cylindrical and spherical meshes only.)
@@ -400,8 +403,14 @@ attributes/sub-elements:
   :theta_grid:
     The mesh divisions along the theta-axis. (For spherical mesh only.)
 
+  :num_rings:
+    The number of hexagonal rings. (For hexagonal mesh only.)
+
+  :orientation:
+    The orientation of the hexagonal mesh, either "x" or "y". (For hexagonal mesh only.)
+
   :origin:
-    The origin in cartesian coordinates. (For cylindrical and spherical meshes only.)
+    The origin in cartesian coordinates. (For cylindrical, spherical and hexagonal meshes only.)
 
   :library:
     The mesh library used to represent an unstructured mesh. This can be either

docs/source/pythonapi/base.rst

@@ -169,6 +169,7 @@ Meshes
    openmc.CylindricalMesh
    openmc.SphericalMesh
    openmc.UnstructuredMesh
+   openmc.HexagonalMesh
 
 Geometry Plotting
 -----------------

include/openmc/mesh.h

@@ -302,10 +302,10 @@ class StructuredMesh : public Mesh {
 
   struct MeshDistance {
     MeshDistance() = default;
-    MeshDistance(int _index, bool _max_surface, double _distance)
-      : next_index {_index}, max_surface {_max_surface}, distance {_distance}
+    MeshDistance(MeshIndex _offset, bool _max_surface, double _distance)
+      : offset {_offset}, max_surface {_max_surface}, distance {_distance}
     {}
-    int next_index {-1};
+    MeshIndex offset {0, 0, 0};
     bool max_surface {true};
     double distance {INFTY};
     bool operator<(const MeshDistance& o) const
@@ -314,6 +314,8 @@ class StructuredMesh : public Mesh {
     }
   };
 
+  virtual void sanitize_index(MeshIndex& idx) const {};
+
   Position sample_element(int32_t bin, uint64_t* seed) const override
   {
     return sample_element(get_indices_from_bin(bin), seed);
@@ -364,6 +366,13 @@ class StructuredMesh : public Mesh {
   //! \return Array of mesh indices
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const;
 
+  //! Check if mesh indices are inside mesh
+  //
+  //! \param[in] Array of mesh indices
+  //! \param[in] k Suspect axis
+  //! \return are indices inside mesh?
+  virtual bool valid_index(const MeshIndex& ijk, int k) const;
+
   //! Get mesh indices corresponding to a mesh bin
   //
   //! \param[in] bin Mesh bin
@@ -416,9 +425,16 @@ class StructuredMesh : public Mesh {
   virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
     const Position& r0, const Direction& u, double l) const = 0;
 
+  virtual double distance_to_mesh(const MeshIndex& ijk,
+    const std::array<MeshDistance, 4>& distances, double traveled_distance,
+    int& k_max) const;
+
   //! Get a label for the mesh bin
   std::string bin_label(int bin) const override;
 
+  //! Get a label for the surface
+  virtual std::string surface_label(int surface) const;
+
   //! Get mesh dimensions as a tensor
   tensor::Tensor<int> get_shape_tensor() const;
 
@@ -452,6 +468,8 @@ class StructuredMesh : public Mesh {
 
   // Data members
   std::array<int, 3> shape_; //!< Number of mesh elements in each dimension
+  std::vector<std::vector<int>> correlated_axes_ = {{0}, {1}, {2}};
+  std::vector<std::string> axes_labels_ = {"x", "y", "z"};
 
 protected:
 };
@@ -580,6 +598,11 @@ class CylindricalMesh : public PeriodicStructuredMesh {
   CylindricalMesh(hid_t group);
 
   // Overridden methods
+  void sanitize_index(MeshIndex& idx) const override
+  {
+    idx[1] = sanitize_phi(idx[1]);
+  }
+
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
 
   int get_index_in_direction(double r, int i) const override;
@@ -645,6 +668,12 @@ class SphericalMesh : public PeriodicStructuredMesh {
   SphericalMesh(hid_t group);
 
   // Overridden methods
+  void sanitize_index(MeshIndex& idx) const override
+  {
+    idx[1] = sanitize_theta(idx[1]);
+    idx[2] = sanitize_phi(idx[2]);
+  }
+
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
 
   int get_index_in_direction(double r, int i) const override;
@@ -706,6 +735,85 @@ class SphericalMesh : public PeriodicStructuredMesh {
   }
 };
 
+class HexagonalMesh : public PeriodicStructuredMesh {
+public:
+  // Constructors
+  HexagonalMesh() = default;
+  HexagonalMesh(pugi::xml_node node);
+  HexagonalMesh(hid_t group);
+
+  // Overridden methods
+
+  //! Get a label for the mesh bin
+  std::string bin_label(int bin) const override;
+
+  std::string surface_label(int surface) const override;
+
+  int n_bins() const override;
+
+  MeshIndex get_indices(Position r, bool& in_mesh) const override;
+
+  //! Get mesh indices corresponding to a mesh bin
+  //
+  //! \param[in] bin Mesh bin
+  //! \return ijk Mesh indices
+  MeshIndex get_indices_from_bin(int bin) const override;
+
+  bool valid_index(const MeshIndex& ijk, int k) const override;
+
+  int get_bin_from_indices(const MeshIndex& ijk) const override;
+
+  std::string get_mesh_type() const override;
+
+  static const std::string mesh_type;
+
+  Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
+
+  MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
+    const Position& r0, const Direction& u, double l) const override;
+
+  double distance_to_mesh(const MeshIndex& ijk,
+    const std::array<MeshDistance, 4>& distances, double traveled_distance,
+    int& k_max) const override;
+
+  std::pair<vector<double>, vector<double>> plot(
+    Position plot_ll, Position plot_ur) const override;
+
+  void to_hdf5_inner(hid_t group) const override;
+
+  int set_grid();
+
+  // Data members
+  int num_rings_;
+  double pitch_;
+  vector<double> grid_;
+  Direction q_;
+  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};
+
+  StructuredMesh::MeshDistance find_z_crossing(
+    const Position& r, const Direction& u, double l, int shell) const;
+
+  double volume(const MeshIndex& ijk) const override;
+
+  int get_index_in_direction(double r, int i) const override
+  {
+    fatal_error("This function is not implemented");
+    return -1;
+  }
+
+  int get_index_in_z_direction(double z) const;
+};
+
 // Abstract class for unstructured meshes
 class UnstructuredMesh : public Mesh {
 

include/openmc/tallies/filter_meshsurface.h

@@ -22,21 +22,6 @@ class MeshSurfaceFilter : public MeshFilter {
   // Accessors
 
   void set_mesh(int32_t mesh) override;
-
-  enum class MeshDir {
-    OUT_LEFT,   // x min
-    IN_LEFT,    // x min
-    OUT_RIGHT,  // x max
-    IN_RIGHT,   // x max
-    OUT_BACK,   // y min
-    IN_BACK,    // y min
-    OUT_FRONT,  // y max
-    IN_FRONT,   // y max
-    OUT_BOTTOM, // z min
-    IN_BOTTOM,  // z min
-    OUT_TOP,    // z max
-    IN_TOP      // z max
-  };
 };
 
 } // namespace openmc

openmc/filter.py

@@ -30,13 +30,6 @@
     'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
 )
 
-_CURRENT_NAMES = (
-    'x-min out', 'x-min in', 'x-max out', 'x-max in',
-    'y-min out', 'y-min in', 'y-max out', 'y-max in',
-    'z-min out', 'z-min in', 'z-max out', 'z-max in'
-)
-
-
 
 class FilterMeta(ABCMeta):
     """Metaclass for filters that ensures class names are appropriate."""
@@ -79,7 +72,7 @@ def __new__(cls, name, bases, namespace, **kwargs):
 
 
 def _repeat_and_tile(bins, repeat_factor, data_size):
-    filter_bins = np.repeat(bins, repeat_factor)
+    filter_bins = np.repeat(bins, repeat_factor, axis=0)
     tile_factor = data_size // len(filter_bins)
     return np.tile(filter_bins, tile_factor)
 
@@ -992,16 +985,9 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         # Append mesh ID as outermost index of multi-index
         mesh_key = f'mesh {self.mesh.id}'
 
-        # Determine index base (0-based for unstructured, 1-based otherwise)
-        idx_start = 0 if isinstance(self.mesh, openmc.UnstructuredMesh) else 1
-
-        # Generate a multi-index sub-column for each axis
-        for label, dim_size in zip(self.mesh._axis_labels, self.mesh.dimension):
-            filter_dict[mesh_key, label] = _repeat_and_tile(
-                np.arange(idx_start, idx_start + dim_size), stride, data_size)
-            stride *= dim_size
-
-        return pd.DataFrame(filter_dict)
+        columns = [(mesh_key, label) for label in self.mesh._axis_labels]
+        indices = _repeat_and_tile(list(self.mesh.indices), stride, data_size)
+        return pd.DataFrame(indices, columns=columns)
 
     def to_xml_element(self):
         """Return XML Element representing the Filter.
@@ -1282,15 +1268,34 @@ class MeshSurfaceFilter(MeshFilter):
     def shape(self):
         return (self.num_bins,)
 
+    @staticmethod
+    def _current_names(mesh):
+        if isinstance(mesh, openmc.HexagonalMesh):
+            names = []
+            ax0, ax1, ax2, ax3 = mesh._axis_labels
+            pairs = [(ax0, ax1), (ax1, ax2), (ax2, ax0)]
+            for ax0, ax1 in pairs:
+                for minmax0, minmax1 in [('max','min'),('min','max')]:
+                    for inout in ('out','in'):
+                        names.append(f"{ax0}-{minmax0} {ax1}-{minmax1} {inout}")
+            for minmax in ('min', 'max'):
+                for inout in ('out','in'):
+                    names.append(f"{ax3}-{minmax} {inout}")
+            return names
+
+        names = []
+        for ax in mesh._axis_labels:
+            for minmax in ('min', 'max'):
+                for inout in ('out','in'):
+                    names.append(f"{ax}-{minmax} {inout}")
+        return names
+
     @MeshFilter.mesh.setter
     def mesh(self, mesh):
         cv.check_type('filter mesh', mesh, openmc.MeshBase)
         self._mesh = mesh
-
-        # Take the product of mesh indices and current names
-        n_dim = mesh.n_dimension
         self.bins = [mesh_tuple + (surf,) for mesh_tuple, surf in
-                     product(mesh.indices, _CURRENT_NAMES[:4*n_dim])]
+                     product(mesh.indices, self._current_names(mesh))]
 
     def get_pandas_dataframe(self, data_size, stride, **kwargs):
         """Builds a Pandas DataFrame for the Filter's bins.
@@ -1320,47 +1325,19 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
 
         """
-        # Initialize Pandas DataFrame
-        df = pd.DataFrame()
-
-        # Initialize dictionary to build Pandas Multi-index column
-        filter_dict = {}
-
         # Append mesh ID as outermost index of multi-index
         mesh_key = f'mesh {self.mesh.id}'
 
         # Find mesh dimensions - use 3D indices for simplicity
-        n_surfs = 4 * len(self.mesh.dimension)
-        if len(self.mesh.dimension) == 3:
-            nx, ny, nz = self.mesh.dimension
-        elif len(self.mesh.dimension) == 2:
-            nx, ny = self.mesh.dimension
-            nz = 1
-        else:
-            nx = self.mesh.dimension
-            ny = nz = 1
+        n_surfs = len(self._current_names(self._mesh))
 
-        # Generate multi-index sub-column for x-axis
-        filter_dict[mesh_key, 'x'] = _repeat_and_tile(
-            np.arange(1, nx + 1), n_surfs * stride, data_size)
-
-        # Generate multi-index sub-column for y-axis
-        if len(self.mesh.dimension) > 1:
-            filter_dict[mesh_key, 'y'] = _repeat_and_tile(
-                np.arange(1, ny + 1), n_surfs * nx * stride, data_size)
-
-        # Generate multi-index sub-column for z-axis
-        if len(self.mesh.dimension) > 2:
-            filter_dict[mesh_key, 'z'] = _repeat_and_tile(
-                np.arange(1, nz + 1), n_surfs * nx * ny * stride, data_size)
-
-        # Generate multi-index sub-column for surface
-        filter_dict[mesh_key, 'surf'] = _repeat_and_tile(
-            _CURRENT_NAMES[:n_surfs], stride, data_size)
-
-        # Initialize a Pandas DataFrame from the mesh dictionary
-        return pd.concat([df, pd.DataFrame(filter_dict)])
+        columns = [(mesh_key, label) for label in self.mesh._axis_labels]
+        indices = _repeat_and_tile(list(self.mesh.indices), stride*n_surfs, data_size)
+        filter_dict = dict(zip(columns,indices.T))
+        surfs = _repeat_and_tile(self._current_names(self._mesh), stride, data_size)
+        filter_dict[mesh_key, 'surf'] = surfs
 
+        return pd.DataFrame(filter_dict)
 
 class CollisionFilter(Filter):
     """Bins tally events based on the number of collisions.