Merge pull request #1260 from pshriwise/overlap_plotting

Overlap plotting

include/openmc/geometry.h

@@ -47,7 +47,7 @@ inline bool coincident(double d1, double d2) {
 //! Check for overlapping cells at a particle's position.
 //==============================================================================
 
-bool check_cell_overlap(Particle* p);
+bool check_cell_overlap(Particle* p, bool error=true);
 
 //==============================================================================
 //! Locate a particle in the geometry tree and set its geometry data fields.

include/openmc/plot.h

@@ -49,10 +49,19 @@ struct RGBColor {
     blue = v[2];
   }
 
+  bool operator ==(const RGBColor& other) {
+    return red == other.red && green == other.green && blue == other.blue;
+  }
+
   // Members
   uint8_t red, green, blue;
 };
 
+// some default colors
+const RGBColor WHITE {255, 255, 255};
+const RGBColor RED {255,   0,   0};
+
+
 typedef xt::xtensor<RGBColor, 2> ImageData;
 
 struct IdData {
@@ -61,6 +70,7 @@ struct IdData {
 
   // Methods
   void set_value(size_t y, size_t x, const Particle& p, int level);
+  void set_overlap(size_t y, size_t x);
 
   // Members
   xt::xtensor<int32_t, 3> data_; //!< 2D array of cell & material ids
@@ -72,6 +82,7 @@ struct PropertyData {
 
   // Methods
   void set_value(size_t y, size_t x, const Particle& p, int level);
+  void set_overlap(size_t y, size_t x);
 
   // Members
   xt::xtensor<double, 3> data_; //!< 2D array of temperature & density data
@@ -106,6 +117,7 @@ public:
   Position width_; //!< Plot width in geometry
   PlotBasis basis_; //!< Plot basis (XY/XZ/YZ)
   std::array<size_t, 3> pixels_; //!< Plot size in pixels
+  bool color_overlaps_; //!< Show overlapping cells?
   int level_; //!< Plot universe level
 };
 
@@ -173,6 +185,9 @@ T PlotBase::get_map() const {
         if (found_cell) {
           data.set_value(y, x, p, j);
         }
+        if (color_overlaps_ && check_cell_overlap(&p, false)) {
+          data.set_overlap(y, x);
+        }
       } // inner for
     } // outer for
   } // omp parallel
@@ -200,6 +215,7 @@ private:
   void set_user_colors(pugi::xml_node plot_node);
   void set_meshlines(pugi::xml_node plot_node);
   void set_mask(pugi::xml_node plot_node);
+  void set_overlap_color(pugi::xml_node plot_node);
 
 // Members
 public:
@@ -207,9 +223,10 @@ public:
   PlotType type_; //!< Plot type (Slice/Voxel)
   PlotColorBy color_by_; //!< Plot coloring (cell/material)
   int meshlines_width_; //!< Width of lines added to the plot
-  int index_meshlines_mesh_; //!< Index of the mesh to draw on the plot
+  int index_meshlines_mesh_ {-1}; //!< Index of the mesh to draw on the plot
   RGBColor meshlines_color_; //!< Color of meshlines on the plot
-  RGBColor not_found_; //!< Plot background color
+  RGBColor not_found_ {WHITE}; //!< Plot background color
+  RGBColor overlap_color_ {RED}; //!< Plot overlap color
   std::vector<RGBColor> colors_; //!< Plot colors
   std::string path_plot_; //!< Plot output filename
 };

openmc/capi/plot.py

@@ -1,4 +1,5 @@
-from ctypes import c_int, c_size_t, c_int32, c_double, Structure, POINTER
+from ctypes import (c_bool, c_int, c_size_t, c_int32,
+                    c_double, Structure, POINTER)
 
 from . import _dll
 from .error import _error_handler
@@ -84,10 +85,12 @@ class _PlotBase(Structure):
                 ('width_', _Position),
                 ('basis_', c_int),
                 ('pixels_', 3*c_size_t),
+                ('color_overlaps_', c_bool),
                 ('level_', c_int)]
 
     def __init__(self):
         self.level_ = -1
+        self.color_overlaps_ = False
 
     @property
     def origin(self):
