Refactor and fix RegularMesh meshline plots

include/openmc/mesh.h

@@ -85,6 +85,9 @@ public:
   //! Get the number of mesh cell surfaces.
   virtual int n_surface_bins() const = 0;
 
+  virtual std::pair<std::vector<double>, std::vector<double>>
+  plot(Position plot_ll, Position plot_ur) const = 0;
+
   //! Write mesh data to an HDF5 group
   //
   //! \param[in] group HDF5 group
@@ -127,6 +130,9 @@ public:
 
   int n_surface_bins() const override;
 
+  std::pair<std::vector<double>, std::vector<double>>
+  plot(Position plot_ll, Position plot_ur) const override;
+
   void to_hdf5(hid_t group) const override;
 
   // New methods
@@ -187,6 +193,9 @@ public:
 
   int n_surface_bins() const override;
 
+  std::pair<std::vector<double>, std::vector<double>>
+  plot(Position plot_ll, Position plot_ur) const override;
+
   void to_hdf5(hid_t group) const override;
 
   // New methods

src/mesh.cpp

@@ -729,6 +729,40 @@ void RegularMesh::surface_bins_crossed(const Particle* p,
   }
 }
 
+std::pair<std::vector<double>, std::vector<double>>
+RegularMesh::plot(Position plot_ll, Position plot_ur) const
+{
+  std::array<int, 2> axes {-1, -1};
+  if (plot_ur.z == plot_ll.z) {
+    axes[0] = 0;
+    if (n_dimension_ > 1) axes[1] = 1;
+  } else if (plot_ur.y == plot_ll.y) {
+    axes[0] = 0;
+    if (n_dimension_ > 2) axes[1] = 2;
+  } else if (plot_ur.x == plot_ll.x) {
+    if (n_dimension_ > 1) axes[0] = 1;
+    if (n_dimension_ > 2) axes[1] = 2;
+  } else {
+    fatal_error("Can only plot mesh lines on an axis-aligned plot");
+  }
+
+  std::array<std::vector<double>, 2> axis_lines;
+  for (int i_ax = 0; i_ax < 2; ++i_ax) {
+    int axis = axes[i_ax];
+    if (axis == -1) continue;
+    std::vector<double>& lines {axis_lines[i_ax]};
+
+    double coord = lower_left_[axis];
+    for (int i = 0; i < shape_[axis] + 1; ++i) {
+      if (coord >= plot_ll[axis] && coord <= plot_ur[axis])
+        lines.push_back(coord);
+      coord += width_[axis];
+    }
+  }
+
+  return {axis_lines[0], axis_lines[1]};
+}
+
 void RegularMesh::to_hdf5(hid_t group) const
 {
   hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
@@ -1125,6 +1159,11 @@ int RectilinearMesh::n_surface_bins() const
   return 4 * n_dimension_ * n_bins();
 }
 
+std::pair<std::vector<double>, std::vector<double>>
+RectilinearMesh::plot(Position plot_ll, Position plot_ur) const
+{
+}
+
 void RectilinearMesh::to_hdf5(hid_t group) const
 {
   hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));

src/plot.cpp

@@ -698,19 +698,19 @@ void draw_mesh_lines(Plot pl, ImageData& data)
   RGBColor rgb;
   rgb = pl.meshlines_color_;
 
