Updates to get reg mesh working as an abstracted class

include/openmc/mesh.h

@@ -81,6 +81,8 @@ public:
 
   virtual double get_volume_frac(int bin = -1) const = 0;
 
+  virtual int num_bins() const = 0;
+
   int id_ {-1};  //!< User-specified ID
   int n_dimension_; //!< Number of dimensions
 };
@@ -228,6 +230,8 @@ public:
 
   //double get_volume_frac(int bin = -1) const;
 
+  int num_bins() 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
@@ -366,6 +370,8 @@ void read_meshes(pugi::xml_node root);
 //! \param[in] group HDF5 group
 void meshes_to_hdf5(hid_t group);
 
+RegularMesh* get_regular_mesh(int32_t index);
+
 void free_memory_mesh();
 
 } // namespace openmc

src/mesh.cpp

@@ -865,6 +865,12 @@ double RegularMesh::get_volume_frac(int bin) const {
   return volume_frac_;
 }
 
+int RegularMesh::num_bins() const {
+  int n_bins = 1;
+  for (auto v : shape_) n_bins *= v;
+  return n_bins;
+}
+
 //==============================================================================
 // RectilinearMesh implementation
 //==============================================================================
@@ -1402,6 +1408,11 @@ check_regular_mesh(int32_t index, RegularMesh** mesh)
 // C API functions
 //==============================================================================
 
+RegularMesh* get_regular_mesh(int32_t index) {
+  return dynamic_cast<RegularMesh*>(model::meshes[index].get());
+}
+
+
 //! Extend the meshes array by n elements
 extern "C" int
 openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
@@ -1451,6 +1462,11 @@ openmc_mesh_set_id(int32_t index, int32_t id)
 extern "C" int
 openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
 {
+  if (index < 0 || index >= model::meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
   RegularMesh* mesh;
   if (int err = check_regular_mesh(index, &mesh)) return err;
   *dims = mesh->shape_.data();
@@ -1469,7 +1485,6 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
   std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
   mesh->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
   mesh->n_dimension_ = mesh->shape_.size();
-
   return 0;
 }
 

src/plot.cpp

@@ -696,6 +696,7 @@ void draw_mesh_lines(Plot pl, ImageData& data)
 
   Position width = ur_plot - ll_plot;
 
+
   // Find the (axis-aligned) lines of the mesh that intersect this plot.
   auto axis_lines = model::meshes[pl.index_meshlines_mesh_]
     ->plot(ll_plot, ur_plot);
@@ -713,6 +714,49 @@ void draw_mesh_lines(Plot pl, ImageData& data)
     ax2_min = 0;
     ax2_max = pl.pixels_[1];
   }
+  // auto m = get_regular_mesh(pl.index_meshlines_mesh_);
+
+  // int ijk_ll[3], ijk_ur[3];
+  // bool in_mesh;
+  // m->get_indices(ll_plot, &(ijk_ll[0]), &in_mesh);
+  // m->get_indices(ur_plot, &(ijk_ur[0]), &in_mesh);
+
+  // // Fortran/C++ index correction
+  // ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++;
+
+  // Position r_ll, r_ur;
+  // // sweep through all meshbins on this plane and draw borders
+  // for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) {
+  //   for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) {
+  //     // check if we're in the mesh for this ijk
+  //     if (i > 0 && i <= m->shape_[outer] && j >0 && j <= m->shape_[inner] ) {
+  //       int outrange[3], inrange[3];
+  //       // get xyz's of lower left and upper right of this mesh cell
+  //       r_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1);
+  //       r_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1);
+  //       r_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i;
+  //       r_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j;
+
+  //       // map the xyz ranges to pixel ranges
+  //       double frac = (r_ll[outer] - ll_plot[outer]) / width[outer];
+  //       outrange[0] = int(frac * double(pl.pixels_[0]));
+  //       frac = (r_ur[outer] - ll_plot[outer]) / width[outer];
+  //       outrange[1] = int(frac * double(pl.pixels_[0]));
+
+  //       frac = (r_ur[inner] - ll_plot[inner]) / width[inner];
+  //       inrange[0] = int((1. - frac) * (double)pl.pixels_[1]);
+  //       frac = (r_ll[inner] - ll_plot[inner]) / width[inner];
+  //       inrange[1] = int((1. - frac) * (double)pl.pixels_[1]);
+
+  //       // draw lines
+  //       for (int out_ = outrange[0]; out_ <= outrange[1]; out_++) {
+  //         for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
+  //           data(out_, inrange[0] + plus) = rgb;
+  //           data(out_, inrange[1] + plus) = rgb;
+  //           data(out_, inrange[0] - plus) = rgb;
+  //           data(out_, inrange[1] - plus) = rgb;
+  //         }
+  //       }
 
   // Iterate across the first axis and draw lines.
   for (auto ax1_val : axis_lines.first) {