From d8fc0578bd6bf5198c8ebba1824b3459e4282bc9 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sun, 25 Feb 2024 03:20:26 -0600 Subject: [PATCH 01/10] implement Phong-shaded ray trace plots --- docs/source/io_formats/plots.rst | 101 +++- docs/source/usersguide/plots.rst | 69 +++ include/openmc/cell.h | 2 +- include/openmc/particle_data.h | 1 - include/openmc/plot.h | 236 +++++++- openmc/model/model.py | 23 +- openmc/plots.py | 428 ++++++++++---- src/particle_data.cpp | 10 +- src/plot.cpp | 532 +++++++++++++----- .../plot_projections/plots.xml | 35 ++ .../plot_projections/results_true.dat | 2 +- .../regression_tests/plot_projections/test.py | 5 +- 12 files changed, 1133 insertions(+), 311 deletions(-) diff --git a/docs/source/io_formats/plots.rst b/docs/source/io_formats/plots.rst index e6b75eafcb6..3fb54521395 100644 --- a/docs/source/io_formats/plots.rst +++ b/docs/source/io_formats/plots.rst @@ -7,13 +7,18 @@ Geometry Plotting Specification -- plots.xml Basic plotting capabilities are available in OpenMC by creating a plots.xml file and subsequently running with the ``--plot`` command-line flag. The root element of the plots.xml is simply ```` and any number output plots can be -defined with ```` sub-elements. Two plot types are currently implemented +defined with ```` sub-elements. Four plot types are currently implemented in openMC: * ``slice`` 2D pixel plot along one of the major axes. Produces a PNG image file. * ``voxel`` 3D voxel data dump. Produces an HDF5 file containing voxel xyz position and cell or material id. +* ``projection`` 2D pixel plot of a three-dimensional view of a geometry + using wireframes around cells or materials and coloring by depth through + each material. +* ``phong`` 2D pixel plot of a three-dimensional view of a geometry + with solid colored surfaces of a set of cells or materials. ------------------ @@ -66,21 +71,22 @@ sub-elements: *Default*: None - Required entry :type: - Keyword for type of plot to be produced. Currently only "slice" and "voxel" + Keyword for type of plot to be produced. Currently "slice", "voxel", + "projection", and "phong" plots are implemented. The "slice" plot type creates 2D pixel maps saved in the PNG file format. The "voxel" plot type produces a binary datafile containing voxel grid positioning and the cell or material (specified by the ``color`` tag) at the center of each voxel. Voxel plot files can be processed into VTK files using the :ref:`scripts_voxel` script provided with OpenMC and subsequently viewed with a 3D viewer such as VISIT or Paraview. - See the :ref:`io_voxel` for information about the datafile structure. + See :ref:`io_voxel` for information about the datafile structure. .. note:: High-resolution voxel files produced by OpenMC can be quite large, but the equivalent VTK files will be significantly smaller. *Default*: "slice" -```` elements of ``type`` "slice" and "voxel" must contain the ``pixels`` +All ```` elements must contain the ``pixels`` attribute or sub-element: :pixels: @@ -96,7 +102,7 @@ attribute or sub-element: ``width``/``pixels`` along that basis direction may not appear in the plot. - *Default*: None - Required entry for "slice" and "voxel" plots + *Default*: None - Required entry for all plots ```` elements of ``type`` "slice" can also contain the following attributes or sub-elements. These are not used in "voxel" plots: @@ -125,6 +131,11 @@ attributes or sub-elements. These are not used in "voxel" plots: Specifies the custom color for the cell or material. Should be 3 integers separated by spaces. + :xs: + Only for plot type "projection", the attenuation coefficient + for volume rendering of color in units of inverse + centimeters. Zero corresponds to transparency. + As an example, if your plot is colored by material and you want material 23 to be blue, the corresponding ``color`` element would look like: @@ -179,3 +190,83 @@ attributes or sub-elements. These are not used in "voxel" plots: *Default*: 0 0 0 (black) *Default*: None + +```` elements of ``type`` "projection" or "phong" can contain the following +attributes or sub-elements. + + :camera_position: + Location in 3D Cartesian space the camera is at. + + + *Default*: None - Required for all phong or projection plots + + :look_at: + Location in 3D Cartesian space the camera is looking at. + + + *Default*: None - Required for all phong or projection plots + + :field_of_view: + The horizontal field of view in degrees. Defaults to roughly + the same value as for the human eye. + + *Default*: 70 + + :orthographic_width: + If set to a nonzero value, an orthographic rather than + perspective projection for the camera is employed. + An orthographic projection puts out parallel rays from + the camera of a width prescribed here in the horizontal + direction, with the width in the vertical direction decided + by the pixel aspect ratio. + + *Default*: 0 + +```` elements of ``type`` "phong" can contain the following +attributes or sub-elements. + + :opaque_ids: + List of integer IDs of cells or materials to be treated + as visible in the plot. Whether the integers are interpreted + as cell or material IDs depends on ``color_by``. + + + *Default*: None - Required for all Phong plots + + :light_position: + Location in 3D Cartesian space of the light. + + + *Default*: Same location as ``camera_position`` + + :diffuse_fraction: + Fraction of light originating from non-directional + sources. If set to one, the coloring is not influenced + by surface curvature, and no shadows appear. If set to + zero, only regions illuminated by the light are not + black. + + + *Default*: 0.1 + +```` elements of ``type`` "projection" can contain the following +attributes or sub-elements. + + :wireframe_color: + RGB value of the wireframe's color + + *Default*: 0, 0, 0 (black) + + :wireframe_thickness: + Integer number of pixels that the wireframe takes up. + The value is a radius of the wireframe. Setting to zero + removes any wireframing. + + *Default*: 0 + + :wireframe_ids: + Integer IDs of cells or materials of regions to draw + wireframes around. Whether the integers are interpreted + as cell or material IDs depends on ``color_by``. + + *Default*: None diff --git a/docs/source/usersguide/plots.rst b/docs/source/usersguide/plots.rst index 2d588519c55..f93bff7a619 100644 --- a/docs/source/usersguide/plots.rst +++ b/docs/source/usersguide/plots.rst @@ -122,6 +122,60 @@ will depend on the 3D viewer, but should be straightforward. million or so). Thus if you want an accurate picture that renders smoothly, consider using only one voxel in a certain direction. +----------- +Phong Plots +----------- + +.. image:: ../_images/phong_triso.png + :width: 300px + +The :class:`openmc.PhongPlot` class allows three dimensional +visualization of detailed geometric features without voxelization. +The plot above visualizes a geometry created by :class:`openmc.TRISO`, +with the materials in the fuel kernel distinguished by color. It was +enclosed in a bounding box such that some kernels are cut off, +revealing the inner structure of the kernel. + +The `Phong reflection model `_ +approximates how light reflects off of a surface. On a diffusely +light-scattering material, the Phong model prescribes the amount of light +reflected from a surface as proportional to the dot product between +the normal vector of the surface and the vector between that point on +the surface and the light. With this assumption, visually appealing +plots of simulation geometries can be created. + +Phong plots use the same ray tracing functions that neutrons +and photons do in OpenMC, so any input that does not leak +particles can be visualized in 3D using a Phong plot. +That being said, these plots are not useful for detecting +overlap or undefined regions, so we recommend the slice +plot approach for geometry debugging. + +Only a few inputs are required for a Phong plot. The camera +location, where it's looking, and a set of opaque material or +cell IDs are required. The colors of materials or cells are +prescribed in the same way as slice plots. The set of IDs +which are opaque in the Phong plot must correspond to materials +if coloring by material, or cells if coloring by cell. + +A minimal Phong plot input could be: :: + + plot = openmc.PhongPlot() + plot.pixels = (600, 600) + plot.camera_position = (10.0, 20.0, -30.0) + plot.look_at = (4.0, 5.0, 1.0) + plot.color_by = 'cell' + + # optional. defaults to camera_position + plot.light_position = (10, 20, 30) + + # controls ambient lighting. Defaults to 10% + plot.diffuse_fraction = 0.1 + plot.opaque_domains = [cell2, cell3] + +These plots are then stored into a :class:`openmc.Plots` instance, +just like the slice plots. + ---------------- Projection Plots ---------------- @@ -131,6 +185,21 @@ Projection Plots .. image:: ../_images/hexlat_anim.gif :width: 200px +The :class:`openmc.ProjectionPlot` class also +produces 3D visualizations of OpenMC geometries without voxelization, +but is intended to show the inside of a model using wireframing +of cell or material boundaries in addition to cell coloring based +on path length of camera rays through the model. The coloring +in these plots is a bit like turning the model into partially +transparent colored glass that can be seen through, without any +refractive effects. This is called volume rendering. The colors +are specified in exactly the same interface employed by slice +plots. + +Similar to Phong plots, these use the native ray tracing capabilities within +OpenMC, so any geometry in which particles successfully run without overlaps +or leaks will work with projection plots. + 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 diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c405f536b5d..53afa15d4cb 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -115,7 +115,7 @@ class Region { //! //! Uses the comobination of half-spaces and binary operators to determine //! if short circuiting can be used. Short cicuiting uses the relative and - //! absolute depth of parenthases in the expression. + //! absolute depth of parentheses in the expression. bool contains_complex(Position r, Direction u, int32_t on_surface) const; //! BoundingBox if the paritcle is in a simple cell. diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index cb68fc45740..06ef822de6a 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -543,7 +543,6 @@ class ParticleData : public GeometryState { int& cell_born() { return cell_born_; } const int& cell_born() const { return cell_born_; } - // index of the current and last material // Total number of collisions suffered by particle int& n_collision() { return n_collision_; } const int& n_collision() const { return n_collision_; } diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 7b66036fb71..954458299e1 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -61,6 +61,14 @@ struct RGBColor { return red == other.red && green == other.green && blue == other.blue; } + RGBColor& operator*=(const double x) + { + red *= x; + green *= x; + blue *= x; + return *this; + } + // Members uint8_t red, green, blue; }; @@ -267,32 +275,83 @@ class Plot : public PlottableInterface, public SlicePlotBase { RGBColor meshlines_color_; //!< Color of meshlines on the plot }; -class ProjectionPlot : public PlottableInterface { - +/* + * Base class for plots which create their image by tracing + * images from a camera through the problem geometry. + */ +class RayTracePlot : public PlottableInterface { public: - ProjectionPlot(pugi::xml_node plot); + RayTracePlot(pugi::xml_node plot); + + // Standard getters. No setting since it's done from XML. + const Position& camera_position() const { return camera_position_; } + const Position& look_at() const { return look_at_; } + const double& horizontal_field_of_view() const + { + return horizontal_field_of_view_; + } - virtual void create_output() const; virtual void print_info() const; -private: + /* 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(GeometryState& p); + +protected: void set_output_path(pugi::xml_node plot_node); + + /* + * Gets the starting position and direction for the pixel corresponding + * to this horizontal and vertical position. + */ + std::pair get_pixel_ray(int horiz, int vert) const; + + // Max intersections before we assume ray tracing is caught in an infinite + // loop: + std::array pixels_; // pixel dimension of resulting image + +private: 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); + + 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 + + // TODO let up_ be set through XML + Direction up_ {0.0, 0.0, 1.0}; // which way is up + + /* The horizontal thickness, if using an orthographic projection. + * If set to zero, we assume using a perspective projection. + */ + double orthographic_width_ {0.0}; + + // Cached camera-to-model matrix + std::vector camera_to_model_; +}; + +class ProjectionRay; +class ProjectionPlot : public RayTracePlot { + + friend class ProjectionRay; + +public: + ProjectionPlot(pugi::xml_node plot); + + virtual void create_output() const; + virtual void print_info() const; + +private: + void set_opacities(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(GeometryState& 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. @@ -319,24 +378,9 @@ class ProjectionPlot : public PlottableInterface { {} }; - // 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}; @@ -344,6 +388,144 @@ class ProjectionPlot : public PlottableInterface { vector xs_; // macro cross section values for cell volume rendering }; +/* + * Plots a geometry with single-scattered Phong lighting + * plus a diffuse lighting contribution. Makes for a visually + * appealing 3D view of a geometry + */ +class PhongPlot : public RayTracePlot { + friend class PhongRay; + +public: + PhongPlot(pugi::xml_node plot); + + virtual void create_output() const; + virtual void print_info() const; + + /* All materials are completely transparent unless explicitly listed + * to be otherwise + */ + bool is_id_opaque(int id) const; + +private: + void set_opaque_ids(pugi::xml_node node); + void set_light_position(pugi::xml_node node); + void set_diffuse_fraction(pugi::xml_node node); + + // TODO + // Right now we use a sorted vector and do binary search, + // but it might be worth trying std::unordered_set, which + // uses a hash function. + std::vector opaque_ids_; + + double diffuse_fraction_ {0.1}; + + // By default, the light is at the camera unless otherwise specified. + Position light_location_; +}; + +// Base class that implements ray tracing logic, not necessarily through +// defined regions of the geometry but also outside of it. +class Ray : public GeometryState { + +public: + Ray(Position r, Direction u) { init_from_r_u(r, u); } + + // Called at every surface intersection within the model + virtual void on_intersection() = 0; + + /* + * Traces the ray through the geometry, calling on_intersection + * at every surface boundary. + */ + void trace(); + + /* + * Implementation code has read-only access to the distance + * between the ray and the next cell it's intersecting + */ + const BoundaryInfo& dist() { return dist_; } + + const int& first_surface() { return first_surface_; } + + const bool& first_inside_model() { return first_inside_model_; } + + const int& i_surface() { return i_surface_; } + + // Stops the ray and exits tracing when called from on_intersection + void stop() { stop_ = true; } + + // Sets the dist_ variable + void compute_distance(); + +private: + static const int MAX_INTERSECTIONS = 1000000; + + bool hit_something_ {false}; + bool intersection_found_ {true}; + bool first_inside_model_ {false}; + bool stop_ {false}; + + unsigned event_counter_ {0}; + + // Records the first intersected surface on the model + int first_surface_ {-1}; + int i_surface_ {-1}; + + BoundaryInfo dist_; +}; + +class ProjectionRay : public Ray { +public: + ProjectionRay(Position r, Direction u, const ProjectionPlot& plot, + vector& line_segments) + : Ray(r, u), plot_(plot), line_segments_(line_segments) + {} + + virtual void on_intersection() override; + +private: + /* Store a reference to the plot object which is running this ray, in order + * to access some of the plot settings which influence the behavior where + * intersections are. + */ + const ProjectionPlot& plot_; + + /* The ray runs through the geometry, and records the lengths of ray segments + * and cells they lie in along the way. + */ + vector& line_segments_; +}; + +class PhongRay : public Ray { +public: + PhongRay(Position r, Direction u, const PhongPlot& plot) + : Ray(r, u), plot_(plot) + { + result_color_ = plot_.not_found_; + } + + virtual void on_intersection() override; + + const RGBColor& result_color() { return result_color_; } + +private: + const PhongPlot& plot_; + + /* After the ray is reflected, it is moving towards the + * camera. It does that in order to see if the exposed surface + * is shadowed by something else. + */ + bool reflected_ {false}; + + // Have to record the first hit ID, so that if the region + // does get shadowed, we recall what its color should be + // when tracing from the surface to the light. + int orig_hit_id_ {-1}; + + RGBColor result_color_; +}; + //=============================================================================== // Non-member functions //=============================================================================== diff --git a/openmc/model/model.py b/openmc/model/model.py index 5194fec1daa..3aee4d1680f 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -146,7 +146,7 @@ def plots(self) -> Optional[openmc.Plots]: @plots.setter def plots(self, plots): - check_type('plots', plots, Iterable, openmc.Plot) + check_type('plots', plots, Iterable, openmc.PlotBase) if isinstance(plots, openmc.Plots): self._plots = plots else: @@ -238,7 +238,8 @@ def from_xml(cls, geometry='geometry.xml', materials='materials.xml', materials = openmc.Materials.from_xml(materials) geometry = openmc.Geometry.from_xml(geometry, materials) settings = openmc.Settings.from_xml(settings) - tallies = openmc.Tallies.from_xml(tallies) if Path(tallies).exists() else None + tallies = openmc.Tallies.from_xml( + tallies) if Path(tallies).exists() else None plots = openmc.Plots.from_xml(plots) if Path(plots).exists() else None return cls(geometry, materials, settings, tallies, plots) @@ -260,12 +261,16 @@ def from_model_xml(cls, path='model.xml'): model = cls() meshes = {} - model.settings = openmc.Settings.from_xml_element(root.find('settings'), meshes) - model.materials = openmc.Materials.from_xml_element(root.find('materials')) - model.geometry = openmc.Geometry.from_xml_element(root.find('geometry'), model.materials) + model.settings = openmc.Settings.from_xml_element( + root.find('settings'), meshes) + model.materials = openmc.Materials.from_xml_element( + root.find('materials')) + model.geometry = openmc.Geometry.from_xml_element( + root.find('geometry'), model.materials) if root.find('tallies') is not None: - model.tallies = openmc.Tallies.from_xml_element(root.find('tallies'), meshes) + model.tallies = openmc.Tallies.from_xml_element( + root.find('tallies'), meshes) if root.find('plots') is not None: model.plots = openmc.Plots.from_xml_element(root.find('plots')) @@ -531,11 +536,13 @@ def export_to_model_xml(self, path='model.xml', remove_surfs=False): if self.tallies: tallies_element = self.tallies.to_xml_element(mesh_memo) - xml.clean_indentation(tallies_element, level=1, trailing_indent=self.plots) + xml.clean_indentation( + tallies_element, level=1, trailing_indent=self.plots) fh.write(ET.tostring(tallies_element, encoding="unicode")) if self.plots: plots_element = self.plots.to_xml_element() - xml.clean_indentation(plots_element, level=1, trailing_indent=False) + xml.clean_indentation( + plots_element, level=1, trailing_indent=False) fh.write(ET.tostring(plots_element, encoding="unicode")) fh.write("\n") diff --git a/openmc/plots.py b/openmc/plots.py index c73e3cdf73e..f19462a6771 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -471,6 +471,15 @@ def colorize(self, geometry, seed=1): for domain in domains: self.colors[domain] = np.random.randint(0, 256, (3,)) + def _colors_to_xml(self, element): + for domain, color in sorted(self._colors.items(), + key=lambda x: self._get_id(x[0])): + subelement = ET.SubElement(element, "color") + subelement.set("id", str(self._get_id(domain))) + if isinstance(color, str): + color = _SVG_COLORS[color.lower()] + subelement.set("rgb", ' '.join(str(x) for x in color)) + def to_xml_element(self): """Save common plot attributes to XML element @@ -797,13 +806,7 @@ def to_xml_element(self): subelement.text = ' '.join(map(str, self._width)) if self._colors: - for domain, color in sorted(self._colors.items(), - key=lambda x: PlotBase._get_id(x[0])): - subelement = ET.SubElement(element, "color") - 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)) + self._colors_to_xml(element) if self._show_overlaps: subelement = ET.SubElement(element, "show_overlaps") @@ -961,7 +964,8 @@ def to_vtk(self, output: Optional[PathLike] = None, """ if self.type != 'voxel': - raise ValueError('Generating a VTK file only works for voxel plots') + raise ValueError( + 'Generating a VTK file only works for voxel plots') # Create plots.xml Plots([self]).export_to_xml(cwd) @@ -977,19 +981,16 @@ def to_vtk(self, output: Optional[PathLike] = None, return voxel_to_vtk(h5_voxel_file, output) -class ProjectionPlot(PlotBase): +class RayTracePlot(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. + This is an abstract base class that ProjectionPlot and PhongPlot finish + the implementation of. + .. versionadded:: 0.14.0 Parameters @@ -1018,21 +1019,6 @@ class ProjectionPlot(PlotBase): 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=''): @@ -1043,10 +1029,6 @@ def __init__(self, plot_id=None, name=''): 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): @@ -1100,6 +1082,151 @@ def orthographic_width(self, orthographic_width): assert orthographic_width >= 0.0 self._orthographic_width = orthographic_width + def _check_domains_consistent_with_color_by(self, domains): + """Check domains are the same as the type we are coloring by + """ + + for region in domains: + + # if an integer is passed, we have to assume + # it was a valid ID + if isinstance(region, int): + continue + + if self._color_by == 'material': + if not isinstance(region, openmc.Material): + raise Exception('Domain list must be materials if ' + 'color_by=material') + else: + if not isinstance(region, openmc.Cell): + raise Exception('Domain list must be materials if ' + 'color_by=cell') + + def to_xml_element(self): + """Return XML representation of the ray trace plot + + Returns + ------- + element : lxml.etree._Element + XML element containing plot data + + """ + + element = super().to_xml_element() + element.set("id", str(self._id)) + + 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, "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 = '' + 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('\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('\tLevel', self._level) + return string + + def _read_xml_attributes(self, elem): + """Helper function called by from_xml_element + of child classes. These are common vaues to be + read by any ray traced plot. + + Returns + ------- + None + """ + + if "filename" in elem.keys(): + self.filename = elem.get("filename") + self.color_by = elem.get("color_by") + + horizontal_fov = elem.find("horizontal_field_of_view") + if horizontal_fov is not None: + self.horizontal_field_of_view = float(horizontal_fov.text) + + tmp = elem.find("orthographic_width") + if tmp is not None: + self.orthographic_width = float(tmp) + + self.pixels = get_elem_tuple(elem, "pixels") + self.camera_position = get_elem_tuple(elem, "camera_position", float) + self.look_at = get_elem_tuple(elem, "look_at", float) + + # 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: + self.mask_background = tuple( + [int(x) for x in background.split()]) + + # Set universe level + level = elem.find("level") + if level is not None: + self.level = int(level.text) + + +class ProjectionPlot(RayTracePlot): + """Plots wireframes of geometry with volume rendered colors + + 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. + + 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=''): + super().__init__(plot_id, name) + self._wireframe_thickness = 1 + self._wireframe_color = _SVG_COLORS['black'] + self._wireframe_domains = [] + self._xs = {} + @property def wireframe_thickness(self): return self._wireframe_thickness @@ -1126,15 +1253,6 @@ def wireframe_domains(self): @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 @property @@ -1171,6 +1289,18 @@ def set_transparent(self, geometry): for domain in domains: self.xs[domain] = 0.0 + def __repr__(self): + string = 'Projection Plot\n' + string += super().__repr__() + 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('\tTransparencies', self._xs) + return string + def to_xml_element(self): """Return XML representation of the projection plot @@ -1180,16 +1310,9 @@ def to_xml_element(self): XML element containing plot data """ - element = super().to_xml_element() 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) @@ -1199,6 +1322,8 @@ def to_xml_element(self): color = _SVG_COLORS[color.lower()] subelement.text = ' '.join(str(x) for x in color) + self._check_domains_consistent_with_color_by(self._wireframe_domains) + if self._wireframe_domains: id_list = [x.id for x in self._wireframe_domains] subelement = ET.SubElement(element, "wireframe_ids") @@ -1216,43 +1341,8 @@ def to_xml_element(self): 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 @@ -1268,24 +1358,12 @@ def from_xml_element(cls, elem): 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: - plot.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) + plot._read_xml_attributes(elem) # Attempt to get wireframe thickness. May not be present wireframe_thickness = elem.get("wireframe_thickness") @@ -1296,30 +1374,128 @@ def from_xml_element(cls, elem): plot.wireframe_color = [int(item) for item in wireframe_color] # Set plot colors - colors = {} - xs = {} + plot._colors = {} + plot._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")) + plot._colors[uid] = get_elem_tuple(color_elem, "rgb") + plot._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()]) + return plot - # Set universe level - level = elem.find("level") - if level is not None: - plot.level = int(level.text) - return plot +class PhongPlot(RayTracePlot): + + def __init__(self, plot_id=None, name=''): + super().__init__(plot_id, name) + self._light_position = None + self._diffuse_fraction = 0.1 + self._opaque_domains = [] + + @property + def light_position(self): + return self._light_position + + @light_position.setter + def light_position(self, x): + cv.check_type('plot light position', x, Iterable, Real) + cv.check_length('plot light position', x, 3) + self._light_position = x + + @property + def diffuse_fraction(self): + return self._diffuse_fraction + + @diffuse_fraction.setter + def diffuse_fraction(self, x): + cv.check_type('diffuse fraction', x, Real) + cv.check_greater_than('diffuse fraction', x, 0.0, equality=True) + cv.check_less_than('diffuse fraction', x, 1.0, equality=True) + self._diffuse_fraction = x + + @property + def opaque_domains(self): + return self._opaque_domains + + @opaque_domains.setter + def opaque_domains(self, x): + cv.check_type('opaque domains', x, Iterable) + self._opaque_domains = x + + def __repr__(self): + string = 'Phong Plot\n' + string += super().__repr__() + string += '{: <16}=\t{}\n'.format('\tDiffuse Fraction', + self._diffuse_fraction) + string += '{: <16}=\t{}\n'.format('\tLight position', + self._light_position) + string += '{: <16}=\t{}\n'.format('\tOpaque domains', + self._opaque_domains) + return string + + def to_xml_element(self): + """Return XML representation of the Phong plot + + Returns + ------- + element : lxml.etree._Element + XML element containing plot data + + """ + element = super().to_xml_element() + element.set("type", "phong") + + # no light position means put it at the camera + if self._light_position: + subelement = ET.SubElement(element, "light_position") + subelement.text = ' '.join(map(str, self._light_position)) + + # no diffuse fraction defaults to 0.1 + if self._diffuse_fraction: + subelement = ET.SubElement(element, "diffuse_fraction") + subelement.text = str(self._diffuse_fraction) + + self._check_domains_consistent_with_color_by(self._opaque_domains) + subelement = ET.SubElement(element, "opaque_ids") + + # Extract all IDs, or use the integer value passed in + # explicitly if that was given + subelement.text = ' '.join( + [str(domain) if isinstance(domain, int) else + str(domain.id) for domain in self._opaque_domains]) + + if self._colors: + self._colors_to_xml(element) + + return element + + @classmethod + def from_xml_element(cls, elem): + """Generate plot object from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + + Returns + ------- + openmc.ProjectionPlot + ProjectionPlot object + + """ + + plot_id = int(elem.get("id")) + plot = cls(plot_id) + plot.type = "phong" + + plot._read_xml_attributes(elem) + + # Set plot colors + plot._colors = {} + for color_elem in elem.findall("color"): + uid = color_elem.get("id") + plot._colors[uid] = get_elem_tuple(color_elem, "rgb") class Plots(cv.CheckedList): @@ -1339,13 +1515,14 @@ class Plots(cv.CheckedList): Parameters ---------- - plots : Iterable of openmc.Plot or openmc.ProjectionPlot + plots : Iterable of openmc.Plot, openmc.ProjectionPlot, + or openmc.PhongPlot plots to add to the collection """ def __init__(self, plots=None): - super().__init__((Plot, ProjectionPlot), 'plots collection') + super().__init__((Plot, ProjectionPlot, PhongPlot), 'plots collection') self._plots_file = ET.Element("plots") if plots is not None: self += plots @@ -1355,7 +1532,8 @@ def append(self, plot): Parameters ---------- - plot : openmc.Plot or openmc.ProjectionPlot + plot : openmc.Plot, openmc.ProjectionPlot, or + openmc.PhongPlot Plot to append """ @@ -1489,8 +1667,14 @@ def from_xml_element(cls, elem): plot_type = e.get('type') if plot_type == 'projection': plots.append(ProjectionPlot.from_xml_element(e)) - else: + elif plot_type == 'phong': + plots.append(PhongPlot.from_xml_element(e)) + elif plot_type == 'voxel': plots.append(Plot.from_xml_element(e)) + elif plot_type == 'slice': + plots.append(Plot.from_xml_element(e)) + else: + raise ValueError("Unknown plot type: {}".format(plot_type)) return plots @classmethod diff --git a/src/particle_data.cpp b/src/particle_data.cpp index d8a5bb9d99c..065f8be4265 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -15,6 +15,11 @@ namespace openmc { +void GeometryState::mark_as_lost(const char* message) +{ + fatal_error(message); +} + void GeometryState::mark_as_lost(const std::string& message) { mark_as_lost(message.c_str()); @@ -25,11 +30,6 @@ void GeometryState::mark_as_lost(const std::stringstream& message) mark_as_lost(message.str()); } -void GeometryState::mark_as_lost(const char* message) -{ - fatal_error(message); -} - void LocalCoord::rotate(const vector& rotation) { r = r.rotate(rotation); diff --git a/src/plot.cpp b/src/plot.cpp index bf733cff48f..d303d780940 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -210,6 +210,8 @@ void read_plots_xml(pugi::xml_node root) std::make_unique(node, Plot::PlotType::voxel)); else if (type_str == "projection") model::plots.emplace_back(std::make_unique(node)); + else if (type_str == "phong") + model::plots.emplace_back(std::make_unique(node)); else fatal_error( fmt::format("Unsupported plot type '{}' in plot {}", type_str, id)); @@ -1029,23 +1031,45 @@ RGBColor random_color(void) int(prn(&model::plotter_seed) * 255), int(prn(&model::plotter_seed) * 255)}; } -ProjectionPlot::ProjectionPlot(pugi::xml_node node) : PlottableInterface(node) +RayTracePlot::RayTracePlot(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); + set_output_path(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."); + + // 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(); + + // Cache the camera-to-model matrix + 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}; +} + +ProjectionPlot::ProjectionPlot(pugi::xml_node node) : RayTracePlot(node) +{ + set_opacities(node); + set_wireframe_thickness(node); + set_wireframe_ids(node); + set_wireframe_color(node); } void ProjectionPlot::set_wireframe_color(pugi::xml_node plot_node) @@ -1061,7 +1085,7 @@ void ProjectionPlot::set_wireframe_color(pugi::xml_node plot_node) } } -void ProjectionPlot::set_output_path(pugi::xml_node node) +void RayTracePlot::set_output_path(pugi::xml_node node) { // Set output file path std::string filename; @@ -1085,10 +1109,9 @@ void ProjectionPlot::set_output_path(pugi::xml_node node) // 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(GeometryState& p) +int RayTracePlot::advance_to_boundary_from_void(GeometryState& p) { - constexpr double scoot = 1e-5; - double min_dist = {INFINITY}; + double min_dist {INFINITY}; auto coord = p.coord(0); Universe* uni = model::universes[model::root_universe].get(); int intersected_surface = -1; @@ -1103,7 +1126,7 @@ int ProjectionPlot::advance_to_boundary_from_void(GeometryState& p) 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; + p.coord(j).r += (min_dist + TINY_BIT) * p.coord(j).u; return std::abs(intersected_surface); } } @@ -1169,26 +1192,9 @@ bool ProjectionPlot::trackstack_equivalent( } } -void ProjectionPlot::create_output() const +std::pair RayTracePlot::get_pixel_ray( + int horiz, int vert) 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. @@ -1197,9 +1203,47 @@ void ProjectionPlot::create_output() const 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; + // focal_plane_dist can be changed to alter the perspective distortion + // effect. This is in units of cm. This seems to look good most of the + // time. TODO let this variable be set through XML. + constexpr double focal_plane_dist = 10.0; + const double dx = 2.0 * focal_plane_dist * std::tan(0.5 * horiz_fov_radians); + const double dy = p1 / p0 * dx; + + std::pair result; + + // Generate the starting position/direction of the ray + if (orthographic_width_ == 0.0) { // perspective projection + Direction camera_local_vec; + camera_local_vec.x = focal_plane_dist; + camera_local_vec.y = -0.5 * dx + horiz * dx / p0; + camera_local_vec.z = 0.5 * dy - vert * dy / p1; + camera_local_vec /= camera_local_vec.norm(); + + result.first = camera_position_; + result.second = camera_local_vec.rotate(camera_to_model_); + } else { // orthographic projection + + double x_pix_coord = (static_cast(horiz) - p0 / 2.0) / p0; + double y_pix_coord = (static_cast(vert) - p1 / 2.0) / p1; + + Direction yaxis = { + camera_to_model_[1], camera_to_model_[4], camera_to_model_[7]}; + Direction zaxis = { + camera_to_model_[2], camera_to_model_[5], camera_to_model_[8]}; + result.first = camera_position_ + + yaxis * x_pix_coord * orthographic_width_ + + zaxis * y_pix_coord * orthographic_width_; + result.second = { + camera_to_model_[0], camera_to_model_[3], camera_to_model_[6]}; + } + + return result; +} + +void ProjectionPlot::create_output() const +{ size_t width = pixels_[0]; size_t height = pixels_[1]; ImageData data({width, height}, not_found_); @@ -1235,9 +1279,6 @@ void ProjectionPlot::create_output() const const int n_threads = num_threads(); const int tid = thread_num(); - GeometryState p; - p.u() = {1.0, 0.0, 0.0}; - int vert = tid; for (int iter = 0; iter <= pixels_[1] / n_threads; iter++) { @@ -1252,105 +1293,14 @@ void ProjectionPlot::create_output() const for (int horiz = 0; horiz < pixels_[0]; ++horiz) { - // Projection mode below decides ray starting conditions - Position init_r; - Direction init_u; - - // 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); - init_u = camera_local_vec.rotate(camera_to_model); - init_r = camera_position_; - } else { // orthographic projection - init_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; - - init_r = camera_position_; - init_r += cam_yaxis * x_pix_coord * orthographic_width_; - init_r += cam_zaxis * y_pix_coord * orthographic_width_; - } - - // Resets internal geometry state of particle - p.init_from_r_u(init_r, init_u); - - bool hitsomething = false; - bool intersection_found = true; - int loop_counter = 0; + // RayTracePlot implements camera ray generation + std::pair ru = get_pixel_ray(horiz, vert); this_line_segments[tid][horiz].clear(); + ProjectionRay ray( + ru.first, ru.second, *this, this_line_segments[tid][horiz]); - 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.lowest_coord().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.lowest_coord().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"); - } + ray.trace(); // Now color the pixel based on what we have intersected... // Loops backwards over intersections. @@ -1444,9 +1394,8 @@ void ProjectionPlot::create_output() const #endif } -void ProjectionPlot::print_info() const +void RayTracePlot::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); @@ -1455,6 +1404,12 @@ void ProjectionPlot::print_info() const fmt::print("Pixels: {} {}\n", pixels_[0], pixels_[1]); } +void ProjectionPlot::print_info() const +{ + fmt::print("Plot Type: Projection\n"); + RayTracePlot::print_info(); +} + void ProjectionPlot::set_opacities(pugi::xml_node node) { xs_.resize(colors_.size(), 1e6); // set to large value for opaque by default @@ -1485,7 +1440,7 @@ void ProjectionPlot::set_opacities(pugi::xml_node node) } } -void ProjectionPlot::set_orthographic_width(pugi::xml_node node) +void RayTracePlot::set_orthographic_width(pugi::xml_node node) { if (check_for_node(node, "orthographic_width")) { double orthographic_width = @@ -1522,7 +1477,7 @@ void ProjectionPlot::set_wireframe_ids(pugi::xml_node node) std::sort(wireframe_ids_.begin(), wireframe_ids_.end()); } -void ProjectionPlot::set_pixels(pugi::xml_node node) +void RayTracePlot::set_pixels(pugi::xml_node node) { vector pxls = get_node_array(node, "pixels"); if (pxls.size() != 2) @@ -1532,19 +1487,19 @@ void ProjectionPlot::set_pixels(pugi::xml_node node) pixels_[1] = pxls[1]; } -void ProjectionPlot::set_camera_position(pugi::xml_node node) +void RayTracePlot::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")); + fatal_error(fmt::format( + "camera_position 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) +void RayTracePlot::set_look_at(pugi::xml_node node) { vector look_at = get_node_array(node, "look_at"); if (look_at.size() != 3) { @@ -1555,7 +1510,7 @@ void ProjectionPlot::set_look_at(pugi::xml_node node) look_at_.z = look_at[2]; } -void ProjectionPlot::set_field_of_view(pugi::xml_node node) +void RayTracePlot::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")) { @@ -1569,6 +1524,303 @@ void ProjectionPlot::set_field_of_view(pugi::xml_node node) } } +PhongPlot::PhongPlot(pugi::xml_node node) : RayTracePlot(node) +{ + set_opaque_ids(node); + set_diffuse_fraction(node); + set_light_position(node); +} + +void PhongPlot::print_info() const +{ + fmt::print("Plot Type: Phong\n"); + RayTracePlot::print_info(); +} + +void PhongPlot::create_output() const +{ + size_t width = pixels_[0]; + size_t height = pixels_[1]; + ImageData data({width, height}, not_found_); + +#pragma omp parallel for schedule(dynamic) collapse(2) + for (int horiz = 0; horiz < pixels_[0]; ++horiz) { + for (int vert = 0; vert < pixels_[1]; ++vert) { + + // RayTracePlot implements camera ray generation + std::pair ru = get_pixel_ray(horiz, vert); + PhongRay ray(ru.first, ru.second, *this); + + ray.trace(); + data(horiz, vert) = ray.result_color(); + } + } + +#ifdef USE_LIBPNG + output_png(path_plot(), data); +#else + output_ppm(path_plot(), data); +#endif +} + +void PhongPlot::set_opaque_ids(pugi::xml_node node) +{ + if (check_for_node(node, "opaque_ids")) { + opaque_ids_ = get_node_array(node, "opaque_ids"); + // It is read in as actual ID values, but we have to convert to indices in + // mat/cell array + for (auto& x : opaque_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(opaque_ids_.begin(), opaque_ids_.end()); +} + +bool PhongPlot::is_id_opaque(int id) const +{ + return std::binary_search(opaque_ids_.begin(), opaque_ids_.end(), id); +} + +void PhongPlot::set_light_position(pugi::xml_node node) +{ + if (check_for_node(node, "light_position")) { + auto light_pos_tmp = get_node_array(node, "light_position"); + + if (light_pos_tmp.size() != 3) + fatal_error("Light position must be given as 3D coordinates"); + + light_location_.x = light_pos_tmp[0]; + light_location_.y = light_pos_tmp[1]; + light_location_.z = light_pos_tmp[2]; + } else { + light_location_ = camera_position(); + } +} + +void PhongPlot::set_diffuse_fraction(pugi::xml_node node) +{ + if (check_for_node(node, "diffuse_fraction")) { + diffuse_fraction_ = std::stod(get_node_value(node, "diffuse_fraction")); + if (diffuse_fraction_ < 0.0 || diffuse_fraction_ > 1.0) { + fatal_error("Must have 0<=diffuse fraction<= 1"); + } + } +} + +void Ray::compute_distance() +{ + dist_ = distance_to_boundary(*this); +} + +void Ray::trace() +{ + int first_surface_ = -1; // surface first passed when entering the model + first_inside_model_ = true; // false after entering the model + while (intersection_found_) { + bool inside_cell = false; + + i_surface_ = std::abs(surface()) - 1; + + // This means no surface was intersected. See cell.cpp + // and search for numeric_limits to see where we return it. + if (surface() == std::numeric_limits::max()) { + warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u())); + return; + } + + if (i_surface_ > 0 && + model::surfaces[i_surface_]->geom_type_ == GeometryType::DAG) { +#ifdef DAGMC + int32_t i_cell = next_cell( + i_surface_, cell_last(n_coord() - 1), 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(*this, settings::verbosity >= 10); + } + + if (inside_cell) { + + hit_something_ = true; + intersection_found_ = true; + + compute_distance(); + + if (first_inside_model_) { + i_surface_ = first_surface_ - 1; + first_inside_model_ = false; + } + + // Call the specialized logic for this type of ray + on_intersection(); + if (stop_) + return; + + // There are no more intersections to process + // if we hit the edge of the model, so stop + // the particle in that case: + if (dist_.distance == INFTY || dist_.distance == INFINITY) { + return; + } + + // Advance particle, prepare for next intersection + for (int lev = 0; lev < n_coord(); ++lev) { + coord(lev).r += dist_.distance * coord(lev).u; + } + surface() = dist_.surface_index; + n_coord_last() = n_coord(); + n_coord() = dist_.coord_level; + if (dist_.lattice_translation[0] != 0 || + dist_.lattice_translation[1] != 0 || + dist_.lattice_translation[2] != 0) { + cross_lattice(*this, dist_, settings::verbosity >= 10); + } + + } else { + first_surface_ = RayTracePlot::advance_to_boundary_from_void(*this); + intersection_found_ = first_surface_ != -1; // -1 if no surface found + } + event_counter_++; + if (event_counter_ > MAX_INTERSECTIONS) + fatal_error("Infinite loop in ray traced plot"); + } +} + +void ProjectionRay::on_intersection() +{ + // This allows drawing wireframes with surface intersection + // edges on the model boundary for the same cell. + if (first_inside_model()) { + line_segments_.emplace_back( + plot_.color_by_ == PlottableInterface::PlotColorBy::mats + ? material() + : lowest_coord().cell, + 0.0, first_surface()); + } + + // This records a tuple with the following info + // + // 1) ID (material or cell depending on color_by_) + // 2) Distance traveled by the ray through that ID + // 3) Index of the intersected surface (starting from 1) + line_segments_.emplace_back( + plot_.color_by_ == PlottableInterface::PlotColorBy::mats + ? material() + : lowest_coord().cell, + dist().distance, std::abs(dist().surface_index)); +} + +void PhongRay::on_intersection() +{ + // Check if we hit an opaque material or cell + int hit_id = plot_.color_by_ == PlottableInterface::PlotColorBy::mats + ? material() + : lowest_coord().cell; + + // If we are reflected and have advanced beyond the camera, + // the ray is done. This is checked here because we should + // kill the ray even if the material is not opaque. + if (reflected_ && (r() - plot_.camera_position()).dot(u()) >= 0.0) { + stop(); + return; + } + + // Anything that's not opaque has zero impact on the plot. + if (!plot_.is_id_opaque(hit_id)) + return; + + if (!reflected_) { + + // reflect the particle and set the color to be colored by + // the normal or the diffuse lighting contribution + reflected_ = true; + result_color_ = plot_.colors_[hit_id]; + Direction to_light = plot_.light_location_ - r(); + to_light /= to_light.norm(); + + // Not sure what can cause this error. Color as an error and proceed + // from there. It seems to happen only for a few pixels on the outer + // boundary of a hex lattice. TODO + // + // We cannot detect it in the outer loop, and it only matters here, so + // that's why the error handling is a little different than for a lost + // ray. + if (i_surface() == -1) { + result_color_ = plot_.overlap_color_; + stop(); + return; + } + + Direction normal = model::surfaces.at(i_surface())->normal(r_local()); + normal /= normal.norm(); + + // Need to apply translations to find the normal vector in + // the base level universe's coordinate system. Uses fact + // that rotation matrices are orthonormal. + auto inverse_rotate = [](const Direction u, + const std::vector& rotation) { + return Direction { + u.x * rotation[0] + u.y * rotation[3] + u.z * rotation[6], + u.x * rotation[1] + u.y * rotation[4] + u.z * rotation[7], + u.x * rotation[2] + u.y * rotation[5] + u.z * rotation[8]}; + }; + for (int lev = n_coord() - 2; lev >= 0; --lev) { + if (coord(lev + 1).rotated) { + const Cell& c {*model::cells[coord(lev).cell]}; + normal = inverse_rotate(normal, c.rotation_); + } + } + + if (surface() > 0) { + normal *= -1.0; + } + + // Facing away from the light means no lighting + double dotprod = normal.dot(to_light); + dotprod = dotprod >= 0.0 ? dotprod : 0.0; + + double modulation = + plot_.diffuse_fraction_ + (1.0 - plot_.diffuse_fraction_) * dotprod; + result_color_ *= modulation; + + // Now point the particle to the camera. We now begin + // checking to see if it's occluded by another surface + u() = to_light; + + orig_hit_id_ = hit_id; + + surface() = -surface(); // go to other side + + // Must fully restart coordinate search. Why? Not sure. + clear(); + + bool found = exhaustive_find_cell(*this); + if (!found) { + fatal_error("Lost particle after reflection."); + } + + // Must recalculate distance to boundary due to the + // direction change + compute_distance(); + + } else { + // If it's not facing the light, we color with the diffuse contribution, so + // next we check if we're going to occlude the last reflected surface. if + // so, color by the diffuse contribution instead + + if (orig_hit_id_ == -1) + fatal_error("somehow a ray got reflected but not original ID set?"); + + result_color_ = plot_.colors_[orig_hit_id_]; + result_color_ *= plot_.diffuse_fraction_; + stop(); + } +} + extern "C" int openmc_id_map(const void* plot, int32_t* data_out) { diff --git a/tests/regression_tests/plot_projections/plots.xml b/tests/regression_tests/plot_projections/plots.xml index 85689da1475..6e1311286b5 100644 --- a/tests/regression_tests/plot_projections/plots.xml +++ b/tests/regression_tests/plot_projections/plots.xml @@ -52,4 +52,39 @@ + + 0. 0. 0. + 10. 10. 10. + 200 200 + phong.png + 1 3 + + + + + + + 0. 0. 0. + 10. 10. 10. + 0.5 + 200 200 + phong_diffuse.png + 1 3 + + + + + + + 0. 0. 0. + 10. 10. 10. + 0. 10. 10. + 200 200 + phong_move_light.png + 1 3 + + + + + diff --git a/tests/regression_tests/plot_projections/results_true.dat b/tests/regression_tests/plot_projections/results_true.dat index 1a67da3002a..8fad4a78527 100644 --- a/tests/regression_tests/plot_projections/results_true.dat +++ b/tests/regression_tests/plot_projections/results_true.dat @@ -1 +1 @@ -24fb0f41ee018ea086962dbd6bcd0b536d11d4b34644bfef4f0e74f8b462fe41a84af39c7ff79046d5d7cfe209084eac54712fa0ec01038e97eb43df1abd0334 \ No newline at end of file +c15c592ea100dcb58a4ac7e665b67e34c69b463e65c96dd561ddfa65a7997370cbe875453f697297ec23776dbd5d47af0a789beb07c36397e450d2285021b17b \ No newline at end of file diff --git a/tests/regression_tests/plot_projections/test.py b/tests/regression_tests/plot_projections/test.py index cf18503f4b8..37a6ecf28d8 100644 --- a/tests/regression_tests/plot_projections/test.py +++ b/tests/regression_tests/plot_projections/test.py @@ -2,5 +2,8 @@ from tests.regression_tests import config def test_plot(): - harness = PlotTestHarness(('plot_1.png', 'example1.png', 'example2.png', 'example3.png', 'orthographic_example1.png')) + harness = PlotTestHarness(('plot_1.png', 'example1.png', 'example2.png', + 'example3.png', 'orthographic_example1.png', + 'phong.png', 'phong_diffuse.png', + 'phong_move_light.png')) harness.main() From b02d5863955967f82116d9b309b539599382731c Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 11:12:39 -0600 Subject: [PATCH 02/10] remove outdated TODO --- include/openmc/plot.h | 1 - 1 file changed, 1 deletion(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 954458299e1..eafaaf7987c 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -323,7 +323,6 @@ class RayTracePlot : public PlottableInterface { Position camera_position_; // where camera is Position look_at_; // point camera is centered looking at - // TODO let up_ be set through XML Direction up_ {0.0, 0.0, 1.0}; // which way is up /* The horizontal thickness, if using an orthographic projection. From efadc5eac83e69a3a05fa68ab9045fc7d35ff268 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 11:31:24 -0600 Subject: [PATCH 03/10] use std::set for opaque ids in plot --- include/openmc/plot.h | 12 ++---------- src/plot.cpp | 17 ++++++----------- 2 files changed, 8 insertions(+), 21 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index eafaaf7987c..ad71adc722f 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -2,6 +2,7 @@ #define OPENMC_PLOT_H #include +#include #include #include @@ -401,21 +402,12 @@ class PhongPlot : public RayTracePlot { virtual void create_output() const; virtual void print_info() const; - /* All materials are completely transparent unless explicitly listed - * to be otherwise - */ - bool is_id_opaque(int id) const; - private: void set_opaque_ids(pugi::xml_node node); void set_light_position(pugi::xml_node node); void set_diffuse_fraction(pugi::xml_node node); - // TODO - // Right now we use a sorted vector and do binary search, - // but it might be worth trying std::unordered_set, which - // uses a hash function. - std::vector opaque_ids_; + std::set opaque_ids_; double diffuse_fraction_ {0.1}; diff --git a/src/plot.cpp b/src/plot.cpp index d303d780940..55e33394bb1 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1566,21 +1566,16 @@ void PhongPlot::create_output() const void PhongPlot::set_opaque_ids(pugi::xml_node node) { if (check_for_node(node, "opaque_ids")) { - opaque_ids_ = get_node_array(node, "opaque_ids"); + auto opaque_ids_tmp = get_node_array(node, "opaque_ids"); + // It is read in as actual ID values, but we have to convert to indices in // mat/cell array - for (auto& x : opaque_ids_) + for (auto& x : opaque_ids_tmp) 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(opaque_ids_.begin(), opaque_ids_.end()); -} -bool PhongPlot::is_id_opaque(int id) const -{ - return std::binary_search(opaque_ids_.begin(), opaque_ids_.end(), id); + opaque_ids_.insert(opaque_ids_tmp.begin(), opaque_ids_tmp.end()); + } } void PhongPlot::set_light_position(pugi::xml_node node) @@ -1730,7 +1725,7 @@ void PhongRay::on_intersection() } // Anything that's not opaque has zero impact on the plot. - if (!plot_.is_id_opaque(hit_id)) + if (plot_.opaque_ids_.find(hit_id) == plot_.opaque_ids_.end()) return; if (!reflected_) { From bc43061249af0583a631330893e64a3fce9a9392 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 11:56:01 -0600 Subject: [PATCH 04/10] use doxygen class comments for new plot types --- include/openmc/plot.h | 48 ++++++++++++++++++++++++++++++++----------- src/plot.cpp | 4 ++-- 2 files changed, 38 insertions(+), 14 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index ad71adc722f..bfadb7a8fdb 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -79,9 +79,14 @@ 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 + * @brief Interface for plottable objects. + * + * PlottableInterface classes must have a unique ID in the plots.xml file. + * They guarantee the ability to create output in some form. This interface + * is designed to be implemented by classes that produce plot-relevant data + * which can be visualized. */ class PlottableInterface { private: @@ -240,8 +245,8 @@ T SlicePlotBase::get_map() const data.set_overlap(y, x); } } // inner for - } // outer for - } // omp parallel + } // outer for + } // omp parallel return data; } @@ -276,9 +281,14 @@ class Plot : public PlottableInterface, public SlicePlotBase { RGBColor meshlines_color_; //!< Color of meshlines on the plot }; -/* - * Base class for plots which create their image by tracing - * images from a camera through the problem geometry. +/** + * @class RaytracePlot + * @brief Base class for plots that generate images through ray tracing. + * + * This class serves as a base for plots that create their visuals by tracing rays + * from a camera through the problem geometry. It inherits from PlottableInterface, + * ensuring that it provides an implementation for generating output specific to + * ray-traced visualization. */ class RayTracePlot : public PlottableInterface { public: @@ -336,6 +346,16 @@ class RayTracePlot : public PlottableInterface { }; class ProjectionRay; + +/** + * @class ProjectionPlot + * @brief Creates plots that are like colorful x-ray imaging + * + * ProjectionPlot is a specialized form of RayTracePlot designed for creating + * projection plots. This involves tracing rays from a camera through the problem + * geometry and rendering the results based on depth of penetration through materials + * or cells and their colors. + */ class ProjectionPlot : public RayTracePlot { friend class ProjectionRay; @@ -388,10 +408,14 @@ class ProjectionPlot : public RayTracePlot { vector xs_; // macro cross section values for cell volume rendering }; -/* + +/** + * @class PhongPlot + * @brief Plots 3D objects as the eye might see them. + * * Plots a geometry with single-scattered Phong lighting - * plus a diffuse lighting contribution. Makes for a visually - * appealing 3D view of a geometry + * plus a diffuse lighting contribution. The result is a + * physically reasonable, aesthetic 3D view of a geometry. */ class PhongPlot : public RayTracePlot { friend class PhongRay; @@ -407,7 +431,7 @@ class PhongPlot : public RayTracePlot { void set_light_position(pugi::xml_node node); void set_diffuse_fraction(pugi::xml_node node); - std::set opaque_ids_; + std::set opaque_ids_; double diffuse_fraction_ {0.1}; diff --git a/src/plot.cpp b/src/plot.cpp index 55e33394bb1..37a406ccbda 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -264,8 +264,8 @@ void Plot::create_image() const } data(x, y) = colors_[model::material_map[id]]; } // color_by if-else - } // x for loop - } // y for loop + } // x for loop + } // y for loop // draw mesh lines if present if (index_meshlines_mesh_ >= 0) { From 913571eaaf423f83f9bd8466984fef73609ec703 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 11:57:45 -0600 Subject: [PATCH 05/10] clang format --- include/openmc/plot.h | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index bfadb7a8fdb..61e9b729111 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -285,10 +285,10 @@ class Plot : public PlottableInterface, public SlicePlotBase { * @class RaytracePlot * @brief Base class for plots that generate images through ray tracing. * - * This class serves as a base for plots that create their visuals by tracing rays - * from a camera through the problem geometry. It inherits from PlottableInterface, - * ensuring that it provides an implementation for generating output specific to - * ray-traced visualization. + * This class serves as a base for plots that create their visuals by tracing + * rays from a camera through the problem geometry. It inherits from + * PlottableInterface, ensuring that it provides an implementation for + * generating output specific to ray-traced visualization. */ class RayTracePlot : public PlottableInterface { public: @@ -352,9 +352,9 @@ class ProjectionRay; * @brief Creates plots that are like colorful x-ray imaging * * ProjectionPlot is a specialized form of RayTracePlot designed for creating - * projection plots. This involves tracing rays from a camera through the problem - * geometry and rendering the results based on depth of penetration through materials - * or cells and their colors. + * projection plots. This involves tracing rays from a camera through the + * problem geometry and rendering the results based on depth of penetration + * through materials or cells and their colors. */ class ProjectionPlot : public RayTracePlot { @@ -408,7 +408,6 @@ class ProjectionPlot : public RayTracePlot { vector xs_; // macro cross section values for cell volume rendering }; - /** * @class PhongPlot * @brief Plots 3D objects as the eye might see them. From eb32d8f6ad860d9f77aa95a22ad12b358e81aebf Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 12:10:46 -0600 Subject: [PATCH 06/10] add application of inverse rotation method to Position class --- include/openmc/plot.h | 4 ++-- include/openmc/position.h | 6 +++++- src/plot.cpp | 28 +++++++++++----------------- src/position.cpp | 7 +++++++ 4 files changed, 25 insertions(+), 20 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 61e9b729111..3e1389e1539 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -319,8 +319,6 @@ class RayTracePlot : public PlottableInterface { */ std::pair get_pixel_ray(int horiz, int vert) const; - // Max intersections before we assume ray tracing is caught in an infinite - // loop: std::array pixels_; // pixel dimension of resulting image private: @@ -473,6 +471,8 @@ class Ray : public GeometryState { void compute_distance(); private: + // Max intersections before we assume ray tracing is caught in an infinite + // loop: static const int MAX_INTERSECTIONS = 1000000; bool hit_something_ {false}; diff --git a/include/openmc/position.h b/include/openmc/position.h index 11ea3764792..4e72588a44b 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -94,9 +94,13 @@ struct Position { //! \result Reflected vector Position reflect(Position n) const; - //! Rotate the position based on a rotation matrix + //! Rotate the position by applying a rotation matrix Position rotate(const vector& rotation) const; + //! Rotate the position by aplpying the inverse of a rotation matrix + //! using the fact that rotation matrices are orthonormal. + Position inverse_rotate(const vector& rotation) const; + // Data members double x = 0.; double y = 0.; diff --git a/src/plot.cpp b/src/plot.cpp index 37a406ccbda..f0d8d3fae3e 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1546,11 +1546,9 @@ void PhongPlot::create_output() const #pragma omp parallel for schedule(dynamic) collapse(2) for (int horiz = 0; horiz < pixels_[0]; ++horiz) { for (int vert = 0; vert < pixels_[1]; ++vert) { - // RayTracePlot implements camera ray generation std::pair ru = get_pixel_ray(horiz, vert); PhongRay ray(ru.first, ru.second, *this); - ray.trace(); data(horiz, vert) = ray.result_color(); } @@ -1680,8 +1678,10 @@ void Ray::trace() intersection_found_ = first_surface_ != -1; // -1 if no surface found } event_counter_++; - if (event_counter_ > MAX_INTERSECTIONS) - fatal_error("Infinite loop in ray traced plot"); + if (event_counter_ > MAX_INTERSECTIONS) { + warning("Likely infinite loop in ray traced plot"); + return; + } } } @@ -1737,9 +1737,11 @@ void PhongRay::on_intersection() Direction to_light = plot_.light_location_ - r(); to_light /= to_light.norm(); - // Not sure what can cause this error. Color as an error and proceed - // from there. It seems to happen only for a few pixels on the outer - // boundary of a hex lattice. TODO + // TODO + // Not sure what can cause i_surface()==-1, although it sometimes happens + // for a few pixels. It's very very rare, so proceed by coloring the pixel + // with the overlap color. It seems to happen only for a few pixels on the + // outer boundary of a hex lattice. // // We cannot detect it in the outer loop, and it only matters here, so // that's why the error handling is a little different than for a lost @@ -1754,19 +1756,11 @@ void PhongRay::on_intersection() normal /= normal.norm(); // Need to apply translations to find the normal vector in - // the base level universe's coordinate system. Uses fact - // that rotation matrices are orthonormal. - auto inverse_rotate = [](const Direction u, - const std::vector& rotation) { - return Direction { - u.x * rotation[0] + u.y * rotation[3] + u.z * rotation[6], - u.x * rotation[1] + u.y * rotation[4] + u.z * rotation[7], - u.x * rotation[2] + u.y * rotation[5] + u.z * rotation[8]}; - }; + // the base level universe's coordinate system. for (int lev = n_coord() - 2; lev >= 0; --lev) { if (coord(lev + 1).rotated) { const Cell& c {*model::cells[coord(lev).cell]}; - normal = inverse_rotate(normal, c.rotation_); + normal = normal.inverse_rotate(c.rotation_); } } diff --git a/src/position.cpp b/src/position.cpp index 5b3613b2897..c9a2aa63677 100644 --- a/src/position.cpp +++ b/src/position.cpp @@ -82,6 +82,13 @@ Position Position::rotate(const vector& rotation) const x * rotation[6] + y * rotation[7] + z * rotation[8]}; } +Position Position::inverse_rotate(const vector& rotation) const +{ + return {x * rotation[0] + y * rotation[3] + z * rotation[6], + x * rotation[1] + y * rotation[4] + z * rotation[7], + x * rotation[2] + y * rotation[5] + z * rotation[8]}; +} + std::ostream& operator<<(std::ostream& os, Position r) { os << "(" << r.x << ", " << r.y << ", " << r.z << ")"; From 0ba097ebcbe0d62c8293c3d6b3c7ee42a8574d08 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 12:18:58 -0600 Subject: [PATCH 07/10] review comments --- include/openmc/plot.h | 3 +- openmc/plots.py | 124 ++++++++++++++++++++---------------------- src/plot.cpp | 2 +- 3 files changed, 63 insertions(+), 66 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 3e1389e1539..9fbd1411f4c 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -288,7 +288,8 @@ class Plot : public PlottableInterface, public SlicePlotBase { * This class serves as a base for plots that create their visuals by tracing * rays from a camera through the problem geometry. It inherits from * PlottableInterface, ensuring that it provides an implementation for - * generating output specific to ray-traced visualization. + * generating output specific to ray-traced visualization. ProjectionPlot + * and PhongPlot provide concrete implementations of this class. */ class RayTracePlot : public PlottableInterface { public: diff --git a/openmc/plots.py b/openmc/plots.py index f19462a6771..320dd6bc2a4 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -170,7 +170,7 @@ def _get_plot_image(plot, cwd): from IPython.display import Image - # Make sure .png file was created +#Make sure.png file was created stem = plot.filename if plot.filename is not None else f'plot_{plot.id}' png_file = Path(cwd) / f'{stem}.png' if not png_file.exists(): @@ -199,14 +199,14 @@ def voxel_to_vtk(voxel_file: PathLike, output: PathLike = 'plot.vti'): Path of the .vti file produced """ - # imported vtk only if used as vtk is an option dependency +#imported vtk only if used as vtk is an option dependency import vtk _min_version = (2, 0) - # Read data from voxel file +#Read data from voxel file with h5py.File(voxel_file, "r") as fh: - # check version +#check version version = tuple(fh.attrs["version"]) if version < _min_version: old_version = ".".join(map(str, version)) @@ -229,15 +229,15 @@ def voxel_to_vtk(voxel_file: PathLike, output: PathLike = 'plot.vti'): grid.SetOrigin(*lower_left) grid.SetSpacing(*width) - # transpose data from OpenMC ordering (zyx) to VTK ordering (xyz) - # and flatten to 1-D array +#transpose data from OpenMC ordering(zyx) to VTK ordering(xyz) +#and flatten to 1 - D array h5data = fh["data"][...] data = vtk.vtkIntArray() data.SetName("id") - # set the array using the h5data array +#set the array using the h5data array data.SetArray(h5data, h5data.size, True) - # add data to image grid +#add data to image grid grid.GetCellData().AddArray(data) writer = vtk.vtkXMLImageDataWriter() @@ -297,7 +297,7 @@ class PlotBase(IDManagerMixin): used_ids = set() def __init__(self, plot_id=None, name=''): - # Initialize Plot class attributes +#Initialize Plot class attributes self.id = plot_id self.name = name self._pixels = [400, 400] @@ -433,8 +433,8 @@ def _check_color(err_string, color): 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 +#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 @@ -458,16 +458,16 @@ def colorize(self, geometry, seed=1): cv.check_type('seed', seed, Integral) cv.check_greater_than('seed', seed, 1, equality=True) - # Get collections of the domains which will be plotted +#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 +#Set the seed for the random number generator np.random.seed(seed) - # Generate random colors for each feature +#Generate random colors for each feature for domain in domains: self.colors[domain] = np.random.randint(0, 256, (3,)) @@ -705,7 +705,7 @@ def from_geometry(cls, geometry, basis='xy', slice_coord=0.): cv.check_type('geometry', geometry, openmc.Geometry) cv.check_value('basis', basis, _BASES) - # Decide which axes to keep +#Decide which axes to keep if basis == 'xy': pick_index = (0, 1) slice_index = 2 @@ -716,7 +716,7 @@ def from_geometry(cls, geometry, basis='xy', slice_coord=0.): pick_index = (0, 2) slice_index = 1 - # Get lower-left and upper-right coordinates for desired axes +#Get lower - left and upper - right coordinates for desired axes lower_left, upper_right = geometry.bounding_box lower_left = lower_left[np.array(pick_index)] upper_right = upper_right[np.array(pick_index)] @@ -762,17 +762,17 @@ def highlight_domains(self, geometry, domains, seed=1, cv.check_less_than('alpha', alpha, 1., equality=True) cv.check_type('background', background, Iterable) - # Get a background (R,G,B) tuple to apply in alpha compositing +#Get a background(R, G, B) tuple to apply in alpha compositing if isinstance(background, str): if background.lower() not in _SVG_COLORS: raise ValueError(f"'{background}' is not a valid color.") background = _SVG_COLORS[background.lower()] - # Generate a color scheme +#Generate a color scheme self.colorize(geometry, seed) - # Apply alpha compositing to the colors for all domains - # other than those the user wishes to highlight +#Apply alpha compositing to the colors for all domains +#other than those the user wishes to highlight for domain, color in self.colors.items(): if domain not in domains: if isinstance(color, str): @@ -862,7 +862,7 @@ def from_xml_element(cls, elem): plot.pixels = get_elem_tuple(elem, "pixels") plot._background = get_elem_tuple(elem, "background") - # Set plot colors +#Set plot colors colors = {} for color_elem in elem.findall("color"): uid = int(color_elem.get("id")) @@ -870,7 +870,7 @@ def from_xml_element(cls, elem): for x in color_elem.get("rgb").split()]) plot.colors = colors - # Set masking information +#Set masking information mask_elem = elem.find("mask") if mask_elem is not None: plot.mask_components = [ @@ -880,7 +880,7 @@ def from_xml_element(cls, elem): plot.mask_background = tuple( [int(x) for x in background.split()]) - # show overlaps +#show overlaps overlap_elem = elem.find("show_overlaps") if overlap_elem is not None: plot.show_overlaps = (overlap_elem.text in ('true', '1')) @@ -888,12 +888,12 @@ def from_xml_element(cls, elem): if overlap_color is not None: plot.overlap_color = overlap_color - # Set universe level +#Set universe level level = elem.find("level") if level is not None: plot.level = int(level.text) - # Set meshlines +#Set meshlines mesh_elem = elem.find("meshlines") if mesh_elem is not None: meshlines = {'type': mesh_elem.get('meshtype')} @@ -931,13 +931,13 @@ def to_ipython_image(self, openmc_exec='openmc', cwd='.'): Image generated """ - # Create plots.xml +#Create plots.xml Plots([self]).export_to_xml(cwd) - # Run OpenMC in geometry plotting mode +#Run OpenMC in geometry plotting mode openmc.plot_geometry(False, openmc_exec, cwd) - # Return produced image +#Return produced image return _get_plot_image(self, cwd) def to_vtk(self, output: Optional[PathLike] = None, @@ -967,10 +967,10 @@ def to_vtk(self, output: Optional[PathLike] = None, raise ValueError( 'Generating a VTK file only works for voxel plots') - # Create plots.xml +#Create plots.xml Plots([self]).export_to_xml(cwd) - # Run OpenMC in geometry plotting mode and produces a h5 file +#Run OpenMC in geometry plotting mode and produces a h5 file openmc.plot_geometry(False, openmc_exec, cwd) stem = self.filename if self.filename is not None else f'plot_{self.id}' @@ -1022,7 +1022,7 @@ class RayTracePlot(PlotBase): """ def __init__(self, plot_id=None, name=''): - # Initialize Plot class attributes +#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) @@ -1087,9 +1087,8 @@ def _check_domains_consistent_with_color_by(self, domains): """ for region in domains: - - # if an integer is passed, we have to assume - # it was a valid ID +#if an integer is passed, we have to assume +#it was a valid ID if isinstance(region, int): continue @@ -1124,7 +1123,7 @@ def to_xml_element(self): 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 +#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) @@ -1177,18 +1176,18 @@ def _read_xml_attributes(self, elem): self.camera_position = get_elem_tuple(elem, "camera_position", float) self.look_at = get_elem_tuple(elem, "look_at", float) - # Set masking information +#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) +#TODO : set mask components(needs geometry information) background = mask_elem.get("background") if background is not None: self.mask_background = tuple( [int(x) for x in background.split()]) - # Set universe level +#Set universe level level = elem.find("level") if level is not None: self.level = int(level.text) @@ -1217,7 +1216,7 @@ class ProjectionPlot(RayTracePlot): 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. + values can be obtained using the :meth:`set_transparent` method. """ def __init__(self, plot_id=None, name=''): @@ -1279,13 +1278,13 @@ def set_transparent(self, geometry): cv.check_type('geometry', geometry, openmc.Geometry) - # Get collections of the domains which will be plotted +#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 +#Generate random colors for each feature for domain in domains: self.xs[domain] = 0.0 @@ -1322,15 +1321,15 @@ def to_xml_element(self): color = _SVG_COLORS[color.lower()] subelement.text = ' '.join(str(x) for x in color) - self._check_domains_consistent_with_color_by(self._wireframe_domains) + self._check_domains_consistent_with_color_by(self.wireframe_domains) 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 +#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): @@ -1365,7 +1364,7 @@ def from_xml_element(cls, elem): plot._read_xml_attributes(elem) - # Attempt to get wireframe thickness. May not be present +#Attempt to get wireframe thickness.May not be present wireframe_thickness = elem.get("wireframe_thickness") if wireframe_thickness: plot.wireframe_thickness = int(wireframe_thickness) @@ -1373,13 +1372,11 @@ def from_xml_element(cls, elem): if wireframe_color: plot.wireframe_color = [int(item) for item in wireframe_color] - # Set plot colors - plot._colors = {} - plot._xs = {} +#Set plot colors for color_elem in elem.findall("color"): uid = color_elem.get("id") - plot._colors[uid] = get_elem_tuple(color_elem, "rgb") - plot._xs[uid] = float(color_elem.get("xs")) + plot.colors[uid] = get_elem_tuple(color_elem, "rgb") + plot.xs[uid] = float(color_elem.get("xs")) return plot @@ -1445,21 +1442,21 @@ def to_xml_element(self): element = super().to_xml_element() element.set("type", "phong") - # no light position means put it at the camera +#no light position means put it at the camera if self._light_position: subelement = ET.SubElement(element, "light_position") subelement.text = ' '.join(map(str, self._light_position)) - # no diffuse fraction defaults to 0.1 +#no diffuse fraction defaults to 0.1 if self._diffuse_fraction: subelement = ET.SubElement(element, "diffuse_fraction") subelement.text = str(self._diffuse_fraction) - self._check_domains_consistent_with_color_by(self._opaque_domains) + self._check_domains_consistent_with_color_by(self.opaque_domains) subelement = ET.SubElement(element, "opaque_ids") - # Extract all IDs, or use the integer value passed in - # explicitly if that was given +#Extract all IDs, or use the integer value passed in +#explicitly if that was given subelement.text = ' '.join( [str(domain) if isinstance(domain, int) else str(domain.id) for domain in self._opaque_domains]) @@ -1491,11 +1488,10 @@ def from_xml_element(cls, elem): plot._read_xml_attributes(elem) - # Set plot colors - plot._colors = {} +#Set plot colors for color_elem in elem.findall("color"): uid = color_elem.get("id") - plot._colors[uid] = get_elem_tuple(color_elem, "rgb") + plot.colors[uid] = get_elem_tuple(color_elem, "rgb") class Plots(cv.CheckedList): @@ -1615,14 +1611,14 @@ def to_xml_element(self): XML element containing all plot elements """ - # Reset xml element tree +#Reset xml element tree self._plots_file.clear() self._create_plot_subelements() - # Clean the indentation in the file to be user-readable +#Clean the indentation in the file to be user - readable clean_indentation(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 @@ -1636,13 +1632,13 @@ def export_to_xml(self, path='plots.xml'): Path to file to write. Defaults to 'plots.xml'. """ - # Check if path is a directory +#Check if path is a directory p = Path(path) if p.is_dir(): p /= 'plots.xml' self.to_xml_element() - # Write the XML Tree to the plots.xml file +#Write the XML Tree to the plots.xml file tree = ET.ElementTree(self._plots_file) tree.write(str(p), xml_declaration=True, encoding='utf-8') @@ -1661,7 +1657,7 @@ def from_xml_element(cls, elem): Plots collection """ - # Generate each plot +#Generate each plot plots = cls() for e in elem.findall('plot'): plot_type = e.get('type') diff --git a/src/plot.cpp b/src/plot.cpp index f0d8d3fae3e..7989afbee6c 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1770,7 +1770,7 @@ void PhongRay::on_intersection() // Facing away from the light means no lighting double dotprod = normal.dot(to_light); - dotprod = dotprod >= 0.0 ? dotprod : 0.0; + dotprod = std::max(0.0, dotprod); double modulation = plot_.diffuse_fraction_ + (1.0 - plot_.diffuse_fraction_) * dotprod; From f064a81bc0e3294bb07e2ea8ba624503baabaadd Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 12:25:46 -0600 Subject: [PATCH 08/10] get clang-format to pass --- include/openmc/plot.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 9fbd1411f4c..9ccaf75d81a 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -245,8 +245,8 @@ T SlicePlotBase::get_map() const data.set_overlap(y, x); } } // inner for - } // outer for - } // omp parallel + } + } return data; } From 30140fd38a86f9263e84602104c96caa36f7bce3 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Mar 2024 12:29:25 -0600 Subject: [PATCH 09/10] try to make clangformat happy --- src/plot.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/plot.cpp b/src/plot.cpp index 7989afbee6c..8dc66107642 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -264,8 +264,8 @@ void Plot::create_image() const } data(x, y) = colors_[model::material_map[id]]; } // color_by if-else - } // x for loop - } // y for loop + } + } // draw mesh lines if present if (index_meshlines_mesh_ >= 0) { From 97435eb74953ac5a0a60ade6d8657969b5653287 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Sat, 16 Mar 2024 14:05:01 -0500 Subject: [PATCH 10/10] update comment --- src/plot.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/plot.cpp b/src/plot.cpp index 8dc66107642..338c55213b6 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1195,9 +1195,7 @@ bool ProjectionPlot::trackstack_equivalent( std::pair RayTracePlot::get_pixel_ray( int horiz, int vert) const { - // 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. + // Compute field of view in radians 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]);