@@ -112,32 +115,6 @@ class _PlotBase(Structure):
 
         raise ValueError("Plot basis {} is invalid".format(self.basis_))
 
-    @property
-    def h_res(self):
-        return self.pixels_[0]
-
-    @property
-    def v_res(self):
-        return self.pixels_[1]
-
-    @property
-    def level(self):
-        return int(self.level_)
-
-    @origin.setter
-    def origin(self, origin):
-        self.origin_.x = origin[0]
-        self.origin_.y = origin[1]
-        self.origin_.z = origin[2]
-
-    @width.setter
-    def width(self, width):
-        self.width_.x = width
-
-    @height.setter
-    def height(self, height):
-        self.width_.y = height
-
     @basis.setter
     def basis(self, basis):
         if isinstance(basis, str):
@@ -164,6 +141,40 @@ class _PlotBase(Structure):
         raise ValueError("{} of type {} is an"
                          " invalid plot basis".format(basis, type(basis)))
 
+    @property
+    def h_res(self):
+        return self.pixels_[0]
+
+    @property
+    def v_res(self):
+        return self.pixels_[1]
+
+    @property
+    def level(self):
+        return int(self.level_)
+
+    @property
+    def color_overlaps(self):
+        return self.color_overlaps_
+
+    @color_overlaps.setter
+    def color_overlaps(self, color_overlaps):
+        self.color_overlaps_ = color_overlaps
+
+    @origin.setter
+    def origin(self, origin):
+        self.origin_.x = origin[0]
+        self.origin_.y = origin[1]
+        self.origin_.z = origin[2]
+
+    @width.setter
+    def width(self, width):
+        self.width_.x = width
+
+    @height.setter
+    def height(self, height):
+        self.width_.y = height
+
     @h_res.setter
     def h_res(self, h_res):
         self.pixels_[0] = h_res
@@ -176,6 +187,14 @@ class _PlotBase(Structure):
     def level(self, level):
         self.level_ = level
 
+    @property
+    def color_overlaps(self):
+        return self.color_overlaps_
+
+    @color_overlaps.setter
+    def color_overlaps(self, val):
+        self.color_overlaps_ = val
+
     def __repr__(self):
         out_str = ["-----",
                    "Plot:",
@@ -186,6 +205,7 @@ class _PlotBase(Structure):
                    "Basis: {}".format(self.basis),
                    "HRes: {}".format(self.h_res),
                    "VRes: {}".format(self.v_res),
+                   "Color Overlaps: {}".format(self.color_overlaps),
                    "Level: {}".format(self.level)]
         return '\n'.join(out_str)
 

openmc/plots.py

@@ -208,6 +208,10 @@ class Plot(IDManagerMixin):
         The cells or materials to plot
     mask_background : Iterable of int or str
         Color to apply to all cells/materials not listed in mask_components
+    show_overlaps : bool
+        Inidicate whether or not overlapping regions are shown
+    overlap_color : Iterable of int or str
+        Color to apply to overlapping regions
     colors : dict
         Dictionary indicating that certain cells/materials (keys) should be
         displayed with a particular color.
@@ -236,6 +240,8 @@ class Plot(IDManagerMixin):
         self._background = None
         self._mask_components = None
         self._mask_background = None
+        self._show_overlaps = False
+        self._overlap_color = None
         self._colors = {}
         self._level = None
         self._meshlines = None
@@ -284,6 +290,14 @@ class Plot(IDManagerMixin):
     def mask_background(self):
         return self._mask_background
 
+    @property
+    def show_overlaps(self):
+        return self._show_overlaps
+
+    @property
+    def overlap_color(self):
+        return self._overlap_color
+
     @property
     def colors(self):
         return self._colors
@@ -343,15 +357,7 @@ class Plot(IDManagerMixin):
 
     @background.setter
     def background(self, background):
-        cv.check_type('plot background', background, Iterable)
-        if isinstance(background, str):
-            if background.lower() not in _SVG_COLORS:
-                raise ValueError("'{}' is not a valid color.".format(background))
-        else:
-            cv.check_length('plot background', background, 3)
-            for rgb in background:
-                cv.check_greater_than('plot background', rgb, 0, True)
-                cv.check_less_than('plot background', rgb, 256)
+        self._check_color('plot background', background)
         self._background = background
 
     @colors.setter
