Overlap detection for plotter

include/openmc/geometry.h

@@ -3,6 +3,8 @@
 
 #include <cmath>
 #include <cstdint>
+#include <unordered_map>
+#include <vector>
 
 #include "openmc/array.h"
 #include "openmc/constants.h"
@@ -13,6 +15,19 @@ namespace openmc {
 class BoundaryInfo;
 class GeometryState;
 
+// OverlapKey to store cell and universe data of a single overlap
+struct OverlapKey {
+  int universe_id;
+  int cell1_id;
+  int cell2_id;
+
+  bool operator==(const OverlapKey& other) const
+  {
+    return universe_id == other.universe_id && cell1_id == other.cell1_id &&
+           cell2_id == other.cell2_id;
+  }
+};
+
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -24,6 +39,10 @@ extern "C" int n_coord_levels; //!< Number of CSG coordinate levels
 
 extern vector<int64_t> overlap_check_count;
 
+// Overlap data structures get cleared every slice_data run
+extern vector<OverlapKey> overlap_keys;
+extern std::unordered_map<OverlapKey, size_t> overlap_key_index;
+
 } // namespace model
 
 //==============================================================================
@@ -38,8 +57,7 @@ inline bool coincident(double d1, double d2)
 //==============================================================================
 //! Check for overlapping cells at a particle's position.
 //==============================================================================
-
-bool check_cell_overlap(GeometryState& p, bool error = true);
+size_t check_cell_overlap(GeometryState& p, bool error = true);
 
 //==============================================================================
 //! Get the cell instance for a particle at the specified universe level
@@ -79,4 +97,16 @@ BoundaryInfo distance_to_boundary(GeometryState& p);
 
 } // namespace openmc
 
+// Hash specialization for use in std::unordered_map
+template<>
+struct std::hash<openmc::OverlapKey> {
+  size_t operator()(const openmc::OverlapKey& k) const noexcept
+  {
+    size_t h = std::hash<int> {}(k.universe_id);
+    h ^= std::hash<int> {}(k.cell1_id) + 0x9e3779b9 + (h << 6) + (h >> 2);
+    h ^= std::hash<int> {}(k.cell2_id) + 0x9e3779b9 + (h << 6) + (h >> 2);
+    return h;
+  }
+};
+
 #endif // OPENMC_GEOMETRY_H

include/openmc/plot.h

@@ -5,6 +5,7 @@
 #include <sstream>
 #include <unordered_map>
 #include <unordered_set>
+#include <vector>
 
 #include "openmc/tensor.h"
 #include "pugixml.hpp"
