Introduce new C API function for slice plots

include/openmc/capi.h

@@ -125,6 +125,9 @@ int openmc_nuclide_name(int index, const char** name);
 int openmc_plot_geometry();
 int openmc_id_map(const void* slice, int32_t* data_out);
 int openmc_property_map(const void* slice, double* data_out);
+int openmc_slice_plot(const double origin[3], const double u_span[3],
+  const double v_span[3], const size_t pixels[2], bool color_overlaps,
+  int level, int32_t filter_index, int32_t* geom_data, double* property_data);
 int openmc_get_plot_index(int32_t id, int32_t* index);
 int openmc_plot_get_id(int32_t index, int32_t* id);
 int openmc_plot_set_id(int32_t index, int32_t id);

include/openmc/plot.h

@@ -17,6 +17,8 @@
 #include "openmc/particle.h"
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/filter_match.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -147,36 +149,56 @@ class PlottableInterface {
 
 struct IdData {
   // Constructor
-  IdData(size_t h_res, size_t v_res);
+  IdData(size_t h_res, size_t v_res, bool include_filter = false);
 
   // Methods
-  void set_value(size_t y, size_t x, const GeometryState& p, int level);
+  void set_value(size_t y, size_t x, const Particle& p, int level,
+    Filter* filter = nullptr, FilterMatch* match = nullptr);
   void set_overlap(size_t y, size_t x);
 
   // Members
-  xt::xtensor<int32_t, 3> data_; //!< 2D array of cell & material ids
+  xt::xtensor<int32_t, 3> data_; //!< 3D array of cell ID, cell instance,
+                                 //!< and material ID
 };
 
 struct PropertyData {
   // Constructor
-  PropertyData(size_t h_res, size_t v_res);
+  PropertyData(size_t h_res, size_t v_res, bool include_filter = false);
 
   // Methods
-  void set_value(size_t y, size_t x, const GeometryState& p, int level);
+  void set_value(size_t y, size_t x, const Particle& p, int level,
+    Filter* filter = nullptr, FilterMatch* match = nullptr);
   void set_overlap(size_t y, size_t x);
 
   // Members
   xt::xtensor<double, 3> data_; //!< 2D array of temperature & density data
 };
 
+struct RasterData {
+  // Constructor
+  RasterData(size_t h_res, size_t v_res, bool include_filter = false);
+
+  // Methods
+  void set_value(size_t y, size_t x, const Particle& p, int level,
+    Filter* filter = nullptr, FilterMatch* match = nullptr);
+  void set_overlap(size_t y, size_t x);
+
+  // Members
+  xt::xtensor<int32_t, 3>
+    id_data_; //!< [v_res, h_res, 3 or 4]: cell, instance, mat, [filter_bin]
+  xt::xtensor<double, 3>
+    property_data_;     //!< [v_res, h_res, 2]: temperature, density
+  bool include_filter_; //!< Whether filter bin index is included
+};
+
 //===============================================================================
 // Plot class
 //===============================================================================
 
 class SlicePlotBase {
 public:
   template<class T>
-  T get_map() const;
+  T get_map(int32_t filter_index = -1) const;
 
   enum class PlotBasis { xy = 1, xz = 2, yz = 3 };
 
@@ -188,69 +210,66 @@ class SlicePlotBase {
   // Members
 public:
   Position origin_;           //!< Plot origin in geometry
-  Position width_;            //!< Plot width in geometry
-  PlotBasis basis_;           //!< Plot basis (XY/XZ/YZ)
+  Position u_span_;           //!< Full-width span vector in geometry
+  Position v_span_;           //!< Full-height span vector in geometry
+  Position width_;            //!< Axis-aligned plot width in geometry
+  PlotBasis basis_;           //!< Plot basis (XY/XZ/YZ) for axis-aligned slices
   array<size_t, 3> pixels_;   //!< Plot size in pixels
   bool slice_color_overlaps_; //!< Show overlapping cells?
   int slice_level_ {-1};      //!< Plot universe level
 private:
 };
 
 template<class T>
-T SlicePlotBase::get_map() const
+T SlicePlotBase::get_map(int32_t filter_index) const
 {
 
   size_t width = pixels_[0];
   size_t height = pixels_[1];
 
-  // get pixel size
-  double in_pixel = (width_[0]) / static_cast<double>(width);
-  double out_pixel = (width_[1]) / static_cast<double>(height);
+  // Determine if filter is being used
+  bool include_filter = (filter_index >= 0);
+  Filter* filter = nullptr;
+  if (include_filter) {
+    filter = model::tally_filters[filter_index].get();
+  }
 
   // size data array
-  T data(width, height);
-
-  // setup basis indices and initial position centered on pixel
-  int in_i, out_i;
-  Position xyz = origin_;
-  switch (basis_) {
-  case PlotBasis::xy:
-    in_i = 0;
-    out_i = 1;
-    break;
-  case PlotBasis::xz:
-    in_i = 0;
-    out_i = 2;
-    break;
-  case PlotBasis::yz:
-    in_i = 1;
-    out_i = 2;
-    break;
-  default:
-    UNREACHABLE();
-  }
+  T data(width, height, include_filter);
 
-  // set initial position
-  xyz[in_i] = origin_[in_i] - width_[0] / 2. + in_pixel / 2.;
-  xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.;
+  // compute pixel steps and top-left pixel center
+  Position u_step = u_span_ / static_cast<double>(width);
+  Position v_step = v_span_ / static_cast<double>(height);
+
+  Position start =
+    origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step;
+
+  // Validate that span vectors define a valid plane
+  Position cross = u_span_.cross(v_span_);
+  if (cross.norm() == 0.0) {
+    fatal_error("Slice span vectors are invalid (zero area).");
+  }
 
-  // arbitrary direction
-  Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0};
+  // Use an arbitrary direction that is not aligned with any coordinate axis.
+  // The direction has no physical meaning for plotting but is used by
+  // Surface::sense() to break ties when a pixel is coincident with a surface.
+  Direction dir = {1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 0.0};
 
 #pragma omp parallel
   {
-    GeometryState p;
-    p.r() = xyz;
+    Particle p;
+    p.r() = start;
     p.u() = dir;
     p.coord(0).universe() = model::root_universe;
     int level = slice_level_;
     int j {};
+    FilterMatch match;
 
 #pragma omp for
     for (int y = 0; y < height; y++) {
-      p.r()[out_i] = xyz[out_i] - out_pixel * y;
+      Position row = start - v_step * static_cast<double>(y);
       for (int x = 0; x < width; x++) {
-        p.r()[in_i] = xyz[in_i] + in_pixel * x;
+        p.r() = row + u_step * static_cast<double>(x);
         p.n_coord() = 1;
         // local variables
         bool found_cell = exhaustive_find_cell(p);
@@ -259,7 +278,7 @@ T SlicePlotBase::get_map() const
           j = level;
         }
         if (found_cell) {
-          data.set_value(y, x, p, j);
+          data.set_value(y, x, p, j, filter, &match);
         }
         if (slice_color_overlaps_ && check_cell_overlap(p, false)) {
           data.set_overlap(y, x);