@@ -359,17 +365,7 @@ class Plot(IDManagerMixin):
         cv.check_type('plot colors', colors, Mapping)
         for key, value in colors.items():
             cv.check_type('plot color key', key, (openmc.Cell, openmc.Material))
-            cv.check_type('plot color value', value, Iterable)
-            if isinstance(value, str):
-                if value.lower() not in _SVG_COLORS:
-                    raise ValueError("'{}' is not a valid color.".format(value))
-            else:
-                cv.check_length('plot color (RGB)', value, 3)
-                for component in value:
-                    cv.check_type('RGB component', component, Real)
-                    cv.check_greater_than('RGB component', component, 0, True)
-                    cv.check_less_than('RGB component', component, 255, True)
-
+            self._check_color('plot color value', value)
         self._colors = colors
 
     @mask_components.setter
@@ -380,17 +376,20 @@ class Plot(IDManagerMixin):
 
     @mask_background.setter
     def mask_background(self, mask_background):
-        cv.check_type('plot mask background', mask_background, Iterable)
-        if isinstance(mask_background, str):
-            if mask_background.lower() not in _SVG_COLORS:
-                raise ValueError("'{}' is not a valid color.".format(mask_background))
-        else:
-            cv.check_length('plot mask_background', mask_background, 3)
-            for rgb in mask_background:
-                cv.check_greater_than('plot mask background', rgb, 0, True)
-                cv.check_less_than('plot mask background', rgb, 256)
+        self._check_color('plot mask background', mask_background)
         self._mask_background = mask_background
 
+    @show_overlaps.setter
+    def show_overlaps(self, show_overlaps):
+        cv.check_type('Show overlaps flag for Plot ID="{}"'.format(self.id),
+                      show_overlaps, bool)
+        self._show_overlaps = show_overlaps
+
+    @overlap_color.setter
+    def overlap_color(self, overlap_color):
+        self._check_color('plot overlap color', overlap_color)
+        self._overlap_color = overlap_color
+
     @level.setter
     def level(self, plot_level):
         cv.check_type('plot level', plot_level, Integral)
@@ -421,15 +420,23 @@ class Plot(IDManagerMixin):
                                   0, equality=True)
 
         if 'color' in meshlines:
-            cv.check_type('plot meshlines color', meshlines['color'], Iterable,
-                          Integral)
-            cv.check_length('plot meshlines color', meshlines['color'], 3)
-            for rgb in meshlines['color']:
-                cv.check_greater_than('plot meshlines color', rgb, 0, True)
-                cv.check_less_than('plot meshlines color', rgb, 256)
+            self._check_color('plot meshlines color', meshlines['color'])
 
         self._meshlines = meshlines
 
+    @staticmethod
+    def _check_color(err_string, color):
+        cv.check_type(err_string, color, Iterable)
+        if isinstance(color, str):
+            if color.lower() not in _SVG_COLORS:
+                raise ValueError("'{}' is not a valid color.".format(color))
+        else:
+            cv.check_length(err_string, color, 3)
+            for rgb in color:
+                cv.check_type(err_string, rgb, Real)
+                cv.check_greater_than('RGB component', rgb, 0, True)
+                cv.check_less_than('RGB component', rgb, 256)
+
     def __repr__(self):
         string = 'Plot\n'
         string += '{: <16}=\t{}\n'.format('\tID', self._id)
@@ -446,6 +453,8 @@ class Plot(IDManagerMixin):
                                             self._mask_components)
         string += '{: <16}=\t{}\n'.format('\tMask background',
                                             self._mask_background)
+        string += '{: <16}=\t{}\n'.format('\Overlap Color',
+                                            self._overlap_color)
         string += '{: <16}=\t{}\n'.format('\tColors', self._colors)
         string += '{: <16}=\t{}\n'.format('\tLevel', self._level)
         string += '{: <16}=\t{}\n'.format('\tMeshlines', self._meshlines)