@@ -155,7 +156,7 @@ struct IdData {
   // 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);
+  void set_overlap(size_t y, size_t x, size_t overlap_idx);
 
   // Members
   tensor::Tensor<int32_t> data_; //!< 2D array of cell & material ids
@@ -168,7 +169,7 @@ struct PropertyData {
   // 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);
+  void set_overlap(size_t y, size_t x, size_t overlap_idx);
 
   // Members
   tensor::Tensor<double> data_; //!< 2D array of temperature & density data
@@ -181,7 +182,7 @@ struct RasterData {
   // 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);
+  void set_overlap(size_t y, size_t x, size_t overlap_idx);
 
   // Members
   tensor::Tensor<int32_t>
@@ -278,8 +279,11 @@ T SlicePlotBase::get_map(int32_t filter_index) const
         if (found_cell) {
           data.set_value(y, x, p, j, filter, &match);
         }
-        if (show_overlaps_ && check_cell_overlap(p, false)) {
-          data.set_overlap(y, x);
+        if (show_overlaps_) {
+          auto overlap_idx = check_cell_overlap(p, false);
+          if (overlap_idx != SIZE_MAX) {
+            data.set_overlap(y, x, overlap_idx);
+          }
         }
       } // inner for
     }

openmc/lib/plot.py

@@ -1,6 +1,7 @@
 from collections.abc import Mapping
 from ctypes import (c_bool, c_int, c_size_t, c_int32,
-                    c_double, c_uint8, Structure, POINTER)
+                    c_double, c_uint8, Structure, POINTER,
+                    byref)
 from weakref import WeakValueDictionary
 
 from ..exceptions import AllocationError, InvalidIDError
@@ -255,6 +256,33 @@ def property_map(plot):
     return prop_data
 
 
+# Python wrappings for overlap functions
+
+def slice_data_overlap_count():
+    count = c_size_t()
+    _dll.openmc_slice_data_overlap_count(byref(count))
+    return count.value
+
+def slice_data_overlap_info():
+    n = slice_data_overlap_count()
+    overlap_info = np.empty(n * 3, dtype=np.int32)
+
+    if n > 0:
+        _dll.openmc_slice_data_overlap_info(
+            n,
+            overlap_info.ctypes.data_as(POINTER(c_int32)),
+        )
+
+    return overlap_info, n
+
+_dll.openmc_slice_data_overlap_count.argtypes = [POINTER(c_size_t)]
+_dll.openmc_slice_data_overlap_count.restype = c_int
+_dll.openmc_slice_data_overlap_count.errcheck = _error_handler
+
+_dll.openmc_slice_data_overlap_info.argtypes = [c_size_t, POINTER(c_int32)]
+_dll.openmc_slice_data_overlap_info.restype = c_int
+_dll.openmc_slice_data_overlap_info.errcheck = _error_handler
+
 _dll.openmc_get_plot_index.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_get_plot_index.restype = c_int
 _dll.openmc_get_plot_index.errcheck = _error_handler

src/geometry.cpp

@@ -26,16 +26,22 @@ int n_coord_levels;
 
 vector<int64_t> overlap_check_count;
 
+vector<OverlapKey> overlap_keys;
+std::unordered_map<OverlapKey, size_t> overlap_key_index;
+
 } // namespace model
 
 //==============================================================================
 // Non-member functions
 //==============================================================================
 
-bool check_cell_overlap(GeometryState& p, bool error)
+size_t check_cell_overlap(GeometryState& p, bool error)
 {
   int n_coord = p.n_coord();
 
+  // If no overlap found, return a nonphysical index
+  size_t overlap_index = SIZE_MAX;
+
   // Loop through each coordinate level
   for (int j = 0; j < n_coord; j++) {
     Universe& univ = *model::universes[p.coord(j).universe()];
@@ -44,21 +50,39 @@ bool check_cell_overlap(GeometryState& p, bool error)
     for (auto index_cell : univ.cells_) {
       Cell& c = *model::cells[index_cell];
       if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {
+#pragma omp atomic
+        ++model::overlap_check_count[index_cell];
+
         if (index_cell != p.coord(j).cell()) {
           if (error) {
             fatal_error(
               fmt::format("Overlapping cells detected: {}, {} on universe {}",
                 c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));
           }
-          return true;
+
+          // With no fatal error (plotter is calling), now adds overlaps and
+          // calls them; ensures order does not matter when making overlap key
+          int cell_a = model::cells[index_cell]->id_;
+          int cell_b = model::cells[p.coord(j).cell()]->id_;
+          int a = std::min(cell_a, cell_b);
+          int b = std::max(cell_a, cell_b);
+          OverlapKey key {univ.id_, a, b};
+
+          auto it = model::overlap_key_index.find(key);
+          if (it != model::overlap_key_index.end()) {
+            overlap_index = it->second; // already exists, reuse index
+          } else {
+            size_t idx = model::overlap_keys.size();
+            model::overlap_keys.push_back(key);
+            model::overlap_key_index[key] = idx;
+            overlap_index = idx;
+          }
+          break;
         }
-#pragma omp atomic
-        ++model::overlap_check_count[index_cell];
       }
     }
   }
-
-  return false;
+  return overlap_index;
 }
 
 //==============================================================================

src/plot.cpp

