diff --git a/docs/source/_images/hexlat_anim.gif b/docs/source/_images/hexlat_anim.gif new file mode 100644 index 0000000000..fceedc0104 Binary files /dev/null and b/docs/source/_images/hexlat_anim.gif differ diff --git a/docs/source/usersguide/plots.rst b/docs/source/usersguide/plots.rst index d57917a3df..35a2aff61e 100644 --- a/docs/source/usersguide/plots.rst +++ b/docs/source/usersguide/plots.rst @@ -121,3 +121,91 @@ will depend on the 3D viewer, but should be straightforward. program (Visit, ParaView, etc.) if the number of voxels is large (>10 million or so). Thus if you want an accurate picture that renders smoothly, consider using only one voxel in a certain direction. + +---------------- +Projection Plots +---------------- + +.. image:: ../_images/hexlat_anim.gif + :width: 200px + +The :class:`openmc.ProjectionPlot` class presents an alternative method +of producing 3D visualizations of OpenMC geometries. It was developed to +overcome the primary shortcoming of voxel plots, that an enormous number +of voxels must be employed to capture detailed geometric features. +Projection plots perform volume rendering on material or +cell volumes, with colors specified in the same manner as slice plots. +This is done using the native ray tracing capabilities within OpenMC, +so any geometry in which particles successfully run without overlaps +or leaks will work with projection plots. + +One drawback of projection plots is that particle tracks cannot be overlaid +on them at present. Moreover, checking for overlap regions is not currently possible with projection plots. The image heading this section can +be created by adding the following code to the hexagonal lattice example packaged +with OpenMC, before exporting to plots.xml. + +:: + + r = 5 + import numpy as np + for i in range(100): + phi = 2 * np.pi * i/100 + thisp = openmc.ProjectionPlot(plot_id = 4 + i) + thisp.filename = 'frame%s'%(str(i).zfill(3)) + thisp.look_at = [0, 0, 0] + thisp.camera_position = [r * np.cos(phi), r * np.sin(phi), 6 * np.sin(phi)] + thisp.pixels = [200, 200] + thisp.color_by = 'material' + thisp.colorize(geometry) + thisp.set_transparent(geometry) + thisp.xs[fuel] = 1.0 + thisp.xs[iron] = 1.0 + thisp.wireframe_domains = [fuel] + thisp.wireframe_thickness = 2 + + plot_file.append(thisp) + +This generates a sequence of png files which can be joined to form a gif. +Each image specifies a different camera position using some simple periodic +functions to create a perfectly looped gif. :attr:`ProjectionPlot.look_at` +defines where the camera's centerline should point at. +:attr:`ProjectionPlot.camera_position` similarly defines where the camera +is situated in the universe level we seek to plot. The other settings +resemble those employed by :class:`openmc.Plot`, with the exception of +the :class:`ProjectionPlot.set_transparent` method and :attr:`ProjectionPlot.xs` +dictionary. These are used to control volume rendering of material +volumes. "xs" here stands for cross section, and it defines material +opacities in units of inverse centimeters. Setting this value to a +large number would make a material or cell opaque, and setting it to +zero makes a material transparent. Thus, the :class:`ProjectionPlot.set_transparent` +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. 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/dagmc.h b/include/openmc/dagmc.h index 9930763474..5451474d04 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -168,8 +168,7 @@ private: // Non-member functions //============================================================================== -int32_t next_cell( - DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed); +int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ); } // namespace openmc diff --git a/include/openmc/plot.h b/include/openmc/plot.h index a415b17473..225242c702 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -23,12 +23,13 @@ namespace openmc { // Global variables //=============================================================================== -class Plot; +class PlottableInterface; namespace model { extern std::unordered_map plot_map; //!< map of plot ids to index -extern vector plots; //!< Plot instance container +extern vector> + plots; //!< Plot instance container extern uint64_t plotter_seed; // Stream index used by the plotter @@ -66,6 +67,50 @@ struct RGBColor { // some default colors const RGBColor WHITE {255, 255, 255}; const RGBColor RED {255, 0, 0}; +const RGBColor BLACK {0, 0, 0}; + +/* + * PlottableInterface classes just have to have a unique ID in the plots.xml + * file, and guarantee being able to create output in some way. + */ +class PlottableInterface { +private: + void set_id(pugi::xml_node plot_node); + int id_; // unique plot ID + + void set_bg_color(pugi::xml_node plot_node); + void set_universe(pugi::xml_node plot_node); + void set_default_colors(pugi::xml_node plot_node); + void set_user_colors(pugi::xml_node plot_node); + void set_overlap_color(pugi::xml_node plot_node); + void set_mask(pugi::xml_node plot_node); + +protected: + // Plot output filename, derived classes have logic to set it + std::string path_plot_; + +public: + enum class PlotColorBy { cells = 0, mats = 1 }; + + // Creates the output image named path_plot_ + virtual void create_output() const = 0; + + // Print useful info to the terminal + virtual void print_info() const = 0; + + const std::string& path_plot() const { return path_plot_; } + const int id() const { return id_; } + const int level() const { return level_; } + + // Public color-related data + PlottableInterface(pugi::xml_node plot_node); + int level_; // Universe level to plot + bool color_overlaps_; // Show overlapping cells? + PlotColorBy color_by_; // Plot coloring (cell/material) + RGBColor not_found_ {WHITE}; // Plot background color + RGBColor overlap_color_ {RED}; // Plot overlap color + vector colors_; // Plot colors +}; typedef xt::xtensor ImageData; @@ -93,32 +138,30 @@ struct PropertyData { xt::xtensor data_; //!< 2D array of temperature & density data }; -enum class PlotType { slice = 1, voxel = 2 }; - -enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; - -enum class PlotColorBy { cells = 0, mats = 1 }; - //=============================================================================== // Plot class //=============================================================================== -class PlotBase { + +class SlicePlotBase { public: template T get_map() const; + enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; + // Members public: Position origin_; //!< Plot origin in geometry Position width_; //!< Plot width in geometry PlotBasis basis_; //!< Plot basis (XY/XZ/YZ) array pixels_; //!< Plot size in pixels - bool color_overlaps_; //!< Show overlapping cells? - int level_; //!< Plot universe level + bool slice_color_overlaps_; //!< Show overlapping cells? + int slice_level_ {-1}; //!< Plot universe level +private: }; template -T PlotBase::get_map() const +T SlicePlotBase::get_map() const { size_t width = pixels_[0]; @@ -164,7 +207,7 @@ T PlotBase::get_map() const p.r() = xyz; p.u() = dir; p.coord(0).universe = model::root_universe; - int level = level_; + int level = slice_level_; int j {}; #pragma omp for @@ -182,7 +225,7 @@ T PlotBase::get_map() const if (found_cell) { data.set_value(y, x, p, j); } - if (color_overlaps_ && check_cell_overlap(p, false)) { + if (slice_color_overlaps_ && check_cell_overlap(p, false)) { data.set_overlap(y, x); } } // inner for @@ -192,61 +235,129 @@ T PlotBase::get_map() const return data; } -class Plot : public PlotBase { +// Represents either a voxel or pixel plot +class Plot : public PlottableInterface, public SlicePlotBase { public: - // Constructor - Plot(pugi::xml_node plot); + enum class PlotType { slice = 1, voxel = 2 }; + + Plot(pugi::xml_node plot, PlotType type); - // Methods private: - void set_id(pugi::xml_node plot_node); - void set_type(pugi::xml_node plot_node); void set_output_path(pugi::xml_node plot_node); - void set_bg_color(pugi::xml_node plot_node); void set_basis(pugi::xml_node plot_node); void set_origin(pugi::xml_node plot_node); void set_width(pugi::xml_node plot_node); - void set_universe(pugi::xml_node plot_node); - void set_default_colors(pugi::xml_node plot_node); - 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: - int id_; //!< Plot ID + // Add mesh lines to ImageData + void draw_mesh_lines(ImageData& data) const; + void create_image() const; + void create_voxel() const; + + virtual void create_output() const; + virtual void print_info() const; + 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_ {-1}; //!< Index of the mesh to draw on the plot RGBColor meshlines_color_; //!< Color of meshlines on the plot - RGBColor not_found_ {WHITE}; //!< Plot background color - RGBColor overlap_color_ {RED}; //!< Plot overlap color - vector colors_; //!< Plot colors - std::string path_plot_; //!< Plot output filename +}; + +class ProjectionPlot : public PlottableInterface { + +public: + ProjectionPlot(pugi::xml_node plot); + + virtual void create_output() const; + virtual void print_info() const; + +private: + void set_output_path(pugi::xml_node plot_node); + void set_look_at(pugi::xml_node node); + void set_camera_position(pugi::xml_node node); + void set_field_of_view(pugi::xml_node node); + void set_pixels(pugi::xml_node node); + void set_opacities(pugi::xml_node node); + void set_orthographic_width(pugi::xml_node node); + void set_wireframe_thickness(pugi::xml_node node); + 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. + */ + struct TrackSegment; + 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. + */ + struct TrackSegment { + int id; // material or cell ID (which is being colored) + double length; // length of this track intersection + + /* Recording this allows us to draw edges on the wireframe. For instance + * if two surfaces bound a single cell, it allows drawing that sharp edge + * where the surfaces intersect. + */ + int surface; // last surface ID intersected in this segment + TrackSegment(int id_a, double length_a, int surface_a) + : id(id_a), length(length_a), surface(surface_a) + {} + }; + + // 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 + + // 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. + */ + double orthographic_width_ {0.0}; + + // Thickness of the wireframe lines. Can set to zero for no wireframe. + int wireframe_thickness_ {1}; + + RGBColor wireframe_color_ {BLACK}; // wireframe color + vector xs_; // macro cross section values for cell volume rendering }; //=============================================================================== // Non-member functions //=============================================================================== -//! Add mesh lines to image data of a plot object -//! \param[in] plot object -//! \param[out] image data associated with the plot object -void draw_mesh_lines(Plot const& pl, ImageData& data); - -//! Write a PPM image using a plot object's image data -//! \param[in] plot object -//! \param[out] image data associated with the plot object -void output_ppm(Plot const& pl, const ImageData& data); +/* Write a PPM image + * filename - name of output file + * data - image data to write + */ +void output_ppm(const std::string& filename, const ImageData& data); #ifdef USE_LIBPNG -//! Write a PNG image using a plot object's image data -//! \param[in] plot object -//! \param[out] image data associated with the plot object -void output_png(Plot const& pl, const ImageData& data); +/* Write a PNG image + * filename - name of output file + * data - image data to write + */ +void output_png(const std::string& filename, const ImageData& data); #endif //! Initialize a voxel file @@ -286,14 +397,6 @@ void read_plots_xml(pugi::xml_node root); //! Clear memory void free_memory_plot(); -//! Create an image for a plot object -//! \param[in] plot object -void create_image(Plot const& pl); - -//! Create an hdf5 voxel file for a plot object -//! \param[in] plot object -void create_voxel(Plot const& pl); - //! Create a randomly generated RGB color //! \return RGBColor with random value RGBColor random_color(); diff --git a/include/openmc/position.h b/include/openmc/position.h index 4fde420ccb..0200cda3b7 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -83,6 +83,11 @@ struct Position { return x * other.x + y * other.y + z * other.z; } inline double norm() const { return std::sqrt(x * x + y * y + z * z); } + inline Position cross(Position other) const + { + return {y * other.z - z * other.y, z * other.x - x * other.z, + x * other.y - y * other.x}; + } //! Reflect a direction across a normal vector //! \param[in] other Vector to reflect across diff --git a/openmc/_xml.py b/openmc/_xml.py index 6799a4e2d8..78bcbc264f 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 4b7d58e3bd..c3cd88da18 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 @@ -178,14 +178,8 @@ def _get_plot_image(plot, cwd): return Image(str(png_file)) -class Plot(IDManagerMixin): - """Definition of a finite region of space to be plotted. - - OpenMC is capable of generating two-dimensional slice plots and - three-dimensional voxel plots. Colors that are used in plots can be given as - RGB tuples, e.g. (255, 255, 255) would be white, or by a string indicating a - valid `SVG color `_. - +class PlotBase(IDManagerMixin): + """ Parameters ---------- plot_id : int @@ -199,20 +193,12 @@ class Plot(IDManagerMixin): Unique identifier name : str Name of the plot - width : Iterable of float - Width of the plot in each basis direction pixels : Iterable of int - Number of pixels to use in each basis direction - origin : tuple or list of ndarray - Origin (center) of the plot + Number of pixels to use in each direction filename : Path to write the plot to color_by : {'cell', 'material'} Indicate whether the plot should be colored by cell or by material - type : {'slice', 'voxel'} - The type of the plot - basis : {'xy', 'xz', 'yz'} - The basis directions for the plot background : Iterable of int or str Color of the background mask_components : Iterable of openmc.Cell or openmc.Material or int @@ -230,10 +216,6 @@ class Plot(IDManagerMixin): cell/material). level : int Universe depth to plot at - meshlines : dict - Dictionary defining type, id, linewidth and color of a mesh to be - plotted on top of a plot - """ next_id = 1 @@ -243,13 +225,9 @@ class Plot(IDManagerMixin): # Initialize Plot class attributes self.id = plot_id self.name = name - self._width = [4.0, 4.0] self._pixels = [400, 400] - self._origin = [0., 0., 0.] self._filename = None self._color_by = 'cell' - self._type = 'slice' - self._basis = 'xy' self._background = None self._mask_components = None self._mask_background = None @@ -257,24 +235,15 @@ class Plot(IDManagerMixin): self._overlap_color = None self._colors = {} self._level = None - self._meshlines = None @property def name(self): return self._name - @property - def width(self): - return self._width - @property def pixels(self): return self._pixels - @property - def origin(self): - return self._origin - @property def filename(self): return self._filename @@ -283,14 +252,6 @@ class Plot(IDManagerMixin): def color_by(self): return self._color_by - @property - def type(self): - return self._type - - @property - def basis(self): - return self._basis - @property def background(self): return self._background @@ -319,27 +280,11 @@ class Plot(IDManagerMixin): def level(self): return self._level - @property - def meshlines(self): - return self._meshlines - @name.setter def name(self, name): cv.check_type('plot name', name, str) self._name = name - @width.setter - def width(self, width): - cv.check_type('plot width', width, Iterable, Real) - cv.check_length('plot width', width, 2, 3) - self._width = width - - @origin.setter - def origin(self, origin): - cv.check_type('plot origin', origin, Iterable, Real) - cv.check_length('plot origin', origin, 3) - self._origin = origin - @pixels.setter def pixels(self, pixels): cv.check_type('plot pixels', pixels, Iterable, Integral) @@ -358,16 +303,6 @@ class Plot(IDManagerMixin): cv.check_value('plot color_by', color_by, ['cell', 'material']) self._color_by = color_by - @type.setter - def type(self, plottype): - cv.check_value('plot type', plottype, ['slice', 'voxel']) - self._type = plottype - - @basis.setter - def basis(self, basis): - cv.check_value('plot basis', basis, _BASES) - self._basis = basis - @background.setter def background(self, background): self._check_color('plot background', background) @@ -410,6 +345,181 @@ class Plot(IDManagerMixin): cv.check_greater_than('plot level', plot_level, 0, equality=True) self._level = plot_level + @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(f"'{color}' is not a valid 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) + + # Helper function that returns the domain ID given either a + # Cell/Material object or the domain ID itself + @staticmethod + def _get_id(domain): + return domain if isinstance(domain, Integral) else domain.id + + def colorize(self, geometry, seed=1): + """Generate a color scheme for each domain in the plot. + + This routine may be used to generate random, reproducible color schemes. + The colors generated are based upon cell/material IDs in the geometry. + + Parameters + ---------- + geometry : openmc.Geometry + The geometry for which the plot is defined + seed : Integral + The random number seed used to generate the color scheme + + """ + + cv.check_type('geometry', geometry, openmc.Geometry) + cv.check_type('seed', seed, Integral) + cv.check_greater_than('seed', seed, 1, equality=True) + + # Get collections of the domains which will be plotted + if self.color_by == 'material': + domains = geometry.get_all_materials().values() + else: + domains = geometry.get_all_cells().values() + + # Set the seed for the random number generator + np.random.seed(seed) + + # Generate random colors for each feature + 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") + 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(PlotBase._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. + + OpenMC is capable of generating two-dimensional slice plots, or + three-dimensional voxel or projection plots. Colors that are used in plots can be given as + RGB tuples, e.g. (255, 255, 255) would be white, or by a string indicating a + valid `SVG color `_. + + Parameters + ---------- + plot_id : int + Unique identifier for the plot + name : str + Name of the plot + + Attributes + ---------- + width : Iterable of float + Width of the plot in each basis direction + origin : tuple or list of ndarray + Origin (center) of the plot + type : {'slice', 'voxel'} + The type of the plot + basis : {'xy', 'xz', 'yz'} + The basis directions for the plot + meshlines : dict + Dictionary defining type, id, linewidth and color of a mesh to be + plotted on top of a plot + + ''' + + def __init__(self, plot_id=None, name=''): + super().__init__(plot_id, name) + self._width = [4.0, 4.0] + self._origin = [0., 0., 0.] + self._type = 'slice' + self._basis = 'xy' + self._meshlines = None + + @property + def width(self): + return self._width + + @property + def origin(self): + return self._origin + + @property + def type(self): + return self._type + + @property + def basis(self): + return self._basis + + @property + def meshlines(self): + return self._meshlines + + @width.setter + def width(self, width): + cv.check_type('plot width', width, Iterable, Real) + cv.check_length('plot width', width, 2, 3) + self._width = width + + @origin.setter + def origin(self, origin): + cv.check_type('plot origin', origin, Iterable, Real) + cv.check_length('plot origin', origin, 3) + self._origin = origin + + @type.setter + def type(self, plottype): + cv.check_value('plot type', plottype, ['slice', 'voxel', 'projection']) + self._type = plottype + + @basis.setter + def basis(self, basis): + cv.check_value('plot basis', basis, _BASES) + self._basis = basis + @meshlines.setter def meshlines(self, meshlines): cv.check_type('plot meshlines', meshlines, dict) @@ -429,7 +539,8 @@ class Plot(IDManagerMixin): equality=True) if 'linewidth' in meshlines: - cv.check_type('plot mesh linewidth', meshlines['linewidth'], Integral) + cv.check_type('plot mesh linewidth', + meshlines['linewidth'], Integral) cv.check_greater_than('plot mesh linewidth', meshlines['linewidth'], 0, equality=True) @@ -438,19 +549,6 @@ class Plot(IDManagerMixin): 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(f"'{color}' is not a valid 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) @@ -520,38 +618,6 @@ class Plot(IDManagerMixin): plot.basis = basis return plot - def colorize(self, geometry, seed=1): - """Generate a color scheme for each domain in the plot. - - This routine may be used to generate random, reproducible color schemes. - The colors generated are based upon cell/material IDs in the geometry. - - Parameters - ---------- - geometry : openmc.Geometry - The geometry for which the plot is defined - seed : Integral - The random number seed used to generate the color scheme - - """ - - cv.check_type('geometry', geometry, openmc.Geometry) - cv.check_type('seed', seed, Integral) - cv.check_greater_than('seed', seed, 1, equality=True) - - # Get collections of the domains which will be plotted - if self.color_by == 'material': - domains = geometry.get_all_materials().values() - else: - domains = geometry.get_all_cells().values() - - # Set the seed for the random number generator - np.random.seed(seed) - - # Generate random colors for each feature - for domain in domains: - self.colors[domain] = np.random.randint(0, 256, (3,)) - def highlight_domains(self, geometry, domains, seed=1, alpha=0.5, background='gray'): """Use alpha compositing to highlight one or more domains in the plot. @@ -604,7 +670,7 @@ class Plot(IDManagerMixin): self._colors[domain] = (r, g, b) def to_xml_element(self): - """Return XML representation of the plot + """Return XML representation of the slice/voxel plot Returns ------- @@ -613,11 +679,8 @@ class Plot(IDManagerMixin): """ - element = ET.Element("plot") + element = super().to_xml_element() 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.set("type", self._type) if self._type == 'slice': @@ -629,41 +692,15 @@ class Plot(IDManagerMixin): 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): - return domain if isinstance(domain, Integral) else domain.id - if self._colors: for domain, color in sorted(self._colors.items(), - key=lambda x: get_id(x[0])): + key=lambda x: PlotBase._get_id(x[0])): subelement = ET.SubElement(element, "color") - subelement.set("id", str(get_id(domain))) + subelement.set("id", str(PlotBase._get_id(domain))) if isinstance(color, str): 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" @@ -675,11 +712,6 @@ class Plot(IDManagerMixin): 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) - if self._meshlines is not None: subelement = ET.SubElement(element, "meshlines") subelement.set("meshtype", self._meshlines['type']) @@ -716,37 +748,34 @@ class Plot(IDManagerMixin): plot.type = elem.get("type") plot.basis = elem.get("basis") - # Helper function to get a tuple of values - def get_tuple(elem, name, dtype=int): - subelem = elem.find(name) - if subelem is not None: - return tuple([dtype(x) for x in subelem.text.split()]) - - 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 = {} for color_elem in elem.findall("color"): uid = int(color_elem.get("id")) - colors[uid] = tuple([int(x) for x in color_elem.get("rgb").split()]) + colors[uid] = tuple([int(x) + for x in color_elem.get("rgb").split()]) plot.colors = colors # Set masking information mask_elem = elem.find("mask") if mask_elem is not None: - plot.mask_components = [int(x) for x in mask_elem.get("components").split()] + plot.mask_components = [ + int(x) for x in mask_elem.get("components").split()] background = mask_elem.get("background") if background is not None: - plot.mask_background = tuple([int(x) for x in background.split()]) + plot.mask_background = tuple( + [int(x) for x in background.split()]) # show overlaps 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 @@ -803,12 +832,357 @@ class Plot(IDManagerMixin): return _get_plot_image(self, cwd) +class ProjectionPlot(PlotBase): + """Definition of a camera's view of OpenMC geometry + + Colors are defined in the same manner as the Plot class, but with the addition + of a coloring parameter resembling a macroscopic cross section in units of inverse + centimeters. The volume rendering technique is used to color regions of the model. + An infinite cross section denotes a fully opaque region, and zero represents a + transparent region which will expose the color of the regions behind it. + + The camera projection may either by orthographic or perspective. Perspective + projections are more similar to a pinhole camera, and orthographic projections + preserve parallel lines and distances. + + Parameters + ---------- + plot_id : int + Unique identifier for the plot + name : str + Name of the plot + + Attributes + ---------- + horizontal_field_of_view : float + 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. 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. + 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. 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 at all. + Defaults to one pixel. + wireframe_color : tuple of ints + 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. 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=''): + # Initialize Plot class attributes + super().__init__(plot_id, name) + self._horizontal_field_of_view = 70.0 + self._camera_position = (1.0, 0.0, 0.0) + self._look_at = (0.0, 0.0, 0.0) + self._up = (0.0, 0.0, 1.0) + self._orthographic_width = 0.0 + self._wireframe_thickness = 1 + self._wireframe_color = _SVG_COLORS['black'] + self._wireframe_domains = [] + self._xs = {} + + @property + def horizontal_field_of_view(self): + return self._horizontal_field_of_view + + @property + def camera_position(self): + return self._camera_position + + @property + def look_at(self): + return self._look_at + + @property + def up(self): + return self._up + + @property + def orthographic_width(self): + return self._orthographic_width + + @property + def wireframe_thickness(self): + return self._wireframe_thickness + + @property + def wireframe_color(self): + return self._wireframe_color + + @property + def wireframe_domains(self): + return self._wireframe_domains + + @property + def xs(self): + return self._xs + + @horizontal_field_of_view.setter + def horizontal_field_of_view(self, horizontal_field_of_view): + cv.check_type('plot horizontal field of view', horizontal_field_of_view, + Real) + assert horizontal_field_of_view > 0.0 + assert horizontal_field_of_view < 180.0 + self._horizontal_field_of_view = horizontal_field_of_view + + @camera_position.setter + def camera_position(self, camera_position): + cv.check_type('plot camera position', camera_position, Iterable, Real) + cv.check_length('plot camera position', camera_position, 3) + self._camera_position = camera_position + + @look_at.setter + def look_at(self, look_at): + cv.check_type('plot look at', look_at, Iterable, Real) + cv.check_length('plot look at', look_at, 3) + self._look_at = look_at + + @up.setter + def up(self, up): + cv.check_type('plot up', up, Iterable, Real) + cv.check_length('plot up', up, 3) + self._up = up + + @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 + def wireframe_thickness(self, wireframe_thickness): + cv.check_type('plot wireframe thickness', + wireframe_thickness, Integral) + assert wireframe_thickness >= 0 + self._wireframe_thickness = wireframe_thickness + + @wireframe_color.setter + def wireframe_color(self, wireframe_color): + self._check_color('plot wireframe color', wireframe_color) + self._wireframe_color = wireframe_color + + @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 \ + wireframe_region if color_by=Material') + else: + if not isinstance(region, openmc.Cell): + raise Exception('Must provide a list of cells for \ + wireframe_region if color_by=cell') + self._wireframe_domains = wireframe_domains + + @xs.setter + def xs(self, xs): + cv.check_type('plot xs', xs, Mapping) + for key, value in xs.items(): + cv.check_type('plot xs key', key, (openmc.Cell, openmc.Material)) + cv.check_type('plot xs value', value, Real) + assert value >= 0.0 + self._xs = xs + + def set_transparent(self, geometry): + """Sets all volume rendering XS to zero for the model + + Parameters + ---------- + geometry : openmc.Geometry + The geometry for which the plot is defined + """ + + cv.check_type('geometry', geometry, openmc.Geometry) + + # Get collections of the domains which will be plotted + if self.color_by == 'material': + domains = geometry.get_all_materials().values() + else: + domains = geometry.get_all_cells().values() + + # Generate random colors for each feature + for domain in domains: + self.xs[domain] = 0.0 + + def to_xml_element(self): + """Return XML representation of the projection plot + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing plot data + + """ + + element = super().to_xml_element() + element.set("id", str(self._id)) + element.set("type", "projection") + + subelement = ET.SubElement(element, "camera_position") + subelement.text = ' '.join(map(str, self._camera_position)) + + subelement = ET.SubElement(element, "look_at") + subelement.text = ' '.join(map(str, self._look_at)) + + subelement = ET.SubElement(element, "wireframe_thickness") + subelement.text = str(self._wireframe_thickness) + + subelement = ET.SubElement(element, "wireframe_color") + color = self._wireframe_color + if isinstance(color, str): + color = _SVG_COLORS[color.lower()] + subelement.text = ' '.join(str(x) for x in color) + + 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]) + + # 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): + subelement = ET.SubElement(element, "color") + subelement.set("id", str(domain.id)) + if isinstance(color, str): + color = _SVG_COLORS[color.lower()] + subelement.set("rgb", ' '.join(str(x) for x in color)) + subelement.set("xs", str(self._xs[domain])) + + 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 + + def __repr__(self): + string = 'Projection Plot\n' + string += '{: <16}=\t{}\n'.format('\tID', self._id) + string += '{: <16}=\t{}\n'.format('\tName', self._name) + string += '{: <16}=\t{}\n'.format('\tFilename', self._filename) + string += '{: <16}=\t{}\n'.format('\tHorizontal FOV', + self._horizontal_field_of_view) + string += '{: <16}=\t{}\n'.format('\tOrthographic width', + self._orthographic_width) + string += '{: <16}=\t{}\n'.format('\tWireframe thickness', + self._wireframe_thickness) + string += '{: <16}=\t{}\n'.format('\tWireframe color', + self._wireframe_color) + string += '{: <16}=\t{}\n'.format('\tWireframe domains', + self._wireframe_domains) + string += '{: <16}=\t{}\n'.format('\tCamera position', + self._camera_position) + string += '{: <16}=\t{}\n'.format('\tLook at', self._look_at) + string += '{: <16}=\t{}\n'.format('\tUp', self._up) + string += '{: <16}=\t{}\n'.format('\tPixels', self._pixels) + string += '{: <16}=\t{}\n'.format('\tColor by', self._color_by) + string += '{: <16}=\t{}\n'.format('\tBackground', self._background) + string += '{: <16}=\t{}\n'.format('\tColors', self._colors) + string += '{: <16}=\t{}\n'.format('\tTransparencies', self._xs) + string += '{: <16}=\t{}\n'.format('\tLevel', self._level) + return string + + @classmethod + def from_xml_element(cls, elem): + """Generate plot object from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.ProjectionPlot + ProjectionPlot object + + """ + plot_id = int(elem.get("id")) + plot = cls(plot_id) + if "filename" in elem.keys(): + plot.filename = elem.get("filename") + plot.color_by = elem.get("color_by") + plot.type = "projection" + + horizontal_fov = elem.find("horizontal_field_of_view") + if horizontal_fov is not None: + plot.horizontal_field_of_view = float(horizontal_fov.text) + + tmp = elem.find("orthographic_width") + 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") + if wireframe_thickness: + plot.wireframe_thickness = int(wireframe_thickness) + wireframe_color = elem.get("wireframe_color") + if wireframe_color: + plot.wireframe_color = [int(item) for item in wireframe_color] + + # Set plot colors + colors = {} + xs = {} + for color_elem in elem.findall("color"): + uid = color_elem.get("id") + colors[uid] = get_elem_tuple(color_elem, "rgb") + xs[uid] = float(color_elem.get("xs")) + + # Set masking information + mask_elem = elem.find("mask") + if mask_elem is not None: + mask_components = [int(x) + for x in mask_elem.get("components").split()] + # TODO: set mask components (needs geometry information) + background = mask_elem.get("background") + if background is not None: + plot.mask_background = tuple( + [int(x) for x in background.split()]) + + # Set universe level + level = elem.find("level") + if level is not None: + plot.level = int(level.text) + + return plot + + class Plots(cv.CheckedList): """Collection of Plots used for an OpenMC simulation. This class corresponds directly to the plots.xml input file. It can be - thought of as a normal Python list where each member is a :class:`Plot`. It - behaves like a list as the following example demonstrates: + thought of as a normal Python list where each member is inherits from + :class:`PlotBase`. It behaves like a list as the following example + demonstrates: >>> xz_plot = openmc.Plot() >>> big_plot = openmc.Plot() @@ -819,13 +1193,13 @@ class Plots(cv.CheckedList): Parameters ---------- - plots : Iterable of openmc.Plot - Plots to add to the collection + plots : Iterable of openmc.Plot or openmc.ProjectionPlot + plots to add to the collection """ def __init__(self, plots=None): - super().__init__(Plot, 'plots collection') + super().__init__((Plot, ProjectionPlot), 'plots collection') self._plots_file = ET.Element("plots") if plots is not None: self += plots @@ -835,7 +1209,7 @@ class Plots(cv.CheckedList): Parameters ---------- - plot : openmc.Plot + plot : openmc.Plot or openmc.ProjectionPlot Plot to append """ @@ -873,7 +1247,6 @@ class Plots(cv.CheckedList): for plot in self: plot.colorize(geometry, seed) - def highlight_domains(self, geometry, domains, seed=1, alpha=0.5, background='gray'): """Use alpha compositing to highlight one or more domains in the plot. @@ -925,7 +1298,8 @@ class Plots(cv.CheckedList): # Clean the indentation in the file to be user-readable clean_indentation(self._plots_file) - reorder_attributes(self._plots_file) # TODO: Remove when support is Python 3.8+ + # TODO: Remove when support is Python 3.8+ + reorder_attributes(self._plots_file) return self._plots_file @@ -966,7 +1340,11 @@ class Plots(cv.CheckedList): # Generate each plot plots = cls() for e in elem.findall('plot'): - plots.append(Plot.from_xml_element(e)) + plot_type = e.get('type') + if plot_type == 'projection': + plots.append(ProjectionPlot.from_xml_element(e)) + else: + plots.append(Plot.from_xml_element(e)) return plots @classmethod @@ -987,5 +1365,3 @@ class Plots(cv.CheckedList): tree = ET.parse(path) root = tree.getroot() return cls.from_xml_element(root) - - diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 01ce8ef11e..04c8424aac 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -603,6 +603,7 @@ std::pair DAGCell::distance( // indicate that particle is lost surf_idx = -1; dist = INFINITY; + if (settings::run_mode == RunMode::PLOTTING) return {dist, surf_idx}; if (!dagmc_ptr_->is_implicit_complement(vol) || model::universe_map[dag_univ->id_] == model::root_universe) { std::string material_id = p->material() == MATERIAL_VOID @@ -728,19 +729,22 @@ void check_dagmc_root_univ() } } -int32_t next_cell( - DAGUniverse* dag_univ, DAGCell* cur_cell, DAGSurface* surf_xed) -{ - moab::EntityHandle surf = - surf_xed->dagmc_ptr()->entity_by_index(2, surf_xed->dag_index()); - moab::EntityHandle vol = - cur_cell->dagmc_ptr()->entity_by_index(3, cur_cell->dag_index()); +int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) { + auto surfp = dynamic_cast(model::surfaces[surf - 1].get()); + auto cellp = dynamic_cast(model::cells[curr_cell].get()); + auto univp = static_cast(model::universes[univ].get()); + + moab::EntityHandle surf_handle = + surfp->dagmc_ptr()->entity_by_index(2, surfp->dag_index()); + moab::EntityHandle curr_vol = + cellp->dagmc_ptr()->entity_by_index(3, cellp->dag_index()); moab::EntityHandle new_vol; - cur_cell->dagmc_ptr()->next_vol(surf, vol, new_vol); + moab::ErrorCode rval = cellp->dagmc_ptr()->next_vol(surf_handle, curr_vol, new_vol); + if (rval != moab::MB_SUCCESS) return -1; - return cur_cell->dagmc_ptr()->index_by_handle(new_vol) + - dag_univ->cell_idx_offset_; + return cellp->dagmc_ptr()->index_by_handle(new_vol) + + univp->cell_idx_offset_; } } // namespace openmc @@ -756,8 +760,11 @@ void read_dagmc_universes(pugi::xml_node node) "with DAGMC"); } }; + void check_dagmc_root_univ() {}; +int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ); + } // namespace openmc #endif // DAGMC diff --git a/src/output.cpp b/src/output.cpp index 04e15b82db..666ae15b34 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -221,56 +221,11 @@ void print_plot() if (settings::verbosity < 5) return; - for (auto pl : model::plots) { - // Plot id - fmt::print("Plot ID: {}\n", pl.id_); - // Plot filename - fmt::print("Plot file: {}\n", pl.path_plot_); - // Plot level - fmt::print("Universe depth: {}\n", pl.level_); - - // Plot type - if (PlotType::slice == pl.type_) { - fmt::print("Plot Type: Slice\n"); - } else if (PlotType::voxel == pl.type_) { - fmt::print("Plot Type: Voxel\n"); - } - - // Plot parameters - fmt::print( - "Origin: {} {} {}\n", pl.origin_[0], pl.origin_[1], pl.origin_[2]); - - if (PlotType::slice == pl.type_) { - fmt::print("Width: {:4} {:4}\n", pl.width_[0], pl.width_[1]); - } else if (PlotType::voxel == pl.type_) { - fmt::print( - "Width: {:4} {:4} {:4}\n", pl.width_[0], pl.width_[1], pl.width_[2]); - } - - if (PlotColorBy::cells == pl.color_by_) { - fmt::print("Coloring: Cells\n"); - } else if (PlotColorBy::mats == pl.color_by_) { - fmt::print("Coloring: Materials\n"); - } - - if (PlotType::slice == pl.type_) { - switch (pl.basis_) { - case PlotBasis::xy: - fmt::print("Basis: XY\n"); - break; - case PlotBasis::xz: - fmt::print("Basis: XZ\n"); - break; - case PlotBasis::yz: - fmt::print("Basis: YZ\n"); - break; - } - fmt::print("Pixels: {} {}\n", pl.pixels_[0], pl.pixels_[1]); - } else if (PlotType::voxel == pl.type_) { - fmt::print( - "Voxels: {} {} {}\n", pl.pixels_[0], pl.pixels_[1], pl.pixels_[2]); - } - + for (const auto& pl : model::plots) { + fmt::print("Plot ID: {}\n", pl->id()); + fmt::print("Plot file: {}\n", pl->path_plot()); + fmt::print("Universe depth: {}\n", pl->level()); + pl->print_info(); // prints type-specific plot info fmt::print("\n"); } } diff --git a/src/particle.cpp b/src/particle.cpp index cb959e14d7..1ec95a4ae6 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -433,13 +433,7 @@ void Particle::cross_surface() #ifdef DAGMC // in DAGMC, we know what the next cell should be if (surf->geom_type_ == GeometryType::DAG) { - auto surfp = dynamic_cast(surf); - auto cellp = - dynamic_cast(model::cells[cell_last(n_coord() - 1)].get()); - auto univp = static_cast( - model::universes[coord(n_coord() - 1).universe].get()); - // determine the next cell for this crossing - int32_t i_cell = next_cell(univp, cellp, surfp) - 1; + int32_t i_cell = next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) - 1; // save material and temp material_last() = material(); sqrtkT_last() = sqrtkT(); diff --git a/src/plot.cpp b/src/plot.cpp index fac13a1f4a..1f629e9013 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -14,6 +14,7 @@ #endif #include "openmc/constants.h" +#include "openmc/dagmc.h" #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/geometry.h" @@ -97,7 +98,7 @@ void PropertyData::set_overlap(size_t y, size_t x) namespace model { std::unordered_map plot_map; -vector plots; +vector> plots; uint64_t plotter_seed = 1; } // namespace model @@ -110,20 +111,66 @@ extern "C" int openmc_plot_geometry() { for (auto& pl : model::plots) { - write_message(5, "Processing plot {}: {}...", pl.id_, pl.path_plot_); - - if (PlotType::slice == pl.type_) { - // create 2D image - create_image(pl); - } else if (PlotType::voxel == pl.type_) { - // create voxel file for 3D viewing - create_voxel(pl); - } + write_message(5, "Processing plot {}: {}...", pl->id(), pl->path_plot()); + pl->create_output(); } return 0; } +void Plot::create_output() const +{ + if (PlotType::slice == type_) { + // create 2D image + create_image(); + } else if (PlotType::voxel == type_) { + // create voxel file for 3D viewing + create_voxel(); + } +} + +void Plot::print_info() const +{ + // Plot type + if (PlotType::slice == type_) { + fmt::print("Plot Type: Slice\n"); + } else if (PlotType::voxel == type_) { + fmt::print("Plot Type: Voxel\n"); + } + + // Plot parameters + fmt::print("Origin: {} {} {}\n", origin_[0], origin_[1], origin_[2]); + + if (PlotType::slice == type_) { + fmt::print("Width: {:4} {:4}\n", width_[0], width_[1]); + } else if (PlotType::voxel == type_) { + fmt::print("Width: {:4} {:4} {:4}\n", width_[0], width_[1], width_[2]); + } + + if (PlotColorBy::cells == color_by_) { + fmt::print("Coloring: Cells\n"); + } else if (PlotColorBy::mats == color_by_) { + fmt::print("Coloring: Materials\n"); + } + + if (PlotType::slice == type_) { + switch (basis_) { + case PlotBasis::xy: + fmt::print("Basis: XY\n"); + break; + case PlotBasis::xz: + fmt::print("Basis: XZ\n"); + break; + case PlotBasis::yz: + fmt::print("Basis: YZ\n"); + break; + } + fmt::print("Pixels: {} {}\n", pixels_[0], pixels_[1]); + } else if (PlotType::voxel == type_) { + fmt::print("Voxels: {} {} {}\n", pixels_[0], pixels_[1], pixels_[2]); + } +} + void read_plots_xml() { // Check if plots.xml exists; this is only necessary when the plot runmode is @@ -148,8 +195,26 @@ void read_plots_xml() void read_plots_xml(pugi::xml_node root) { for (auto node : root.children("plot")) { - model::plots.emplace_back(node); - model::plot_map[model::plots.back().id_] = model::plots.size() - 1; + std::string id_string = get_node_value(node, "id", true); + int id = std::stoi(id_string); + if (check_for_node(node, "type")) { + std::string type_str = get_node_value(node, "type", true); + if (type_str == "slice") + model::plots.emplace_back( + std::make_unique(node, Plot::PlotType::slice)); + else if (type_str == "voxel") + model::plots.emplace_back( + std::make_unique(node, Plot::PlotType::voxel)); + else if (type_str == "projection") + model::plots.emplace_back(std::make_unique(node)); + else + fatal_error( + fmt::format("Unsupported plot type '{}' in plot {}", type_str, id)); + + model::plot_map[model::plots.back()->id()] = model::plots.size() - 1; + } else { + fatal_error(fmt::format("Must specify plot type in plot {}", id)); + } } } @@ -159,61 +224,58 @@ void free_memory_plot() model::plot_map.clear(); } -//============================================================================== -// CREATE_IMAGE creates an image based on user input from a plots.xml +// creates an image based on user input from a plots.xml // specification in the PNG/PPM format -//============================================================================== - -void create_image(Plot const& pl) +void Plot::create_image() const { - size_t width = pl.pixels_[0]; - size_t height = pl.pixels_[1]; + size_t width = pixels_[0]; + size_t height = pixels_[1]; - ImageData data({width, height}, pl.not_found_); + ImageData data({width, height}, not_found_); // generate ids for the plot - auto ids = pl.get_map(); + auto ids = get_map(); // assign colors for (size_t y = 0; y < height; y++) { for (size_t x = 0; x < width; x++) { - int idx = pl.color_by_ == PlotColorBy::cells ? 0 : 2; + int idx = color_by_ == PlotColorBy::cells ? 0 : 2; auto id = ids.data_(y, x, idx); // no setting needed if not found if (id == NOT_FOUND) { continue; } if (id == OVERLAP) { - data(x, y) = pl.overlap_color_; + data(x, y) = 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_) { + if (PlotColorBy::cells == color_by_) { + data(x, y) = colors_[model::cell_map[id]]; + } else if (PlotColorBy::mats == color_by_) { if (id == MATERIAL_VOID) { data(x, y) = WHITE; continue; } - data(x, y) = pl.colors_[model::material_map[id]]; + data(x, y) = colors_[model::material_map[id]]; } // color_by if-else } // x for loop } // y for loop // draw mesh lines if present - if (pl.index_meshlines_mesh_ >= 0) { - draw_mesh_lines(pl, data); + if (index_meshlines_mesh_ >= 0) { + draw_mesh_lines(data); } // create image file #ifdef USE_LIBPNG - output_png(pl, data); + output_png(path_plot(), data); #else - output_ppm(pl, data); + output_ppm(path_plot(), data); #endif } -void Plot::set_id(pugi::xml_node plot_node) +void PlottableInterface::set_id(pugi::xml_node plot_node) { // Copy data into plots if (check_for_node(plot_node, "id")) { @@ -229,25 +291,15 @@ void Plot::set_id(pugi::xml_node plot_node) } } -void Plot::set_type(pugi::xml_node plot_node) +// Checks if png or ppm is already present +bool file_extension_present( + const std::string& filename, const std::string& extension) { - // Copy plot type - // Default is slice - type_ = PlotType::slice; - // check type specified on plot node - if (check_for_node(plot_node, "type")) { - std::string type_str = get_node_value(plot_node, "type", true); - // set type using node value - if (type_str == "slice") { - type_ = PlotType::slice; - } else if (type_str == "voxel") { - type_ = PlotType::voxel; - } else { - // if we're here, something is wrong - fatal_error( - fmt::format("Unsupported plot type '{}' in plot {}", type_str, id_)); - } - } + std::string file_extension_if_present = + filename.substr(filename.find_last_of(".") + 1); + if (file_extension_if_present == extension) + return true; + return false; } void Plot::set_output_path(pugi::xml_node plot_node) @@ -258,19 +310,22 @@ void Plot::set_output_path(pugi::xml_node plot_node) if (check_for_node(plot_node, "filename")) { filename = get_node_value(plot_node, "filename"); } else { - filename = fmt::format("plot_{}", id_); + filename = fmt::format("plot_{}", id()); } // add appropriate file extension to name switch (type_) { case PlotType::slice: #ifdef USE_LIBPNG - filename.append(".png"); + if (!file_extension_present(filename, "png")) + filename.append(".png"); #else - filename.append(".ppm"); + if (!file_extension_present(filename, "ppm")) + filename.append(".ppm"); #endif break; case PlotType::voxel: - filename.append(".h5"); + if (!file_extension_present(filename, "h5")) + filename.append(".h5"); break; } @@ -284,7 +339,7 @@ void Plot::set_output_path(pugi::xml_node plot_node) pixels_[1] = pxls[1]; } else { fatal_error( - fmt::format(" must be length 2 in slice plot {}", id_)); + fmt::format(" must be length 2 in slice plot {}", id())); } } else if (PlotType::voxel == type_) { if (pxls.size() == 3) { @@ -293,25 +348,20 @@ void Plot::set_output_path(pugi::xml_node plot_node) pixels_[2] = pxls[2]; } else { fatal_error( - fmt::format(" must be length 3 in voxel plot {}", id_)); + fmt::format(" must be length 3 in voxel plot {}", id())); } } } -void Plot::set_bg_color(pugi::xml_node plot_node) +void PlottableInterface::set_bg_color(pugi::xml_node plot_node) { // Copy plot background color if (check_for_node(plot_node, "background")) { vector bg_rgb = get_node_array(plot_node, "background"); - if (PlotType::voxel == type_) { - if (mpi::master) { - warning(fmt::format("Background color ignored in voxel plot {}", id_)); - } - } if (bg_rgb.size() == 3) { not_found_ = bg_rgb; } else { - fatal_error(fmt::format("Bad background RGB in plot {}", id_)); + fatal_error(fmt::format("Bad background RGB in plot {}", id())); } } } @@ -332,7 +382,7 @@ void Plot::set_basis(pugi::xml_node plot_node) basis_ = PlotBasis::yz; } else { fatal_error( - fmt::format("Unsupported plot basis '{}' in plot {}", pl_basis, id_)); + fmt::format("Unsupported plot basis '{}' in plot {}", pl_basis, id())); } } } @@ -344,7 +394,7 @@ void Plot::set_origin(pugi::xml_node plot_node) if (pl_origin.size() == 3) { origin_ = pl_origin; } else { - fatal_error(fmt::format("Origin must be length 3 in plot {}", id_)); + fatal_error(fmt::format("Origin must be length 3 in plot {}", id())); } } @@ -358,7 +408,7 @@ void Plot::set_width(pugi::xml_node plot_node) width_.y = pl_width[1]; } else { fatal_error( - fmt::format(" must be length 2 in slice plot {}", id_)); + fmt::format(" must be length 2 in slice plot {}", id())); } } else if (PlotType::voxel == type_) { if (pl_width.size() == 3) { @@ -366,25 +416,25 @@ void Plot::set_width(pugi::xml_node plot_node) width_ = pl_width; } else { fatal_error( - fmt::format(" must be length 3 in voxel plot {}", id_)); + fmt::format(" must be length 3 in voxel plot {}", id())); } } } -void Plot::set_universe(pugi::xml_node plot_node) +void PlottableInterface::set_universe(pugi::xml_node plot_node) { // Copy plot universe level if (check_for_node(plot_node, "level")) { level_ = std::stoi(get_node_value(plot_node, "level")); if (level_ < 0) { - fatal_error(fmt::format("Bad universe level in plot {}", id_)); + fatal_error(fmt::format("Bad universe level in plot {}", id())); } } else { level_ = PLOT_LEVEL_LOWEST; } } -void Plot::set_default_colors(pugi::xml_node plot_node) +void PlottableInterface::set_default_colors(pugi::xml_node plot_node) { // Copy plot color type and initialize all colors randomly std::string pl_color_by = "cell"; @@ -399,7 +449,7 @@ void Plot::set_default_colors(pugi::xml_node plot_node) colors_.resize(model::materials.size()); } else { fatal_error(fmt::format( - "Unsupported plot color type '{}' in plot {}", pl_color_by, id_)); + "Unsupported plot color type '{}' in plot {}", pl_color_by, id())); } for (auto& c : colors_) { @@ -411,28 +461,21 @@ void Plot::set_default_colors(pugi::xml_node plot_node) } } -void Plot::set_user_colors(pugi::xml_node plot_node) +void PlottableInterface::set_user_colors(pugi::xml_node plot_node) { - if (!plot_node.select_nodes("color").empty() && PlotType::voxel == type_) { - if (mpi::master) { - warning( - fmt::format("Color specifications ignored in voxel plot {}", id_)); - } - } - for (auto cn : plot_node.children("color")) { // Make sure 3 values are specified for RGB vector user_rgb = get_node_array(cn, "rgb"); if (user_rgb.size() != 3) { - fatal_error(fmt::format("Bad RGB in plot {}", id_)); + fatal_error(fmt::format("Bad RGB in plot {}", id())); } // Ensure that there is an id for this color specification int col_id; if (check_for_node(cn, "id")) { col_id = std::stoi(get_node_value(cn, "id")); } else { - fatal_error( - fmt::format("Must specify id for color specification in plot {}", id_)); + fatal_error(fmt::format( + "Must specify id for color specification in plot {}", id())); } // Add RGB if (PlotColorBy::cells == color_by_) { @@ -441,7 +484,7 @@ void Plot::set_user_colors(pugi::xml_node plot_node) colors_[col_id] = user_rgb; } else { warning(fmt::format( - "Could not find cell {} specified in plot {}", col_id, id_)); + "Could not find cell {} specified in plot {}", col_id, id())); } } else if (PlotColorBy::mats == color_by_) { if (model::material_map.find(col_id) != model::material_map.end()) { @@ -449,7 +492,7 @@ void Plot::set_user_colors(pugi::xml_node plot_node) colors_[col_id] = user_rgb; } else { warning(fmt::format( - "Could not find material {} specified in plot {}", col_id, id_)); + "Could not find material {} specified in plot {}", col_id, id())); } } } // color node loop @@ -462,7 +505,7 @@ void Plot::set_meshlines(pugi::xml_node plot_node) if (!mesh_line_nodes.empty()) { if (PlotType::voxel == type_) { - warning(fmt::format("Meshlines ignored in voxel plot {}", id_)); + warning(fmt::format("Meshlines ignored in voxel plot {}", id())); } if (mesh_line_nodes.size() == 1) { @@ -476,7 +519,7 @@ void Plot::set_meshlines(pugi::xml_node plot_node) } else { fatal_error(fmt::format( "Must specify a meshtype for meshlines specification in plot {}", - id_)); + id())); } // Ensure that there is a linewidth for this meshlines specification @@ -487,7 +530,7 @@ void Plot::set_meshlines(pugi::xml_node plot_node) } else { fatal_error(fmt::format( "Must specify a linewidth for meshlines specification in plot {}", - id_)); + id())); } // Check for color @@ -496,7 +539,7 @@ void Plot::set_meshlines(pugi::xml_node plot_node) vector ml_rgb = get_node_array(meshlines_node, "color"); if (ml_rgb.size() != 3) { fatal_error( - fmt::format("Bad RGB for meshlines color in plot {}", id_)); + fmt::format("Bad RGB for meshlines color in plot {}", id())); } meshlines_color_ = ml_rgb; } @@ -504,7 +547,8 @@ void Plot::set_meshlines(pugi::xml_node plot_node) // Set mesh based on type if ("ufs" == meshtype) { if (!simulation::ufs_mesh) { - fatal_error(fmt::format("No UFS mesh for meshlines on plot {}", id_)); + fatal_error( + fmt::format("No UFS mesh for meshlines on plot {}", id())); } else { for (int i = 0; i < model::meshes.size(); ++i) { if (const auto* m = @@ -520,7 +564,7 @@ void Plot::set_meshlines(pugi::xml_node plot_node) } else if ("entropy" == meshtype) { if (!simulation::entropy_mesh) { fatal_error( - fmt::format("No entropy mesh for meshlines on plot {}", id_)); + fmt::format("No entropy mesh for meshlines on plot {}", id())); } else { for (int i = 0; i < model::meshes.size(); ++i) { if (const auto* m = @@ -542,7 +586,7 @@ void Plot::set_meshlines(pugi::xml_node plot_node) std::stringstream err_msg; fatal_error(fmt::format("Must specify a mesh id for meshlines tally " "mesh specification in plot {}", - id_)); + id())); } // find the tally index int idx; @@ -550,30 +594,24 @@ void Plot::set_meshlines(pugi::xml_node plot_node) if (err != 0) { fatal_error(fmt::format("Could not find mesh {} specified in " "meshlines for plot {}", - tally_mesh_id, id_)); + tally_mesh_id, id())); } index_meshlines_mesh_ = idx; } else { - fatal_error(fmt::format("Invalid type for meshlines on plot {}", id_)); + fatal_error(fmt::format("Invalid type for meshlines on plot {}", id())); } } else { - fatal_error(fmt::format("Mutliple meshlines specified in plot {}", id_)); + fatal_error(fmt::format("Mutliple meshlines specified in plot {}", id())); } } } -void Plot::set_mask(pugi::xml_node plot_node) +void PlottableInterface::set_mask(pugi::xml_node plot_node) { // Deal with masks pugi::xpath_node_set mask_nodes = plot_node.select_nodes("mask"); if (!mask_nodes.empty()) { - if (PlotType::voxel == type_) { - if (mpi::master) { - warning(fmt::format("Mask ignored in voxel plot {}", id_)); - } - } - if (mask_nodes.size() == 1) { // Get pointer to mask pugi::xml_node mask_node = mask_nodes[0].node(); @@ -582,7 +620,7 @@ void Plot::set_mask(pugi::xml_node plot_node) vector iarray = get_node_array(mask_node, "components"); if (iarray.size() == 0) { fatal_error( - fmt::format("Missing in mask of plot {}", id_)); + fmt::format("Missing in mask of plot {}", id())); } // First we need to change the user-specified identifiers to indices @@ -594,7 +632,7 @@ void Plot::set_mask(pugi::xml_node plot_node) } else { fatal_error(fmt::format("Could not find cell {} specified in the " "mask in plot {}", - col_id, id_)); + col_id, id())); } } else if (PlotColorBy::mats == color_by_) { if (model::material_map.find(col_id) != model::material_map.end()) { @@ -602,7 +640,7 @@ void Plot::set_mask(pugi::xml_node plot_node) } else { fatal_error(fmt::format("Could not find material {} specified in " "the mask in plot {}", - col_id, id_)); + col_id, id())); } } } @@ -620,12 +658,12 @@ void Plot::set_mask(pugi::xml_node plot_node) } } else { - fatal_error(fmt::format("Mutliple masks specified in plot {}", id_)); + fatal_error(fmt::format("Mutliple masks specified in plot {}", id())); } } } -void Plot::set_overlap_color(pugi::xml_node plot_node) +void PlottableInterface::set_overlap_color(pugi::xml_node plot_node) { color_overlaps_ = false; if (check_for_node(plot_node, "show_overlaps")) { @@ -635,13 +673,13 @@ void Plot::set_overlap_color(pugi::xml_node plot_node) if (!color_overlaps_) { warning(fmt::format( "Overlap color specified in plot {} but overlaps won't be shown.", - id_)); + id())); } vector olap_clr = get_node_array(plot_node, "overlap_color"); if (olap_clr.size() == 3) { overlap_color_ = olap_clr; } else { - fatal_error(fmt::format("Bad overlap RGB in plot {}", id_)); + fatal_error(fmt::format("Bad overlap RGB in plot {}", id())); } } } @@ -654,32 +692,37 @@ void Plot::set_overlap_color(pugi::xml_node plot_node) } } -Plot::Plot(pugi::xml_node plot_node) - : index_meshlines_mesh_ {-1}, overlap_color_ {RED} +PlottableInterface::PlottableInterface(pugi::xml_node plot_node) { set_id(plot_node); - set_type(plot_node); - set_output_path(plot_node); set_bg_color(plot_node); - set_basis(plot_node); - set_origin(plot_node); - set_width(plot_node); set_universe(plot_node); set_default_colors(plot_node); set_user_colors(plot_node); - set_meshlines(plot_node); set_mask(plot_node); set_overlap_color(plot_node); -} // End Plot constructor +} + +Plot::Plot(pugi::xml_node plot_node, PlotType type) + : PlottableInterface(plot_node), index_meshlines_mesh_ {-1}, type_(type) +{ + set_output_path(plot_node); + set_basis(plot_node); + set_origin(plot_node); + set_width(plot_node); + set_meshlines(plot_node); + slice_level_ = level_; // Copy level employed in SlicePlotBase::get_map + slice_color_overlaps_ = color_overlaps_; +} //============================================================================== // OUTPUT_PPM writes out a previously generated image to a PPM file //============================================================================== -void output_ppm(Plot const& pl, const ImageData& data) +void output_ppm(const std::string& filename, const ImageData& data) { // Open PPM file for writing - std::string fname = pl.path_plot_; + std::string fname = filename; fname = strtrim(fname); std::ofstream of; @@ -687,14 +730,14 @@ void output_ppm(Plot const& pl, const ImageData& data) // Write header of << "P6\n"; - of << pl.pixels_[0] << " " << pl.pixels_[1] << "\n"; + of << data.shape()[0] << " " << data.shape()[1] << "\n"; of << "255\n"; of.close(); of.open(fname, std::ios::binary | std::ios::app); // Write color for each pixel - for (int y = 0; y < pl.pixels_[1]; y++) { - for (int x = 0; x < pl.pixels_[0]; x++) { + for (int y = 0; y < data.shape()[1]; y++) { + for (int x = 0; x < data.shape()[0]; x++) { RGBColor rgb = data(x, y); of << rgb.red << rgb.green << rgb.blue; } @@ -707,10 +750,10 @@ void output_ppm(Plot const& pl, const ImageData& data) //============================================================================== #ifdef USE_LIBPNG -void output_png(Plot const& pl, const ImageData& data) +void output_png(const std::string& filename, const ImageData& data) { // Open PNG file for writing - std::string fname = pl.path_plot_; + std::string fname = filename; fname = strtrim(fname); auto fp = std::fopen(fname.c_str(), "wb"); @@ -726,8 +769,8 @@ void output_png(Plot const& pl, const ImageData& data) png_init_io(png_ptr, fp); // Write header (8 bit colour depth) - int width = pl.pixels_[0]; - int height = pl.pixels_[1]; + int width = data.shape()[0]; + int height = data.shape()[1]; png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE); png_write_info(png_ptr, info_ptr); @@ -760,13 +803,13 @@ void output_png(Plot const& pl, const ImageData& data) // DRAW_MESH_LINES draws mesh line boundaries on an image //============================================================================== -void draw_mesh_lines(Plot const& pl, ImageData& data) +void Plot::draw_mesh_lines(ImageData& data) const { RGBColor rgb; - rgb = pl.meshlines_color_; + rgb = meshlines_color_; int ax1, ax2; - switch (pl.basis_) { + switch (basis_) { case PlotBasis::xy: ax1 = 0; ax2 = 1; @@ -783,45 +826,45 @@ void draw_mesh_lines(Plot const& pl, ImageData& data) UNREACHABLE(); } - Position ll_plot {pl.origin_}; - Position ur_plot {pl.origin_}; + Position ll_plot {origin_}; + Position ur_plot {origin_}; - 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.; + ll_plot[ax1] -= width_[0] / 2.; + ll_plot[ax2] -= width_[1] / 2.; + ur_plot[ax1] += width_[0] / 2.; + ur_plot[ax2] += width_[1] / 2.; Position width = ur_plot - ll_plot; // Find the (axis-aligned) lines of the mesh that intersect this plot. auto axis_lines = - model::meshes[pl.index_meshlines_mesh_]->plot(ll_plot, ur_plot); + model::meshes[index_meshlines_mesh_]->plot(ll_plot, ur_plot); // 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.back() - ll_plot[ax2]) / width[ax2]; - ax2_min = (1.0 - frac) * pl.pixels_[1]; + ax2_min = (1.0 - frac) * pixels_[1]; if (ax2_min < 0) ax2_min = 0; frac = (axis_lines.second.front() - ll_plot[ax2]) / width[ax2]; - ax2_max = (1.0 - frac) * pl.pixels_[1]; - if (ax2_max > pl.pixels_[1]) - ax2_max = pl.pixels_[1]; + ax2_max = (1.0 - frac) * pixels_[1]; + if (ax2_max > pixels_[1]) + ax2_max = pixels_[1]; } else { ax2_min = 0; - ax2_max = pl.pixels_[1]; + ax2_max = pixels_[1]; } // 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]; + int ax1_ind = frac * 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]) + for (int plus = 0; plus <= meshlines_width_; plus++) { + if (ax1_ind + plus >= 0 && ax1_ind + plus < pixels_[0]) data(ax1_ind + plus, ax2_ind) = rgb; - if (ax1_ind - plus >= 0 && ax1_ind - plus < pl.pixels_[0]) + if (ax1_ind - plus >= 0 && ax1_ind - plus < pixels_[0]) data(ax1_ind - plus, ax2_ind) = rgb; } } @@ -831,60 +874,57 @@ void draw_mesh_lines(Plot const& pl, ImageData& data) 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]; + ax1_min = frac * 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]; + ax1_max = frac * pixels_[0]; + if (ax1_max > pixels_[0]) + ax1_max = pixels_[0]; } else { ax1_min = 0; - ax1_max = pl.pixels_[0]; + ax1_max = 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 = (1.0 - frac) * pl.pixels_[1]; + int ax2_ind = (1.0 - frac) * 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]) + for (int plus = 0; plus <= meshlines_width_; plus++) { + if (ax2_ind + plus >= 0 && ax2_ind + plus < pixels_[1]) data(ax1_ind, ax2_ind + plus) = rgb; - if (ax2_ind - plus >= 0 && ax2_ind - plus < pl.pixels_[1]) + if (ax2_ind - plus >= 0 && ax2_ind - plus < pixels_[1]) data(ax1_ind, ax2_ind - plus) = rgb; } } } } -//============================================================================== -// CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D -// geometry visualization. It works the same way as create_image by dragging a -// particle across the geometry for the specified number of voxels. The first 3 -// int's in the binary are the number of x, y, and z voxels. The next 3 -// double's are the widths of the voxels in the x, y, and z directions. The -// next 3 double's are the x, y, and z coordinates of the lower left -// point. Finally the binary is filled with entries of four int's each. Each -// 'row' in the binary contains four int's: 3 for x,y,z position and 1 for -// cell or material id. For 1 million voxels this produces a file of -// approximately 15MB. -// ============================================================================= - -void create_voxel(Plot const& pl) +/* outputs a binary file that can be input into silomesh for 3D geometry + * visualization. It works the same way as create_image by dragging a particle + * across the geometry for the specified number of voxels. The first 3 int's in + * the binary are the number of x, y, and z voxels. The next 3 double's are + * the widths of the voxels in the x, y, and z directions. The next 3 double's + * are the x, y, and z coordinates of the lower left point. Finally the binary + * is filled with entries of four int's each. Each 'row' in the binary contains + * four int's: 3 for x,y,z position and 1 for cell or material id. For 1 + * million voxels this produces a file of approximately 15MB. + */ +void Plot::create_voxel() const { // compute voxel widths in each direction array vox; - vox[0] = pl.width_[0] / (double)pl.pixels_[0]; - vox[1] = pl.width_[1] / (double)pl.pixels_[1]; - vox[2] = pl.width_[2] / (double)pl.pixels_[2]; + vox[0] = width_[0] / static_cast(pixels_[0]); + vox[1] = width_[1] / static_cast(pixels_[1]); + vox[2] = width_[2] / static_cast(pixels_[2]); // initial particle position - Position ll = pl.origin_ - pl.width_ / 2.; + Position ll = origin_ - width_ / 2.; // Open binary plot file for writing std::ofstream of; - std::string fname = std::string(pl.path_plot_); + std::string fname = std::string(path_plot_); fname = strtrim(fname); hid_t file_id = file_open(fname, 'w'); @@ -900,7 +940,7 @@ void create_voxel(Plot const& pl) // Write current date and time write_attribute(file_id, "date_and_time", time_stamp().c_str()); array pixels; - std::copy(pl.pixels_.begin(), pl.pixels_.end(), pixels.begin()); + std::copy(pixels_.begin(), pixels_.end(), pixels.begin()); write_attribute(file_id, "num_voxels", pixels); write_attribute(file_id, "voxel_width", vox); write_attribute(file_id, "lower_left", ll); @@ -908,24 +948,24 @@ void create_voxel(Plot const& pl) // Create dataset for voxel data -- note that the dimensions are reversed // since we want the order in the file to be z, y, x hsize_t dims[3]; - dims[0] = pl.pixels_[2]; - dims[1] = pl.pixels_[1]; - dims[2] = pl.pixels_[0]; + dims[0] = pixels_[2]; + dims[1] = pixels_[1]; + dims[2] = pixels_[0]; hid_t dspace, dset, memspace; voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace); - PlotBase pltbase; - pltbase.width_ = pl.width_; - pltbase.origin_ = pl.origin_; + SlicePlotBase pltbase; + pltbase.width_ = width_; + pltbase.origin_ = origin_; pltbase.basis_ = PlotBasis::xy; - pltbase.pixels_ = pl.pixels_; - pltbase.level_ = -1; // all universes for voxel files - pltbase.color_overlaps_ = pl.color_overlaps_; + pltbase.pixels_ = pixels_; + pltbase.slice_color_overlaps_ = color_overlaps_; ProgressBar pb; - for (int z = 0; z < pl.pixels_[2]; z++) { + for (int z = 0; z < pixels_[2]; z++) { // update progress bar - pb.set_value(100. * (double)z / (double)(pl.pixels_[2] - 1)); + pb.set_value( + 100. * static_cast(z) / static_cast((pixels_[2] - 1))); // update z coordinate pltbase.origin_.z = ll.z + z * vox[2]; @@ -934,7 +974,7 @@ void create_voxel(Plot const& pl) IdData ids = pltbase.get_map(); // select only cell/material ID data and flip the y-axis - int idx = pl.color_by_ == PlotColorBy::cells ? 0 : 2; + int idx = color_by_ == PlotColorBy::cells ? 0 : 2; xt::xtensor data_slice = xt::view(ids.data_, xt::all(), xt::all(), idx); xt::xtensor data_flipped = xt::flip(data_slice, 0); @@ -986,16 +1026,568 @@ RGBColor random_color(void) int(prn(&model::plotter_seed) * 255), int(prn(&model::plotter_seed) * 255)}; } +ProjectionPlot::ProjectionPlot(pugi::xml_node node) : PlottableInterface(node) +{ + set_output_path(node); + set_look_at(node); + set_camera_position(node); + set_field_of_view(node); + set_pixels(node); + set_opacities(node); + set_orthographic_width(node); + set_wireframe_thickness(node); + set_wireframe_ids(node); + set_wireframe_color(node); + + if (check_for_node(node, "orthographic_width") && + check_for_node(node, "field_of_view")) + fatal_error("orthographic_width and field_of_view are mutually exclusive " + "parameters."); +} + +void ProjectionPlot::set_wireframe_color(pugi::xml_node plot_node) +{ + // Copy plot background color + if (check_for_node(plot_node, "wireframe_color")) { + vector w_rgb = get_node_array(plot_node, "wireframe_color"); + if (w_rgb.size() == 3) { + wireframe_color_ = w_rgb; + } else { + fatal_error(fmt::format("Bad wireframe RGB in plot {}", id())); + } + } +} + +void ProjectionPlot::set_output_path(pugi::xml_node node) +{ + // Set output file path + std::string filename; + + if (check_for_node(node, "filename")) { + filename = get_node_value(node, "filename"); + } else { + filename = fmt::format("plot_{}", id()); + } + +#ifdef USE_LIBPNG + if (!file_extension_present(filename, "png")) + filename.append(".png"); +#else + if (!file_extension_present(filename, "ppm")) + filename.append(".ppm"); +#endif + path_plot_ = filename; +} + +// Advances to the next boundary from outside the geometry +// Returns -1 if no intersection found, and the surface index +// if an intersection was found. +int ProjectionPlot::advance_to_boundary_from_void(Particle& p) +{ + constexpr double scoot = 1e-5; + double min_dist = {INFINITY}; + auto coord = p.coord(0); + Universe* uni = model::universes[model::root_universe].get(); + int intersected_surface = -1; + for (auto c_i : uni->cells_) { + auto dist = model::cells.at(c_i)->distance(coord.r, coord.u, 0, &p); + if (dist.first < min_dist) { + min_dist = dist.first; + intersected_surface = dist.second; + } + } + if (min_dist > 1e300) + return -1; + else { // advance the particle + for (int j = 0; j < p.n_coord(); ++j) + p.coord(j).r += (min_dist + scoot) * p.coord(j).u; + return std::abs(intersected_surface); + } +} + +bool ProjectionPlot::trackstack_equivalent( + const std::vector& track1, + const std::vector& track2) const +{ + if (wireframe_ids_.empty()) { + // Draw wireframe for all surfaces/cells/materials + if (track1.size() != track2.size()) + return false; + for (int i = 0; i < track1.size(); ++i) { + if (track1[i].id != track2[i].id || + track1[i].surface != track2[i].surface) { + return false; + } + } + return true; + } else { + // This runs in O(nm) where n is the intersection stack size + // and m is the number of IDs we are wireframing. A simpler + // algorithm can likely be found. + for (const int id : wireframe_ids_) { + int t1_i = 0; + int t2_i = 0; + + // Advance to first instance of the ID + while (t1_i < track1.size() && t2_i < track2.size()) { + while (t1_i < track1.size() && track1[t1_i].id != id) + t1_i++; + while (t2_i < track2.size() && track2[t2_i].id != id) + t2_i++; + + // This one is really important! + if ((t1_i == track1.size() && t2_i != track2.size()) || + (t1_i != track1.size() && t2_i == track2.size())) + return false; + if (t1_i == track1.size() && t2_i == track2.size()) + break; + // Check if surface different + if (track1[t1_i].surface != track2[t2_i].surface) + return false; + + // Pretty sure this should not be used: + // if (t2_i != track2.size() - 1 && + // t1_i != track1.size() - 1 && + // track1[t1_i+1].id != track2[t2_i+1].id) return false; + if (t2_i != 0 && t1_i != 0 && + track1[t1_i - 1].surface != track2[t2_i - 1].surface) + return false; + + // Check if neighboring cells are different + // if (track1[t1_i ? t1_i - 1 : 0].id != track2[t2_i ? t2_i - 1 : 0].id) + // return false; if (track1[t1_i < track1.size() - 1 ? t1_i + 1 : t1_i + // ].id != + // track2[t2_i < track2.size() - 1 ? t2_i + 1 : t2_i].id) return + // false; + t1_i++, t2_i++; + } + } + return true; + } +} + +void ProjectionPlot::create_output() const +{ + // Get centerline vector for camera-to-model. We create vectors around this + // that form a pixel array, and then trace rays along that. + auto up = up_ / up_.norm(); + Direction looking_direction = look_at_ - camera_position_; + looking_direction /= looking_direction.norm(); + if (std::abs(std::abs(looking_direction.dot(up)) - 1.0) < 1e-9) + fatal_error("Up vector cannot align with vector between camera position " + "and look_at!"); + Direction cam_yaxis = looking_direction.cross(up); + cam_yaxis /= cam_yaxis.norm(); + Direction cam_zaxis = cam_yaxis.cross(looking_direction); + cam_zaxis /= cam_zaxis.norm(); + + // Transformation matrix for directions + std::vector camera_to_model = {looking_direction.x, cam_yaxis.x, + cam_zaxis.x, looking_direction.y, cam_yaxis.y, cam_zaxis.y, + looking_direction.z, cam_yaxis.z, cam_zaxis.z}; + + // Now we convert to the polar coordinate system with the polar angle + // measuring the angle from the vector up_. Phi is the rotation about up_. For + // now, up_ is hard-coded to be +z. + constexpr double DEGREE_TO_RADIAN = M_PI / 180.0; + double horiz_fov_radians = horizontal_field_of_view_ * DEGREE_TO_RADIAN; + double p0 = static_cast(pixels_[0]); + double p1 = static_cast(pixels_[1]); + double vert_fov_radians = horiz_fov_radians * p1 / p0; + double dphi = horiz_fov_radians / p0; + double dmu = vert_fov_radians / p1; + + size_t width = pixels_[0]; + size_t height = pixels_[1]; + ImageData data({width, height}, not_found_); + + // This array marks where the initial wireframe was drawn. + // We convolve it with a filter that gets adjusted with the + // wireframe thickness in order to thicken the lines. + xt::xtensor wireframe_initial({width, height}, 0); + + /* Holds all of the track segments for the current rendered line of pixels. + * old_segments holds a copy of this_line_segments from the previous line. + * By holding both we can check if the cell/material intersection stack + * differs from the left or upper neighbor. This allows a robustly drawn + * wireframe. If only checking the left pixel (which requires substantially + * less memory), the wireframe tends to be spotty and be disconnected for + * surface edges oriented horizontally in the rendering. + * + * Note that a vector of vectors is required rather than a 2-tensor, + * since the stack size varies within each column. + */ +#ifdef _OPENMP + const int n_threads = omp_get_max_threads(); +#else + const int n_threads = 1; +#endif + std::vector>> this_line_segments( + n_threads); + for (int t = 0; t < n_threads; ++t) { + this_line_segments[t].resize(pixels_[0]); + } + + // The last thread writes to this, and the first thread reads from it. + std::vector> old_segments(pixels_[0]); + +#pragma omp parallel + { + +#ifdef _OPENMP + const int n_threads = omp_get_max_threads(); + const int tid = omp_get_thread_num(); +#else + int n_threads = 1; + int tid = 0; +#endif + + SourceSite s; // Where particle starts from (camera) + s.E = 1; + s.wgt = 1; + s.delayed_group = 0; + s.particle = ParticleType::photon; // just has to be something reasonable + s.parent_id = 1; + s.progeny_id = 2; + s.r = camera_position_; + + Particle p; + s.u.x = 1.0; + s.u.y = 0.0; + s.u.z = 0.0; + p.from_source(&s); + + int vert = tid; + for (int iter = 0; iter <= pixels_[1] / n_threads; iter++) { + + // Save bottom line of current work chunk to compare against later + // I used to have this inside the below if block, but it causes a + // spurious line to be drawn at the bottom of the image. Not sure + // why, but moving it here fixes things. + if (tid == n_threads - 1) + old_segments = this_line_segments[n_threads - 1]; + + if (vert < pixels_[1]) { + + for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + + // Generate the starting position/direction of the ray + if (orthographic_width_ == 0.0) { // perspective projection + double this_phi = + -horiz_fov_radians / 2.0 + dphi * horiz + 0.5 * dphi; + double this_mu = + -vert_fov_radians / 2.0 + dmu * vert + M_PI / 2.0 + 0.5 * dmu; + Direction camera_local_vec; + camera_local_vec.x = std::cos(this_phi) * std::sin(this_mu); + camera_local_vec.y = std::sin(this_phi) * std::sin(this_mu); + camera_local_vec.z = std::cos(this_mu); + s.u = camera_local_vec.rotate(camera_to_model); + } else { // orthographic projection + s.u = looking_direction; + + double x_pix_coord = (static_cast(horiz) - p0 / 2.0) / p0; + double y_pix_coord = (static_cast(vert) - p1 / 2.0) / p0; + s.r = camera_position_ + + cam_yaxis * x_pix_coord * orthographic_width_ + + cam_zaxis * y_pix_coord * orthographic_width_; + } + + p.from_source(&s); // put particle at camera + bool hitsomething = false; + bool intersection_found = true; + int loop_counter = 0; + + this_line_segments[tid][horiz].clear(); + + int first_surface = + -1; // surface first passed when entering the model + bool first_inside_model = true; // false after entering the model + while (intersection_found) { + bool inside_cell = false; + + int32_t i_surface = std::abs(p.surface()) - 1; + if (i_surface > 0 && model::surfaces[i_surface]->geom_type_ == GeometryType::DAG) { +#ifdef DAGMC + int32_t i_cell = next_cell(i_surface, p.cell_last(p.n_coord() - 1), p.lowest_coord().universe); + inside_cell = i_cell >= 0; +#else + fatal_error("Not compiled for DAGMC, but somehow you have a DAGCell!"); +#endif + } else { + inside_cell = exhaustive_find_cell(p); + } + + if (inside_cell) { + + // This allows drawing wireframes with surface intersection + // edges on the model boundary for the same cell. + if (first_inside_model) { + this_line_segments[tid][horiz].emplace_back( + color_by_ == PlotColorBy::mats + ? p.material() + : p.coord(p.n_coord() - 1).cell, + 0.0, first_surface); + first_inside_model = false; + } + + hitsomething = true; + intersection_found = true; + auto dist = distance_to_boundary(p); + this_line_segments[tid][horiz].emplace_back( + color_by_ == PlotColorBy::mats ? p.material() + : p.coord(p.n_coord() - 1).cell, + dist.distance, std::abs(dist.surface_index)); + + // Advance particle + for (int lev = 0; lev < p.n_coord(); ++lev) { + p.coord(lev).r += dist.distance * p.coord(lev).u; + } + p.surface() = dist.surface_index; + p.n_coord_last() = p.n_coord(); + p.n_coord() = dist.coord_level; + if (dist.lattice_translation[0] != 0 || + dist.lattice_translation[1] != 0 || + dist.lattice_translation[2] != 0) { + cross_lattice(p, dist); + } + + } else { + first_surface = advance_to_boundary_from_void(p); + intersection_found = + first_surface != -1; // -1 if no surface found + } + loop_counter++; + if (loop_counter > MAX_INTERSECTIONS) + fatal_error("Infinite loop in projection plot"); + } + + // Now color the pixel based on what we have intersected... + // Loops backwards over intersections. + Position current_color( + not_found_.red, not_found_.green, not_found_.blue); + const auto& segments = this_line_segments[tid][horiz]; + for (unsigned i = segments.size(); i-- > 0;) { + int colormap_idx = segments[i].id; + RGBColor seg_color = colors_[colormap_idx]; + Position seg_color_vec( + seg_color.red, seg_color.green, seg_color.blue); + double mixing = std::exp(-xs_[colormap_idx] * segments[i].length); + current_color = + current_color * mixing + (1.0 - mixing) * seg_color_vec; + RGBColor result; + result.red = static_cast(current_color.x); + result.green = static_cast(current_color.y); + result.blue = static_cast(current_color.z); + data(horiz, vert) = result; + } + + // Check to draw wireframe in horizontal direction. No inter-thread + // comm. + if (horiz > 0) { + if (!trackstack_equivalent(this_line_segments[tid][horiz], + this_line_segments[tid][horiz - 1])) { + wireframe_initial(horiz, vert) = 1; + } + } + } + } // end "if" vert in correct range + + // We require a barrier before comparing vertical neighbors' intersection + // stacks. i.e. all threads must be done with their line. +#pragma omp barrier + + // Now that the horizontal line has finished rendering, we can fill in + // wireframe entries that require comparison among all the threads. Hence + // the omp barrier being used. It has to be OUTSIDE any if blocks! + if (vert < pixels_[1]) { + // Loop over horizontal pixels, checking intersection stack of upper + // neighbor + + const std::vector>* top_cmp = nullptr; + if (tid == 0) + top_cmp = &old_segments; + else + top_cmp = &this_line_segments[tid - 1]; + + for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + if (!trackstack_equivalent( + this_line_segments[tid][horiz], (*top_cmp)[horiz])) { + wireframe_initial(horiz, vert) = 1; + } + } + } + + // We need another barrier to ensure threads don't proceed to modify their + // intersection stacks on that horizontal line while others are + // potentially still working on the above. +#pragma omp barrier + vert += n_threads; + } + } // end omp parallel + + // Now thicken the wireframe lines and apply them to our image + for (int vert = 0; vert < pixels_[1]; ++vert) { + for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + if (wireframe_initial(horiz, vert)) { + if (wireframe_thickness_ == 1) + data(horiz, vert) = wireframe_color_; + for (int i = -wireframe_thickness_ / 2; i < wireframe_thickness_ / 2; + ++i) + 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 + int w_i = std::max(std::min(horiz + i, pixels_[0] - 1), 0); + int w_j = std::max(std::min(vert + j, pixels_[1] - 1), 0); + data(w_i, w_j) = wireframe_color_; + } + } + } + } + +#ifdef USE_LIBPNG + output_png(path_plot(), data); +#else + output_ppm(path_plot(), data); +#endif +} + +void ProjectionPlot::print_info() const +{ + fmt::print("Plot Type: Projection\n"); + fmt::print("Camera position: {} {} {}\n", camera_position_.x, + camera_position_.y, camera_position_.z); + fmt::print("Look at: {} {} {}\n", look_at_.x, look_at_.y, look_at_.z); + fmt::print( + "Horizontal field of view: {} degrees\n", horizontal_field_of_view_); + fmt::print("Pixels: {} {}\n", pixels_[0], pixels_[1]); +} + +void ProjectionPlot::set_opacities(pugi::xml_node node) +{ + xs_.resize(colors_.size(), 1e6); // set to large value for opaque by default + + for (auto cn : node.children("color")) { + // Make sure 3 values are specified for RGB + double user_xs = std::stod(get_node_value(cn, "xs")); + int col_id = std::stoi(get_node_value(cn, "id")); + + // Add RGB + if (PlotColorBy::cells == color_by_) { + if (model::cell_map.find(col_id) != model::cell_map.end()) { + col_id = model::cell_map[col_id]; + xs_[col_id] = user_xs; + } else { + warning(fmt::format( + "Could not find cell {} specified in plot {}", col_id, id())); + } + } else if (PlotColorBy::mats == color_by_) { + if (model::material_map.find(col_id) != model::material_map.end()) { + col_id = model::material_map[col_id]; + xs_[col_id] = user_xs; + } else { + warning(fmt::format( + "Could not find material {} specified in plot {}", col_id, id())); + } + } + } +} + +void ProjectionPlot::set_orthographic_width(pugi::xml_node node) +{ + if (check_for_node(node, "orthographic_width")) { + double orthographic_width = + std::stod(get_node_value(node, "orthographic_width", true)); + if (orthographic_width < 0.0) + fatal_error("Requires positive orthographic_width"); + orthographic_width_ = orthographic_width; + } +} + +void ProjectionPlot::set_wireframe_thickness(pugi::xml_node node) +{ + if (check_for_node(node, "wireframe_thickness")) { + int wireframe_thickness = + std::stoi(get_node_value(node, "wireframe_thickness", true)); + if (wireframe_thickness < 0) + fatal_error("Requires non-negative wireframe thickness"); + wireframe_thickness_ = wireframe_thickness; + } +} + +void ProjectionPlot::set_wireframe_ids(pugi::xml_node node) +{ + if (check_for_node(node, "wireframe_ids")) { + wireframe_ids_ = get_node_array(node, "wireframe_ids"); + // It is read in as actual ID values, but we have to convert to indices in + // mat/cell array + for (auto& x : wireframe_ids_) + x = color_by_ == PlotColorBy::mats ? model::material_map[x] + : model::cell_map[x]; + } + // We make sure the list is sorted in order to later use + // std::binary_search. + std::sort(wireframe_ids_.begin(), wireframe_ids_.end()); +} + +void ProjectionPlot::set_pixels(pugi::xml_node node) +{ + vector pxls = get_node_array(node, "pixels"); + if (pxls.size() != 2) + fatal_error( + fmt::format(" must be length 2 in projection plot {}", id())); + pixels_[0] = pxls[0]; + pixels_[1] = pxls[1]; +} + +void ProjectionPlot::set_camera_position(pugi::xml_node node) +{ + vector camera_pos = get_node_array(node, "camera_position"); + if (camera_pos.size() != 3) { + fatal_error( + fmt::format("look_at element must have three floating point values")); + } + camera_position_.x = camera_pos[0]; + camera_position_.y = camera_pos[1]; + camera_position_.z = camera_pos[2]; +} + +void ProjectionPlot::set_look_at(pugi::xml_node node) +{ + vector look_at = get_node_array(node, "look_at"); + if (look_at.size() != 3) { + fatal_error("look_at element must have three floating point values"); + } + look_at_.x = look_at[0]; + look_at_.y = look_at[1]; + look_at_.z = look_at[2]; +} + +void ProjectionPlot::set_field_of_view(pugi::xml_node node) +{ + // Defaults to 70 degree horizontal field of view (see .h file) + if (check_for_node(node, "field_of_view")) { + double fov = std::stod(get_node_value(node, "field_of_view", true)); + if (fov < 180.0 && fov > 0.0) { + horizontal_field_of_view_ = fov; + } else { + fatal_error(fmt::format( + "Field of view for plot {} out-of-range. Must be in (0, 180).", id())); + } + } +} + extern "C" int openmc_id_map(const void* plot, int32_t* data_out) { - auto plt = reinterpret_cast(plot); + auto plt = reinterpret_cast(plot); if (!plt) { set_errmsg("Invalid slice pointer passed to openmc_id_map"); return OPENMC_E_INVALID_ARGUMENT; } - if (plt->color_overlaps_ && model::overlap_check_count.size() == 0) { + if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) { model::overlap_check_count.resize(model::cells.size()); } @@ -1010,13 +1602,13 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out) extern "C" int openmc_property_map(const void* plot, double* data_out) { - auto plt = reinterpret_cast(plot); + auto plt = reinterpret_cast(plot); if (!plt) { set_errmsg("Invalid slice pointer passed to openmc_id_map"); return OPENMC_E_INVALID_ARGUMENT; } - if (plt->color_overlaps_ && model::overlap_check_count.size() == 0) { + if (plt->slice_color_overlaps_ && model::overlap_check_count.size() == 0) { model::overlap_check_count.resize(model::cells.size()); } diff --git a/tests/regression_tests/distribmat/inputs_true.dat b/tests/regression_tests/distribmat/inputs_true.dat index f6c43c19a5..adc9bf3fbe 100644 --- a/tests/regression_tests/distribmat/inputs_true.dat +++ b/tests/regression_tests/distribmat/inputs_true.dat @@ -48,14 +48,14 @@ + 400 400 0 0 0 7 7 - 400 400 + 400 400 0 0 0 7 7 - 400 400 diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index 2fd6a5b795..30d9ab0d08 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -1,60 +1,7 @@ -import glob -import hashlib -import os - -import h5py -import openmc - -from tests.testing_harness import TestHarness +from tests.testing_harness import PlotTestHarness 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 *.png 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('.png'): - # Add PNG 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'].tobytes() - outstr += fh.attrs['lower_left'].tobytes() - outstr += fh.attrs['voxel_width'].tobytes() - outstr += fh['data'][()].tobytes() - - # Hash the information and return. - sha512 = hashlib.sha512() - sha512.update(outstr) - outstr = sha512.hexdigest() - - return outstr - - def test_plot(): harness = PlotTestHarness(('plot_1.png', 'plot_2.png', 'plot_3.png', 'plot_4.h5')) diff --git a/tests/regression_tests/plot_overlaps/test.py b/tests/regression_tests/plot_overlaps/test.py index cf23abc27f..0828c6e255 100644 --- a/tests/regression_tests/plot_overlaps/test.py +++ b/tests/regression_tests/plot_overlaps/test.py @@ -1,60 +1,7 @@ -import glob -import hashlib -import os - -import h5py -import openmc - -from tests.testing_harness import TestHarness +from tests.testing_harness import PlotTestHarness 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 *.png 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('.png'): - # Add PNG 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'].tobytes() - outstr += fh.attrs['lower_left'].tobytes() - outstr += fh.attrs['voxel_width'].tobytes() - outstr += fh['data'][()].tobytes() - - # Hash the information and return. - sha512 = hashlib.sha512() - sha512.update(outstr) - outstr = sha512.hexdigest() - - return outstr - - def test_plot_overlap(): harness = PlotTestHarness(('plot_1.png', 'plot_2.png', 'plot_3.png', 'plot_4.h5')) diff --git a/tests/regression_tests/plot_projections/__init__.py b/tests/regression_tests/plot_projections/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/plot_projections/geometry.xml b/tests/regression_tests/plot_projections/geometry.xml new file mode 100644 index 0000000000..a648dfe92f --- /dev/null +++ b/tests/regression_tests/plot_projections/geometry.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/tests/regression_tests/plot_projections/materials.xml b/tests/regression_tests/plot_projections/materials.xml new file mode 100644 index 0000000000..35f0de3128 --- /dev/null +++ b/tests/regression_tests/plot_projections/materials.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/tests/regression_tests/plot_projections/plots.xml b/tests/regression_tests/plot_projections/plots.xml new file mode 100644 index 0000000000..85689da147 --- /dev/null +++ b/tests/regression_tests/plot_projections/plots.xml @@ -0,0 +1,55 @@ + + + + + 0. 0. 0. + 20. 20. 20. + 200 200 + + + + 70 + + + + 0. 0. 0. + 10. 10. 0. + 25 25 + 200 200 + + 90 + 4 + example1 + + + + 0. 0. 0. + 20. 20. 20. + 200 200 + 240 240 240 + example2.png + 2 + + + + 0. 0. 0. + 0. 10.0 20. + 200 200 + 110 240 240 + example3.png + + + + 0. 0. 0. + 10. 10. 10. + 25 25 + 200 200 + 25.0 + 2 + orthographic_example1 + + + + + + diff --git a/tests/regression_tests/plot_projections/results_true.dat b/tests/regression_tests/plot_projections/results_true.dat new file mode 100644 index 0000000000..a2cb57b1e9 --- /dev/null +++ b/tests/regression_tests/plot_projections/results_true.dat @@ -0,0 +1 @@ +6448093f8acaa2af02452ec1b29461ffaf89628bb002cc8873a3a6831fe5d613645a197033597e5a064896388f19d0782f05e78ce474354ae2a76c35ff4551d7 \ No newline at end of file diff --git a/tests/regression_tests/plot_projections/settings.xml b/tests/regression_tests/plot_projections/settings.xml new file mode 100644 index 0000000000..adf256d2d4 --- /dev/null +++ b/tests/regression_tests/plot_projections/settings.xml @@ -0,0 +1,13 @@ + + + + plot + + + 5 4 3 + -10 -10 -10 + 10 10 10 + + 1 + + diff --git a/tests/regression_tests/plot_projections/tallies.xml b/tests/regression_tests/plot_projections/tallies.xml new file mode 100644 index 0000000000..b7e678ca0f --- /dev/null +++ b/tests/regression_tests/plot_projections/tallies.xml @@ -0,0 +1,20 @@ + + + + + -10 10 + -10 10 + -10 0 5 7.5 8.75 10 + + + + mesh + 2 + + + + 1 + total + + + diff --git a/tests/regression_tests/plot_projections/test.py b/tests/regression_tests/plot_projections/test.py new file mode 100644 index 0000000000..cf18503f4b --- /dev/null +++ b/tests/regression_tests/plot_projections/test.py @@ -0,0 +1,6 @@ +from tests.testing_harness import PlotTestHarness +from tests.regression_tests import config + +def test_plot(): + harness = PlotTestHarness(('plot_1.png', 'example1.png', 'example2.png', 'example3.png', 'orthographic_example1.png')) + harness.main() diff --git a/tests/regression_tests/plot_voxel/test.py b/tests/regression_tests/plot_voxel/test.py index c4fa20c80f..e034decf79 100644 --- a/tests/regression_tests/plot_voxel/test.py +++ b/tests/regression_tests/plot_voxel/test.py @@ -1,62 +1,11 @@ -import glob -import hashlib -import os -from subprocess import check_call - -import h5py -import openmc import pytest -from tests.testing_harness import TestHarness +from tests.testing_harness import PlotTestHarness from tests.regression_tests import config vtk = pytest.importorskip('vtk') -class PlotVoxelTestHarness(TestHarness): - """Specialized TestHarness for running OpenMC voxel plot 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']) - - check_call(['../../../scripts/openmc-voxel-to-vtk'] + - glob.glob('plot_4.h5')) - - def _test_output_created(self): - """Make sure plots have 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('.h5'): - # Add voxel data to results - with h5py.File(fname, 'r') as fh: - outstr += fh.attrs['filetype'] - outstr += fh.attrs['num_voxels'].tobytes() - outstr += fh.attrs['lower_left'].tobytes() - outstr += fh.attrs['voxel_width'].tobytes() - outstr += fh['data'][()].tobytes() - - # Hash the information and return. - sha512 = hashlib.sha512() - sha512.update(outstr) - outstr = sha512.hexdigest() - - return outstr - def test_plot_voxel(): - harness = PlotVoxelTestHarness(('plot_4.h5', 'plot.vti')) + harness = PlotTestHarness(('plot_4.h5', 'plot.vti'), voxel_convert_checks=['plot_4.h5']) harness.main() diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 157f6e6402..f0edf4ea44 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -1,6 +1,8 @@ from difflib import unified_diff +from subprocess import check_call import filecmp import glob +import h5py import hashlib import os import shutil @@ -380,3 +382,55 @@ class HashedPyAPITestHarness(PyAPITestHarness): def _get_results(self): """Digest info in the statepoint and return as a string.""" return super()._get_results(True) + + +class PlotTestHarness(TestHarness): + """Specialized TestHarness for running OpenMC plotting tests.""" + def __init__(self, plot_names, voxel_convert_checks=[]): + super().__init__(None) + self._plot_names = plot_names + self._voxel_convert_checks = voxel_convert_checks + + def _run_openmc(self): + openmc.plot_geometry(openmc_exec=config['exe']) + + # Check that voxel h5 can be converted to vtk + for voxel_h5_filename in self._voxel_convert_checks: + check_call(['../../../scripts/openmc-voxel-to-vtk'] + + glob.glob(voxel_h5_filename)) + + def _test_output_created(self): + """Make sure *.png 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('.png'): + # Add PNG 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'].tobytes() + outstr += fh.attrs['lower_left'].tobytes() + outstr += fh.attrs['voxel_width'].tobytes() + outstr += fh['data'][()].tobytes() + + # Hash the information and return. + sha512 = hashlib.sha512() + sha512.update(outstr) + outstr = sha512.hexdigest() + + return outstr diff --git a/tests/unit_tests/test_plots.py b/tests/unit_tests/test_plots.py index 2db396cd5b..183341c26e 100644 --- a/tests/unit_tests/test_plots.py +++ b/tests/unit_tests/test_plots.py @@ -37,15 +37,47 @@ def myplot(): } return plot +@pytest.fixture(scope='module') +def myprojectionplot(): + plot = openmc.ProjectionPlot(name='myprojectionplot') + plot.look_at = (0.0, 0.0, 0.0) + plot.camera_position = (4.0, 3.0, 0.0) + plot.pixels = (500, 500) + plot.filename = 'myprojectionplot' + plot.background = (0, 0, 0) + plot.background = 'black' + + plot.color_by = 'material' + m1, m2 = openmc.Material(), openmc.Material() + plot.colors = {m1: (0, 255, 0), m2: (0, 0, 255)} + plot.colors = {m1: 'green', m2: 'blue'} + plot.xs = {m1: 1.0, m2: 0.01} + + plot.mask_components = [openmc.Material()] + plot.mask_background = (255, 255, 255) + plot.mask_background = 'white' + + plot.overlap_color = (255, 211, 0) + plot.overlap_color = 'yellow' + + plot.wireframe_thickness = 2 + + plot.level = 1 + return plot def test_attributes(myplot): assert myplot.name == 'myplot' +def test_attributes_proj(myprojectionplot): + assert myprojectionplot.name == 'myprojectionplot' def test_repr(myplot): r = repr(myplot) assert isinstance(r, str) +def test_repr_proj(myprojectionplot): + r = repr(myprojectionplot) + assert isinstance(r, str) def test_from_geometry(): width = 25. @@ -89,6 +121,18 @@ def test_xml_element(myplot): assert getattr(newplot, attr) == getattr(myplot, attr), attr +def test_to_xml_element_proj(myprojectionplot): + elem = myprojectionplot.to_xml_element() + assert 'id' in elem.attrib + assert 'color_by' in elem.attrib + assert 'type' in elem.attrib + assert elem.find('camera_position') is not None + assert elem.find('wireframe_thickness') is not None + assert elem.find('look_at') is not None + assert elem.find('pixels') is not None + assert elem.find('background').text == '0 0 0' + + def test_plots(run_in_tmpdir): p1 = openmc.Plot(name='plot1') p1.origin = (5., 5., 5.) @@ -99,10 +143,14 @@ def test_plots(run_in_tmpdir): plots = openmc.Plots([p1, p2]) assert len(plots) == 2 - p3 = openmc.Plot(name='plot3') - plots.append(p3) + p3 = openmc.ProjectionPlot(name='plot3') + plots = openmc.Plots([p1, p2, p3]) assert len(plots) == 3 + p4 = openmc.Plot(name='plot4') + plots.append(p4) + assert len(plots) == 4 + plots.export_to_xml() # from_xml