@@ -635,6 +644,18 @@ class Plot(IDManagerMixin):
                 subelement.set("background", ' '.join(
                     str(x) for x in color))
 
+        if self._show_overlaps:
+            subelement = ET.SubElement(element, "show_overlaps")
+            subelement.text = "true"
+
+            if self._overlap_color is not None:
+                color = self._overlap_color
+                if isinstance(color, str):
+                    color = _SVG_COLORS[color.lower()]
+                subelement = ET.SubElement(element, "overlap_color")
+                subelement.text = ' '.join(str(x) for x in color)
+
+
         if self._level is not None:
             subelement = ET.SubElement(element, "level")
             subelement.text = str(self._level)

src/geometry.cpp

@@ -32,7 +32,7 @@ std::vector<int64_t> overlap_check_count;
 // Non-member functions
 //==============================================================================
 
-bool check_cell_overlap(Particle* p)
+bool check_cell_overlap(Particle* p, bool error)
 {
   int n_coord = p->n_coord_;
 
@@ -45,11 +45,14 @@ bool check_cell_overlap(Particle* p)
       Cell& c = *model::cells[index_cell];
       if (c.contains(p->coord_[j].r, p->coord_[j].u, p->surface_)) {
         if (index_cell != p->coord_[j].cell) {
-          std::stringstream err_msg;
-          err_msg << "Overlapping cells detected: " << c.id_ << ", "
-                  << model::cells[p->coord_[j].cell]->id_ << " on universe "
-                  << univ.id_;
-          fatal_error(err_msg);
+          if (error) {
+            std::stringstream err_msg;
+            err_msg << "Overlapping cells detected: " << c.id_ << ", "
+                    << model::cells[p->coord_[j].cell]->id_ << " on universe "
+                    << univ.id_;
+            fatal_error(err_msg);
+          }
+          return true;
         }
         ++model::overlap_check_count[index_cell];
       }

src/plot.cpp

@@ -28,9 +28,10 @@ namespace openmc {
 // Constants
 //==============================================================================
 
-const RGBColor WHITE {255, 255, 255};
+
 constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level
 constexpr int32_t NOT_FOUND {-2};
+constexpr int32_t OVERLAP {-3};
 
 IdData::IdData(size_t h_res, size_t v_res)
   : data_({v_res, h_res, 2}, NOT_FOUND)
@@ -49,6 +50,10 @@ IdData::set_value(size_t y, size_t x, const Particle& p, int level) {
   }
 }
 
+void IdData::set_overlap(size_t y, size_t x) {
+  xt::view(data_, y, x, xt::all()) = OVERLAP;
+}
+
 PropertyData::PropertyData(size_t h_res, size_t v_res)
   : data_({v_res, h_res, 2}, NOT_FOUND)
 { }
@@ -63,6 +68,10 @@ PropertyData::set_value(size_t y, size_t x, const Particle& p, int level) {
   }
 }
 
+void PropertyData::set_overlap(size_t y, size_t x) {
+  data_(y, x) = OVERLAP;
+}
+
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -143,6 +152,10 @@ void create_ppm(Plot pl)
       auto id = ids.data_(y, x, pl.color_by_);
       // no setting needed if not found
       if (id == NOT_FOUND) { continue; }
+      if (id == OVERLAP) {
+        data(x,y) = pl.overlap_color_;
+        continue;
+      }
       if (PlotColorBy::cells == pl.color_by_) {
         data(x,y) = pl.colors_[model::cell_map[id]];
       } else if (PlotColorBy::mats == pl.color_by_) {
@@ -276,9 +289,6 @@ Plot::set_bg_color(pugi::xml_node plot_node)
               << id_;
       fatal_error(err_msg);
     }
-  } else {
-    // default to a white background
-    not_found_ = WHITE;
   }
 }
 
@@ -388,6 +398,10 @@ Plot::set_default_colors(pugi::xml_node plot_node)
 
   for (auto& c : colors_) {
     c = random_color();
+    // make sure we don't interfere with some default colors
+    while (c == RED || c == WHITE) {
+      c = random_color();
+    }
   }
 }
 