@@ -73,7 +73,7 @@ void IdData::set_value(size_t y, size_t x, const Particle& p, int level,
   }
 }
 
-void IdData::set_overlap(size_t y, size_t x)
+void IdData::set_overlap(size_t y, size_t x, size_t /*overlap_idx*/)
 {
   for (size_t k = 0; k < data_.shape(2); ++k)
     data_(y, x, k) = OVERLAP;
@@ -94,7 +94,7 @@ void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level,
   }
 }
 
-void PropertyData::set_overlap(size_t y, size_t x)
+void PropertyData::set_overlap(size_t y, size_t x, size_t /*overlap_idx*/)
 {
   data_(y, x) = OVERLAP;
 }
@@ -154,12 +154,12 @@ void RasterData::set_value(size_t y, size_t x, const Particle& p, int level,
   }
 }
 
-void RasterData::set_overlap(size_t y, size_t x)
+void RasterData::set_overlap(size_t y, size_t x, size_t overlap_idx)
 {
   // Set cell, instance, and material to OVERLAP, but preserve filter bin
   id_data_(y, x, 0) = OVERLAP;
   id_data_(y, x, 1) = OVERLAP;
-  id_data_(y, x, 2) = OVERLAP;
+  id_data_(y, x, 2) = OVERLAP - overlap_idx - 1;
   // Note: id_data_(y, x, 3) is NOT overwritten - preserves filter bin for tally
   // plotting
 
@@ -177,6 +177,8 @@ std::unordered_map<int, int> plot_map;
 vector<std::unique_ptr<PlottableInterface>> plots;
 uint64_t plotter_seed = 1;
 
+std::unique_ptr<RasterData> last_slice_data;
+
 } // namespace model
 
 //==============================================================================
@@ -1980,6 +1982,10 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3],
     model::overlap_check_count.resize(model::cells.size());
   }
 
+  if (color_overlaps) {
+    settings::check_overlaps = true;
+  }
+
   try {
     // Create a temporary SlicePlotBase object to reuse get_map logic
     SlicePlotBase plot_params;
@@ -1991,17 +1997,20 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3],
     plot_params.show_overlaps_ = color_overlaps;
     plot_params.slice_level_ = level;
 
+    // Clear overlap data structures on new slice call
+    model::overlap_keys.clear();
+    model::overlap_key_index.clear();
+
     // Use get_map<RasterData> to generate data
     auto data = plot_params.get_map<RasterData>(filter_index);
-
-    // Copy geometry data
     std::copy(data.id_data_.begin(), data.id_data_.end(), geom_data);
 
     // Copy property data if requested
     if (property_data != nullptr) {
       std::copy(
         data.property_data_.begin(), data.property_data_.end(), property_data);
     }
+
   } catch (const std::exception& e) {
     set_errmsg(e.what());
     return OPENMC_E_UNASSIGNED;
@@ -2010,6 +2019,32 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3],
   return 0;
 }
 
+// Gets the number of overlaps that we need data for
+extern "C" int openmc_slice_data_overlap_count(size_t* count)
+{
+  if (!count) {
+    set_errmsg("Null pointer passed for overlap count.");
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+  *count = model::overlap_keys.size();
+
+  return 0;
+}
+
+// Plotter pre-allocates array size based on what is returned with
+// overlap_count; populates an array of size 3*count
+extern "C" int openmc_slice_data_overlap_info(
+  size_t count, int32_t* overlap_info)
+{
+  for (size_t i = 0; i < count; ++i) {
+    overlap_info[i * 3] = model::overlap_keys[i].universe_id;
+    overlap_info[i * 3 + 1] = model::overlap_keys[i].cell1_id;
+    overlap_info[i * 3 + 2] = model::overlap_keys[i].cell2_id;
+  }
+
+  return 0;
+}
+
 extern "C" int openmc_get_plot_index(int32_t id, int32_t* index)
 {
   auto it = model::plot_map.find(id);