-  int outer, inner;
+  int ax1, ax2;
   switch(pl.basis_) {
   case PlotBasis::xy :
-    outer = 0;
-    inner = 1;
+    ax1 = 0;
+    ax2 = 1;
     break;
   case PlotBasis::xz :
-    outer = 0;
-    inner = 2;
+    ax1 = 0;
+    ax2 = 2;
     break;
   case PlotBasis::yz :
-    outer = 1;
-    inner = 2;
+    ax1 = 1;
+    ax2 = 2;
     break;
   default:
     UNREACHABLE();
@@ -719,70 +719,79 @@ void draw_mesh_lines(Plot pl, ImageData& data)
   Position ll_plot {pl.origin_};
   Position ur_plot {pl.origin_};
 
-  ll_plot[outer] -= pl.width_[0] / 2.;
-  ll_plot[inner] -= pl.width_[1] / 2.;
-  ur_plot[outer] += pl.width_[0] / 2.;
-  ur_plot[inner] += pl.width_[1] / 2.;
+  ll_plot[ax1] -= pl.width_[0] / 2.;
+  ll_plot[ax2] -= pl.width_[1] / 2.;
+  ur_plot[ax1] += pl.width_[0] / 2.;
+  ur_plot[ax2] += pl.width_[1] / 2.;
 
   Position width = ur_plot - ll_plot;
 
-  auto& m = *dynamic_cast<RegularMesh*>(model::meshes[pl.index_meshlines_mesh_].get());
+  // Find the (axis-aligned) lines of the mesh that intersect this plot.
+  std::pair<std::vector<double>, std::vector<double>> axis_lines;
+  axis_lines = model::meshes[pl.index_meshlines_mesh_]->plot(ll_plot, ur_plot);
 
-  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;
-          }
-        }
-
-        for (int in_ = inrange[0]; in_ <= inrange[1]; in_++) {
-          for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
-            data(outrange[0] + plus, in_) = rgb;
-            data(outrange[1] + plus, in_) = rgb;
-            data(outrange[0] - plus, in_) = rgb;
-            data(outrange[1] - plus, in_) = rgb;
-          }
-        }
-
-      } // end if(in mesh)
+  // Draw lines perpendicular to the first axis.
+  {
+    // Find the bounds along the second axis (accounting for low-D meshes).
+    int ax2_min, ax2_max;
+    if (axis_lines.second.size() > 0) {
+      double frac = (axis_lines.second.front() - ll_plot[ax2]) / width[ax2];
+      ax2_min = frac * pl.pixels_[1];
+      if (ax2_min < 0) ax2_min = 0;
+      frac = (axis_lines.second.back() - ll_plot[ax2]) / width[ax2];
+      ax2_max = frac * pl.pixels_[1];
+      if (ax2_max > pl.pixels_[1]) ax2_max = pl.pixels_[1];
+    } else {
+      ax2_min = 0;
+      ax2_max = pl.pixels_[1];
     }
-  } // end outer loops
-} // end draw_mesh_lines
+
+    // Iterate across the first axis and draw lines.
+    for (auto ax1_val : axis_lines.first) {
+      double frac = (ax1_val - ll_plot[ax1]) / width[ax1];
+      int ax1_ind = frac * pl.pixels_[0];
+      for (int ax2_ind = ax2_min; ax2_ind < ax2_max; ++ax2_ind) {
+        for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
+          if (ax1_ind+plus >= 0 && ax1_ind+plus < pl.pixels_[0])
+            data(ax1_ind+plus, ax2_ind) = rgb;
+          if (ax1_ind-plus >= 0 && ax1_ind-plus < pl.pixels_[0])
+            data(ax1_ind-plus, ax2_ind) = rgb;
+        }
+      }
+    }
+  }
+
+  // Draw lines perpendicular to the second axis.
+  {
+    // Find the bounds along the first axis.
+    int ax1_min, ax1_max;
+    if (axis_lines.first.size() > 0) {
+      double frac = (axis_lines.first.front() - ll_plot[ax1]) / width[ax1];
+      ax1_min = frac * pl.pixels_[0];
+      if (ax1_min < 0) ax1_min = 0;
+      frac = (axis_lines.first.back() - ll_plot[ax1]) / width[ax1];
+      ax1_max = frac * pl.pixels_[0];
+      if (ax1_max > pl.pixels_[0]) ax1_max = pl.pixels_[0];
+    } else {
+      ax1_min = 0;
+      ax1_max = pl.pixels_[0];
+    }
+
+    // Iterate across the second axis and draw lines.
+    for (auto ax2_val : axis_lines.second) {
+      double frac = (ax2_val - ll_plot[ax2]) / width[ax2];
+      int ax2_ind = frac * pl.pixels_[1];
+      for (int ax1_ind = ax1_min; ax1_ind < ax1_max; ++ax1_ind) {
+        for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
+          if (ax2_ind+plus >= 0 && ax2_ind+plus < pl.pixels_[1])
+            data(ax1_ind, ax2_ind+plus) = rgb;
+          if (ax2_ind-plus >= 0 && ax2_ind-plus < pl.pixels_[1])
+            data(ax1_ind, ax2_ind-plus) = rgb;
+        }
+      }
+    }
+  }
+}
 
 //==============================================================================
 // CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D

tests/regression_tests/plot/plots.xml

@@ -6,7 +6,10 @@
     <width>25 25</width>
     <pixels>200 200</pixels>
     <color id="1" rgb="255 0 0" /> <!-- Red -->
+    <!--
     <meshlines meshtype="entropy" linewidth="0" />
+    -->
+    <meshlines meshtype="tally" id="2" linewidth="2" />
   </plot>
 
   <plot id="2" basis="xz">

tests/regression_tests/plot/results_true.dat

@@ -1 +1 @@
-6b3eb36488b995b42233e6b7c0164cf6c92a1823b3c91985b97fd076f86c0aad93fbdb74e70026f5f12c61a6aee52a09588171a7fd839511ee7d9efc3952beb2
+ae51cbf147d69804b03f9d88c03d98533665a9822438d1d8b70b40928a1db35d20e18d12ee697b5e1e217799d6e8c352524e8d2227328e719e85c4e6d96f825c