@@ -638,8 +652,39 @@ Plot::set_mask(pugi::xml_node plot_node)
   }
 }
 
+void Plot::set_overlap_color(pugi::xml_node plot_node) {
+  color_overlaps_ = false;
+  if (check_for_node(plot_node, "show_overlaps")) {
+    color_overlaps_ = get_node_value_bool(plot_node, "show_overlaps");
+    // check for custom overlap color
+    if (check_for_node(plot_node, "overlap_color")) {
+      if (!color_overlaps_) {
+        std::stringstream wrn_msg;
+        wrn_msg << "Overlap color specified in plot " << id_
+                << " but overlaps won't be shown.";
+        warning(wrn_msg);
+      }
+      std::vector<int> olap_clr = get_node_array<int>(plot_node, "overlap_color");
+      if (olap_clr.size() == 3) {
+        overlap_color_ = olap_clr;
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Bad overlap RGB in plot " << id_;
+        fatal_error(err_msg);
+      }
+    }
+  }
+
+  // make sure we allocate the vector for counting overlap checks if
+  // they're going to be plotted
+  if (color_overlaps_ && settings::run_mode == RUN_MODE_PLOTTING) {
+    settings::check_overlaps = true;
+    model::overlap_check_count.resize(model::cells.size(), 0);
+  }
+}
+
 Plot::Plot(pugi::xml_node plot_node)
-  : index_meshlines_mesh_{-1}
+  : index_meshlines_mesh_{-1}, overlap_color_{RED}
 {
   set_id(plot_node);
   set_type(plot_node);
@@ -653,6 +698,7 @@ Plot::Plot(pugi::xml_node plot_node)
   set_user_colors(plot_node);
   set_meshlines(plot_node);
   set_mask(plot_node);
+  set_overlap_color(plot_node);
 } // End Plot constructor
 
 //==============================================================================
@@ -849,6 +895,7 @@ void create_voxel(Plot pl)
   pltbase.basis_ = PlotBasis::xy;
   pltbase.pixels_ = pl.pixels_;
   pltbase.level_ = -1; // all universes for voxel files
+  pltbase.color_overlaps_ = pl.color_overlaps_;
 
   ProgressBar pb;
   for (int z = 0; z < pl.pixels_[2]; z++) {
@@ -922,6 +969,10 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
     return OPENMC_E_INVALID_ARGUMENT;
   }
 
+  if (plt->color_overlaps_ && model::overlap_check_count.size() == 0) {
+    model::overlap_check_count.resize(model::cells.size());
+  }
+
   auto ids = plt->get_map<IdData>();
 
   // write id data to array
@@ -938,6 +989,10 @@ extern "C" int openmc_property_map(const void* plot, double* data_out) {
     return OPENMC_E_INVALID_ARGUMENT;
   }
 
+  if (plt->color_overlaps_ && model::overlap_check_count.size() == 0) {
+    model::overlap_check_count.resize(model::cells.size());
+  }
+
   auto props = plt->get_map<PropertyData>();
 
   // write id data to array

tests/regression_tests/plot_overlaps/geometry.xml

@@ -0,0 +1,14 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <surface id="1" type="z-cylinder" coeffs="0 0 2" />
+  <surface id="20" type="z-cylinder" coeffs="0 0 5.5" />
+  <surface id="2" type="z-cylinder" coeffs="0 0 5" />
+  <surface id="3" type="z-cylinder" coeffs="0 0 10" boundary="vacuum" />
+  <surface id="4" type="y-plane" coeffs="5" boundary="vacuum" />
+
+  <cell id="1" material="1" region="  -1 -4" />
+  <cell id="2" material="3" region="1 -20 -4" />
+  <cell id="3" material="2" region="2 -3 -4" />
+
+</geometry>

tests/regression_tests/plot_overlaps/materials.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+  <material id="2">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+  <material id="3">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+</materials>

tests/regression_tests/plot_overlaps/plots.xml

