From 61cab3009ab817fe76f5eb9dcf36cb49603ea919 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Fri, 17 Mar 2023 13:20:16 -0400 Subject: [PATCH] PR comments addressed --- docs/source/usersguide/plots.rst | 23 +++- include/openmc/plot.h | 43 +++---- openmc/_xml.py | 22 ++++ openmc/plots.py | 200 ++++++++++++++----------------- src/plot.cpp | 7 +- 5 files changed, 157 insertions(+), 138 deletions(-) diff --git a/docs/source/usersguide/plots.rst b/docs/source/usersguide/plots.rst index ae386668e..35a2aff61 100644 --- a/docs/source/usersguide/plots.rst +++ b/docs/source/usersguide/plots.rst @@ -160,7 +160,7 @@ with OpenMC, before exporting to plots.xml. thisp.set_transparent(geometry) thisp.xs[fuel] = 1.0 thisp.xs[iron] = 1.0 - thisp.wireframe_regions = [fuel] + thisp.wireframe_domains = [fuel] thisp.wireframe_thickness = 2 plot_file.append(thisp) @@ -182,17 +182,30 @@ can be used to make all materials in the geometry transparent. From there, individual material or cell opacities can be tuned to produce the desired result. +Two camera projections are available when using these plots, perspective +and orthographic. The default, perspective projection, +is a cone of rays passing through each pixel which radiate from the camera +position and span the field of view in the x and y positions. The horizontal +field of view can be set with the :attr: `ProjectionPlot.horizontal_field_of_view` attribute, +which is to be specified in units of degrees. The field of view only influences +behavior in perspective projection mode. + +In the orthographic projection, rays follow the same angle but originate from +different points. The horizontal width of this plane of ray starting points +may be set with the :attr: `ProjectionPlot.orthographic_width` element. If this element +is nonzero, the orthographic projection is employed. Left to its default value +of zero, the perspective projection is employed. + Lastly, projection plots come packaged with wireframe generation that can target either all surface/cell/material boundaries in the geometry, or only wireframing around specific regions. In the above example, we have set only the fuel region from the hexagonal lattice example to have -a wireframe drawn around it. The :attr:`ProjectionPlot.wireframe_thickness` +a wireframe drawn around it. This is accomplished by setting the +:attr: `ProjectionPlot.wireframe_domains`, which may be set to either material +IDs or cell IDs. The :attr:`ProjectionPlot.wireframe_thickness` attribute sets the wireframe thickness in units of pixels. .. note:: When setting specific material or cell regions to have wireframes drawn around them, the plot must be colored by materials if wireframing around specific materials and similarly colored by cell instance if wireframing around specific cells. - - - diff --git a/include/openmc/plot.h b/include/openmc/plot.h index b09727835..fc78c250b 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -28,7 +28,7 @@ class PlottableInterface; namespace model { extern std::unordered_map plot_map; //!< map of plot ids to index -extern std::vector> +extern vector> plots; //!< Plot instance container extern uint64_t plotter_seed; // Stream index used by the plotter @@ -285,6 +285,19 @@ private: void set_wireframe_ids(pugi::xml_node node); void set_wireframe_color(pugi::xml_node node); + /* If starting the particle from outside the geometry, we have to + * find a distance to the boundary in a non-standard surface intersection + * check. It's an exhaustive search over surfaces in the top-level universe. + */ + static int advance_to_boundary_from_void(Particle& p); + + /* Checks if a vector of two TrackSegments is equivalent. We define this + * to mean not having matching intersection lengths, but rather having + * a matching sequence of surface/cell/material intersections. + */ + bool trackstack_equivalent(const vector& track1, + const vector& track2) const; + /* Used for drawing wireframe and colors. We record the list of * surface/cell/material intersections and the corresponding lengths as a ray * traverses the geometry, then color by iterating in reverse. @@ -303,13 +316,18 @@ private: {} }; + // Max intersections before we assume ray tracing is caught in an infinite + // loop: + static const int MAX_INTERSECTIONS = 1000000; + std::array pixels_; // pixel dimension of resulting image double horizontal_field_of_view_ {70.0}; // horiz. f.o.v. in degrees Position camera_position_; // where camera is Position look_at_; // point camera is centered looking at Direction up_ {0.0, 0.0, 1.0}; // which way is up - std::vector - wireframe_ids_; // which color IDs should be wireframed. If empty, all + + // which color IDs should be wireframed. If empty, all cells are wireframed. + vector wireframe_ids_; /* The horizontal thickness, if using an orthographic projection. * If set to zero, we assume using a perspective projection. @@ -320,24 +338,7 @@ private: int wireframe_thickness_ {1}; RGBColor wireframe_color_ {BLACK}; // wireframe color - std::vector - xs_; // macro cross section values for cell volume rendering - - /* If starting the particle from outside the geometry, we have to - * find a distance to the boundary in a non-standard surface intersection - * check. It's an exhaustive search over surfaces in the top-level universe. - */ - static int advance_to_boundary_from_void(Particle& p); - - /* Checks if a vector of two TrackSegments is equivalent. We define this - * to mean not having matching intersection lengths, but rather having - * a matching sequence of surface/cell/material intersections. - */ - bool trackstack_equivalent(const std::vector& track1, - const std::vector& track2) const; - - // Closed form 3x3 matrix inversion - static std::vector invert_matrix(const std::vector& input); + vector xs_; // macro cross section values for cell volume rendering }; //=============================================================================== diff --git a/openmc/_xml.py b/openmc/_xml.py index 6799a4e2d..78bcbc264 100644 --- a/openmc/_xml.py +++ b/openmc/_xml.py @@ -79,3 +79,25 @@ def reorder_attributes(root): attribs = sorted(attrib.items()) attrib.clear() attrib.update(attribs) + + +def get_elem_tuple(elem, name, dtype=int): + '''Helper function to get a tuple of values from an elem + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element that should contain a tuple + name : str + Name of the subelement to obtain tuple from + dtype : data-type + The type of each element in the tuple + + Returns + ------- + tuple of dtype + Data read from the tuple + ''' + subelem = elem.find(name) + if subelem is not None: + return tuple([dtype(x) for x in subelem.text.split()]) diff --git a/openmc/plots.py b/openmc/plots.py index 91a03a5cd..508fb6c23 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -7,7 +7,7 @@ import numpy as np import openmc import openmc.checkvalue as cv -from ._xml import clean_indentation, reorder_attributes +from ._xml import clean_indentation, reorder_attributes, get_elem_tuple from .mixin import IDManagerMixin @@ -177,27 +177,6 @@ def _get_plot_image(plot, cwd): return Image(str(png_file)) -def get_tuple(elem, name, dtype=int): - '''Helper function to get a tuple of values - - Parameters - ---------- - elem : xml.etree.ElementTree.Element - XML element that should contain a tuple - name : str - Name of the subelement to obtain tuple from - dtype : data-type - The type of each element in the tuple - - Returns - ------- - tuple of dtype - Data read from the tuple - ''' - subelem = elem.find(name) - if subelem is not None: - return tuple([dtype(x) for x in subelem.text.split()]) - class PlotBase(IDManagerMixin): """ @@ -411,6 +390,49 @@ class PlotBase(IDManagerMixin): for domain in domains: self.colors[domain] = np.random.randint(0, 256, (3,)) + def to_xml_element(self): + """Save common plot attributes to XML element + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing plot data + + """ + + element = ET.Element("plot") + element.set("id", str(self._id)) + if self._filename is not None: + element.set("filename", self._filename) + element.set("color_by", self._color_by) + + subelement = ET.SubElement(element, "pixels") + subelement.text = ' '.join(map(str, self._pixels)) + + if self._background is not None: + subelement = ET.SubElement(element, "background") + color = self._background + if isinstance(color, str): + color = _SVG_COLORS[color.lower()] + subelement.text = ' '.join(str(x) for x in color) + + if self._mask_components is not None: + subelement = ET.SubElement(element, "mask") + subelement.set("components", ' '.join( + str(get_id(d)) for d in self._mask_components)) + color = self._mask_background + if color is not None: + if isinstance(color, str): + color = _SVG_COLORS[color.lower()] + subelement.set("background", ' '.join( + str(x) for x in color)) + + if self._level is not None: + subelement = ET.SubElement(element, "level") + subelement.text = str(self._level) + + return element + class Plot(PlotBase): '''Definition of a finite region of space to be plotted. @@ -651,11 +673,7 @@ class Plot(PlotBase): """ - element = ET.Element("plot") - element.set("id", str(self._id)) - if self._filename is not None: - element.set("filename", self._filename) - element.set("color_by", self._color_by) + element = super().to_xml_element() element.set("type", self._type) if self._type == 'slice': @@ -667,16 +685,6 @@ class Plot(PlotBase): subelement = ET.SubElement(element, "width") subelement.text = ' '.join(map(str, self._width)) - subelement = ET.SubElement(element, "pixels") - subelement.text = ' '.join(map(str, self._pixels)) - - if self._background is not None: - subelement = ET.SubElement(element, "background") - color = self._background - if isinstance(color, str): - color = _SVG_COLORS[color.lower()] - subelement.text = ' '.join(str(x) for x in color) - # Helper function that returns the domain ID given either a # Cell/Material object or the domain ID itself def get_id(domain): @@ -691,17 +699,6 @@ class Plot(PlotBase): color = _SVG_COLORS[color.lower()] subelement.set("rgb", ' '.join(str(x) for x in color)) - if self._mask_components is not None: - subelement = ET.SubElement(element, "mask") - subelement.set("components", ' '.join( - str(get_id(d)) for d in self._mask_components)) - color = self._mask_background - if color is not None: - if isinstance(color, str): - color = _SVG_COLORS[color.lower()] - subelement.set("background", ' '.join( - str(x) for x in color)) - if self._show_overlaps: subelement = ET.SubElement(element, "show_overlaps") subelement.text = "true" @@ -714,10 +711,6 @@ class Plot(PlotBase): 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) - if self._meshlines is not None: subelement = ET.SubElement(element, "meshlines") subelement.set("meshtype", self._meshlines['type']) @@ -754,10 +747,10 @@ class Plot(PlotBase): plot.type = elem.get("type") plot.basis = elem.get("basis") - plot.origin = get_tuple(elem, "origin", float) - plot.width = get_tuple(elem, "width", float) - plot.pixels = get_tuple(elem, "pixels") - plot._background = get_tuple(elem, "background") + plot.origin = get_elem_tuple(elem, "origin", float) + plot.width = get_elem_tuple(elem, "width", float) + plot.pixels = get_elem_tuple(elem, "pixels") + plot._background = get_elem_tuple(elem, "background") # Set plot colors colors = {} @@ -778,7 +771,7 @@ class Plot(PlotBase): overlap_elem = elem.find("show_overlaps") if overlap_elem is not None: plot.show_overlaps = (overlap_elem.text in ('true', '1')) - overlap_color = get_tuple(elem, "overlap_color") + overlap_color = get_elem_tuple(elem, "overlap_color") if overlap_color is not None: plot.overlap_color = overlap_color @@ -857,35 +850,37 @@ class ProjectionPlot(PlotBase): Attributes ---------- horizontal_field_of_view : float - Field of view horizontally, in units of degrees. + Field of view horizontally, in units of degrees, defaults to 70. camera_position : tuple or list of ndarray - Position of the camera in 3D space + Position of the camera in 3D space. Defaults to (1, 0, 0). look_at : tuple or list of ndarray - The center of the camera's image points to this place in 3D space + The center of the camera's image points to this place in 3D space. + Set to (0, 0, 0) by default. up : tuple or list of ndarray Which way is up for the camera. Must not be parallel to the - line between look_at and camera_position + line between look_at and camera_position. Set to (0, 0, 1) by default. orthographic_width : float If set to a nonzero value, an orthographic projection is used. All rays traced from the orthographic pixel array travel in the same direction. The width of the starting array must be specified, unlike with the default perspective projection. The height of the array is deduced from the ratio of pixel dimensions for the image. - + Defaults to zero, i.e. using perspective projection. wireframe_thickness : int Line thickness employed for drawing wireframes around cells or - material regions. Can be set to zero for no wireframes + material regions. Can be set to zero for no wireframes at all. + Defaults to one pixel. wireframe_color : tuple of ints - RGB color of the wireframe lines - wireframe_regions : iterable of either Material or Cells + RGB color of the wireframe lines. Defaults to black. + wireframe_domains : iterable of either Material or Cells If provided, the wireframe is only drawn around these. If color_by is by material, it must be a list of materials, else cells. - xs : dict A mapping from cell/material IDs to floats. The floating point values are macroscopic cross sections influencing the volume rendering opacity of each geometric region. Zero corresponds to perfect transparency, and - infinity equivalent to opaque. + infinity equivalent to opaque. These must be set by the user, but default + values can be obtained using the set_transparent method. """ def __init__(self, plot_id=None, name=''): @@ -898,7 +893,7 @@ class ProjectionPlot(PlotBase): self._orthographic_width = 0.0 self._wireframe_thickness = 1 self._wireframe_color = _SVG_COLORS['black'] - self._wireframe_regions = [] + self._wireframe_domains = [] self._xs = {} @property @@ -930,8 +925,8 @@ class ProjectionPlot(PlotBase): return self._wireframe_color @property - def wireframe_regions(self): - return self._wireframe_regions + def wireframe_domains(self): + return self._wireframe_domains @property def xs(self): @@ -966,6 +961,7 @@ class ProjectionPlot(PlotBase): @orthographic_width.setter def orthographic_width(self, orthographic_width): cv.check_type('plot orthographic width', orthographic_width, Real) + assert orthographic_width >= 0.0 self._orthographic_width = orthographic_width @wireframe_thickness.setter @@ -979,9 +975,9 @@ class ProjectionPlot(PlotBase): self._check_color('plot wireframe color', wireframe_color) self._wireframe_color = wireframe_color - @wireframe_regions.setter - def wireframe_regions(self, wireframe_regions): - for region in wireframe_regions: + @wireframe_domains.setter + def wireframe_domains(self, wireframe_domains): + for region in wireframe_domains: if self._color_by == 'material': if not isinstance(region, openmc.Material): raise Exception('Must provide a list of materials for \ @@ -990,7 +986,7 @@ class ProjectionPlot(PlotBase): if not isinstance(region, openmc.Cell): raise Exception('Must provide a list of cells for \ wireframe_region if color_by=cell') - self._wireframe_regions = wireframe_regions + self._wireframe_domains = wireframe_domains @xs.setter def xs(self, xs): @@ -1032,16 +1028,9 @@ class ProjectionPlot(PlotBase): """ - element = ET.Element("plot") - element.set("id", str(self._id)) - if self._filename is not None: - element.set("filename", self._filename) - element.set("color_by", self._color_by) + element = super().to_xml_element() element.set("type", "projection") - subelement = ET.SubElement(element, "pixels") - subelement.text = ' '.join(map(str, self._pixels)) - subelement = ET.SubElement(element, "camera_position") subelement.text = ' '.join(map(str, self._camera_position)) @@ -1057,29 +1046,13 @@ class ProjectionPlot(PlotBase): color = _SVG_COLORS[color.lower()] subelement.text = ' '.join(str(x) for x in color) - if self._wireframe_regions: - id_list = [x.id for x in self._wireframe_regions] + if self._wireframe_domains: + id_list = [x.id for x in self._wireframe_domains] subelement = ET.SubElement(element, "wireframe_ids") subelement.text = ' '.join([str(x) for x in id_list]) - if self._background is not None: - subelement = ET.SubElement(element, "background") - color = self._background - if isinstance(color, str): - color = _SVG_COLORS[color.lower()] - subelement.text = ' '.join(str(x) for x in color) - - if self._mask_components is not None: - subelement = ET.SubElement(element, "mask") - subelement.set("components", ' '.join( - str(d.id) for d in self._mask_components)) - color = self._mask_background - if color is not None: - if isinstance(color, str): - color = _SVG_COLORS[color.lower()] - subelement.set("background", ' '.join( - str(x) for x in color)) - + # note that this differs from the slice plot colors + # in that "xs" must also be specified if self._colors: for domain, color in sorted(self._colors.items(), key=lambda x: x[0].id): @@ -1090,9 +1063,13 @@ class ProjectionPlot(PlotBase): subelement.set("rgb", ' '.join(str(x) for x in color)) subelement.set("xs", str(self._xs[domain])) - if self._level is not None: - subelement = ET.SubElement(element, "level") - subelement.text = str(self._level) + subelement = ET.SubElement(element, "horizontal_field_of_view") + subelement.text = str(self._horizontal_field_of_view) + + # do not need to write if orthographic_width == 0.0 + if self._orthographic_width > 0.0: + subelement = ET.SubElement(element, "orthographic_width") + subelement.text = str(self._orthographic_width) return element @@ -1117,10 +1094,15 @@ class ProjectionPlot(PlotBase): plot.filename = elem.get("filename") plot.color_by = elem.get("color_by") plot.type = "projection" + plot.horizontal_field_of_view = elem.get("horizontal_field_of_view") - plot.pixels = get_tuple(elem, "pixels") - plot.camera_position = get_tuple(elem, "camera_position", float) - plot.look_at = get_tuple(elem, "look_at", float) + tmp = elem.get("horizontal_field_of_view") + if tmp is not None: + self.orthographic_width = float(tmp) + + plot.pixels = get_elem_tuple(elem, "pixels") + plot.camera_position = get_elem_tuple(elem, "camera_position", float) + plot.look_at = get_elem_tuple(elem, "look_at", float) # Attempt to get wireframe thickness. May not be present wireframe_thickness = elem.get("wireframe_thickness") @@ -1135,7 +1117,7 @@ class ProjectionPlot(PlotBase): xs = {} for color_elem in elem.findall("color"): uid = color_elem.get("id") - colors[uid] = get_tuple(color_elem, "rgb") + colors[uid] = get_elem_tuple(color_elem, "rgb") xs[uid] = float(color_elem.get("xs")) # Set masking information diff --git a/src/plot.cpp b/src/plot.cpp index f3262185e..e2102b710 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -97,7 +97,7 @@ void PropertyData::set_overlap(size_t y, size_t x) namespace model { std::unordered_map plot_map; -std::vector> plots; +vector> plots; uint64_t plotter_seed = 1; } // namespace model @@ -1293,7 +1293,6 @@ void ProjectionPlot::create_output() const bool hitsomething = false; bool intersection_found = true; int loop_counter = 0; - const int max_intersections = 1000000; this_line_segments[tid][horiz].clear(); @@ -1342,7 +1341,7 @@ void ProjectionPlot::create_output() const first_surface != -1; // -1 if no surface found } loop_counter++; - if (loop_counter > max_intersections) + if (loop_counter > MAX_INTERSECTIONS) fatal_error("Infinite loop in projection plot"); } @@ -1421,6 +1420,8 @@ void ProjectionPlot::create_output() const for (int j = -wireframe_thickness_ / 2; j < wireframe_thickness_ / 2; ++j) if (i * i + j * j < wireframe_thickness_ * wireframe_thickness_) { + + // Check if wireframe pixel is out of bounds if (horiz + i >= 0 && horiz + i < pixels_[0] && vert + j >= 0 && vert + j < pixels_[1]) data(horiz + i, vert + j) = wireframe_color_;