Adding to mesh abstraction to support other operations external to the mesh source files.

include/openmc/mesh.h

@@ -27,13 +27,11 @@ namespace openmc {
 // Global variables
 //==============================================================================
 
-
 class Mesh;
 
 namespace model {
 
 extern std::vector<std::unique_ptr<Mesh>> meshes;
-
 extern std::unordered_map<int32_t, int32_t> mesh_map;
 
 } // namespace model
@@ -58,12 +56,31 @@ public:
   //! \param[in] group HDF5 group
   virtual void to_hdf5(hid_t group) const = 0;
 
+  //! Determine which surface bins were crossed by a particle
+  //
+  //! \param[in] p Particle to check
+  //! \param[out] bins Surface bins that were crossed
+  virtual void
+  surface_bins_crossed(const Particle* p, std::vector<int>& bins) const = 0;
+
   //! Get bin at a given position in space
   //
   //! \param[in] r Position to get bin for
   //! \return Mesh bin
   virtual int get_bin(Position r) const = 0;
 
+  virtual std::string get_label_for_bin(int bin) const = 0;
+
+  //! Count number of bank sites in each mesh bin / energy bin
+  //
+  //! \param[in] bank Array of bank sites
+  //! \param[out] Whether any bank sites are outside the mesh
+  //! \return Array indicating number of sites in each mesh/energy bin
+  virtual xt::xarray<double>
+  count_sites(const std::vector<Particle::Bank>& bank, bool* outside) const = 0;
+
+  virtual double get_volume_frac(int bin = -1) = 0;
+
   int id_ {-1};  //!< User-specified ID
   int n_dimension_; //!< Number of dimensions
 };
@@ -206,6 +223,11 @@ public:
     bool* outside) const;
 
   // Data members
+
+  //std::string get_label_for_bin(int bin) const;
+
+  //double get_volume_frac(int bin = -1) const;
+
   double volume_frac_; //!< Volume fraction of each mesh element
   xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
   xt::xtensor<double, 1> width_; //!< Width of each mesh element
@@ -282,6 +304,12 @@ public:
 
   bool intersects(Position& r0, Position r1, int* ijk);
 
+  //! Determine which surface bins were crossed by a particle
+  //
+  //! \param[in] p Particle to check
+  //! \param[out] bins Surface bins that were crossed
+  void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
+
   bool point_in_tet(const Position& r, moab::EntityHandle tet) const;
 
   //! Write mesh data to an HDF5 group
@@ -305,6 +333,13 @@ public:
 
   void build_tree(const moab::Range& all_tets);
 
+  std::string get_label_for_bin(int bin) const;
+
+  xt::xarray<double>
+  count_sites(const std::vector<Particle::Bank>& bank, bool* outside) const;
+
+  double get_volume_frac(int bin = -1) const;
+
 private:
   std::string filename_;
   moab::Range ehs_;

src/eigenvalue.cpp

@@ -538,9 +538,12 @@ void ufs_count_sites()
     // distributed so that effectively the production of fission sites is not
     // biased
 
+
     std::size_t n = simulation::ufs_mesh->n_bins();
     double vol_frac = simulation::ufs_mesh->volume_frac_;
     simulation::source_frac = xt::xtensor<double, 1>({n}, vol_frac);
+    // auto s = xt::view(simulation::source_frac, xt::all());
+    // s = m->get_volume_frac();
 
   } else {
     // count number of source sites in each ufs mesh cell
@@ -583,6 +586,7 @@ double ufs_get_weight(const Particle* p)
   if (simulation::source_frac(mesh_bin) != 0.0) {
     return simulation::ufs_mesh->volume_frac_
       / simulation::source_frac(mesh_bin);
+    // return m->get_volume_frac() / simulation::source_frac(mesh_bin);
   } else {
     return 1.0;
   }

src/mesh.cpp

@@ -5,6 +5,7 @@
 #include <cmath>  // for ceil
 #include <memory> // for allocator
 #include <string>
+#include <sstream>
 
 #ifdef OPENMC_MPI
 #include "mpi.h"
@@ -32,7 +33,6 @@ namespace openmc {
 
 namespace model {
 
-
 std::vector<std::unique_ptr<Mesh>> meshes;
 std::unordered_map<int32_t, int32_t> mesh_map;
 
@@ -848,6 +848,23 @@ RegularMesh::count_sites(const Particle::Bank* bank, int64_t length,
   return counts;
 }
 
+std::string RegularMesh::get_label_for_bin(int bin) const {
+  int ijk[n_dimension_];
+  get_indices_from_bin(bin, ijk);
+
+  std::stringstream out;
+  out << "Mesh Index (" << ijk[0];
+  if (n_dimension_ > 1) out << ", " << ijk[1];
+  if (n_dimension_ > 2) out << ", " << ijk[2];
+  out << ")";
+
+  return out.str();
+}
+
+double RegularMesh::get_volume_frac(int bin) const {
+  return volume_frac_;
+}
+
 //==============================================================================
 // RectilinearMesh implementation
 //==============================================================================
@@ -1643,6 +1660,20 @@ UnstructuredMesh::get_tet(Position r) const {
   return 0;
 }
 
+
+//! Determine which surface bins were crossed by a particle
+//
+//! \param[in] p Particle to check
+//! \param[out] bins Surface bins that were crossed
+void UnstructuredMesh::surface_bins_crossed(const Particle* p, std::vector<int>& bins) const {
+  return;
+}
+
+std::string UnstructuredMesh::get_label_for_bin(int bin) {
+  std::string s;
+  return s;
+}
+
 int
 UnstructuredMesh::get_bin(Position r) const {
   moab::EntityHandle tet = get_tet(r);
@@ -1727,6 +1758,17 @@ UnstructuredMesh::get_ent_handle_from_bin(int bin) const {
   return ehs_[bin];
 }
 
+xt::xarray<double>
+UnstructuredMesh::count_sites(const std::vector<Particle::Bank>& bank,
+  bool* outside) const {
+    xt::array<double> out;
+    return out;
+  }
+
+double UnstructuredMesh::get_volume_frac(int bin = -1) const {
+  return 0.0;
+}
+
 #endif
 
 

src/tallies/filter_mesh.cpp

@@ -55,6 +55,7 @@ std::string
 MeshFilter::text_label(int bin) const
 {
   auto& mesh = *model::meshes[mesh_];
+
   int n_dim = mesh.n_dimension_;
 
   std::vector<int> ijk(n_dim);