@@ -0,0 +1,35 @@
+<?xml version="1.0"?>
+<plots>
+
+  <plot id="1" basis="xy">
+    <origin>0. 0. 0.</origin>
+    <width>25 25</width>
+    <pixels>200 200</pixels>
+    <color id="1" rgb="255 0 0" /> <!-- Red -->
+    <meshlines meshtype="entropy" linewidth="0" />
+    <show_overlaps>true</show_overlaps>
+  </plot>
+
+  <plot id="2" basis="xz">
+    <origin>0. 0. 0.</origin>
+    <width>25 25</width>
+    <pixels>200 200</pixels>
+    <mask components="1 3" background="255 255 255" />
+    <show_overlaps>true</show_overlaps>
+    <overlap_color>255 211 0</overlap_color>
+  </plot>
+
+  <plot id="3" basis="yz" color_by="material">
+    <origin>0. 0. 0.</origin>
+    <width>25 25</width>
+    <pixels>200 200</pixels>
+    <background>0 0 0</background>
+  </plot>
+
+  <plot id="4" type="voxel">
+    <pixels>100 100 10</pixels>
+    <origin>0. 0. 0.</origin>
+    <width>20 20 10</width>
+  </plot>
+
+</plots>

tests/regression_tests/plot_overlaps/results_true.dat

@@ -0,0 +1 @@
+e67f7737b49af347a617bfdeebebc3288bd85050f33c9208403f3dfb4625a42f63f3396ecb2517ed9214c3e4e68c9038a8f9a7b3e27a414565c5580827dbb6e6

tests/regression_tests/plot_overlaps/settings.xml

@@ -0,0 +1,13 @@
+<?xml version="1.0"?>
+<settings>
+
+  <run_mode>plot</run_mode>
+
+  <mesh id="1">
+    <dimension>5 4 3</dimension>
+    <lower_left>-10 -10 -10</lower_left>
+    <upper_right>10 10 10</upper_right>
+  </mesh>
+  <entropy_mesh>1</entropy_mesh>
+
+</settings>

tests/regression_tests/plot_overlaps/test.py

@@ -0,0 +1,61 @@
+import glob
+import hashlib
+import os
+
+import h5py
+import openmc
+
+from tests.testing_harness import TestHarness
+from tests.regression_tests import config
+
+
+class PlotTestHarness(TestHarness):
+    """Specialized TestHarness for running OpenMC plotting tests."""
+    def __init__(self, plot_names):
+        super().__init__(None)
+        self._plot_names = plot_names
+
+    def _run_openmc(self):
+        openmc.plot_geometry(openmc_exec=config['exe'])
+
+    def _test_output_created(self):
+        """Make sure *.ppm has been created."""
+        for fname in self._plot_names:
+            assert os.path.exists(fname), 'Plot output file does not exist.'
+
+    def _cleanup(self):
+        super()._cleanup()
+        for fname in self._plot_names:
+            if os.path.exists(fname):
+                os.remove(fname)
+
+    def _get_results(self):
+        """Return a string hash of the plot files."""
+        outstr = bytes()
+
+        for fname in self._plot_names:
+            if fname.endswith('.ppm'):
+                # Add PPM output to results
+                with open(fname, 'rb') as fh:
+                    outstr += fh.read()
+            elif fname.endswith('.h5'):
+                # Add voxel data to results
+                with h5py.File(fname, 'r') as fh:
+                    outstr += fh.attrs['filetype']
+                    outstr += fh.attrs['num_voxels'].tostring()
+                    outstr += fh.attrs['lower_left'].tostring()
+                    outstr += fh.attrs['voxel_width'].tostring()
+                    outstr += fh['data'].value.tostring()
+
+        # Hash the information and return.
+        sha512 = hashlib.sha512()
+        sha512.update(outstr)
+        outstr = sha512.hexdigest()
+
+        return outstr
+
+
+def test_plot_overlap():
+    harness = PlotTestHarness(('plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm',
+                               'plot_4.h5'))
+    harness.main()

tests/unit_tests/test_plots.py

@@ -24,6 +24,10 @@ def myplot():
     plot.mask_background = (255, 255, 255)
     plot.mask_background = 'white'
 
+    plot.overlap_color = (255, 211, 0)
+    plot.overlap_color = 'yellow'
+    plot.show_overlaps = True
+
     plot.level = 1
     plot.meshlines = {
         'type': 'tally',