diff --git a/docs/source/_images/phong_triso.png b/docs/source/_images/phong_triso.png
new file mode 100644
index 00000000000..7f1ac61cab6
Binary files /dev/null and b/docs/source/_images/phong_triso.png differ
diff --git a/docs/source/io_formats/plots.rst b/docs/source/io_formats/plots.rst
index c1fa7833083..1a42a4281d7 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 ``<plots>`` and any number output plots can be
-defined with ``<plot>`` sub-elements. Two plot types are currently implemented
+defined with ``<plot>`` 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.
+* ``wireframe_raytrace`` 2D pixel plot of a three-dimensional view of a
+ geometry using wireframes around cells or materials and coloring by depth
+ through each material.
+* ``solid_raytrace`` 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"
- 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 :func:`openmc.voxel_to_vtk` function and
- subsequently viewed with a 3D viewer such as VISIT or Paraview. See the
- :ref:`io_voxel` for information about the datafile structure.
+ Keyword for type of plot to be produced. Currently "slice", "voxel",
+ "wireframe_raytrace", and "solid_raytrace" 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 :func:`openmc.voxel_to_vtk` function and subsequently viewed with a 3D
+ viewer such as VISIT or Paraview. 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"
-``<plot>`` elements of ``type`` "slice" and "voxel" must contain the ``pixels``
+All ``<plot>`` 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
``<plot>`` 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:
+ The attenuation coefficient for volume rendering of color in units of
+ inverse centimeters. Zero corresponds to transparency. Only for plot type
+ "wireframe_raytrace".
+
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:
@@ -191,3 +202,80 @@ attributes or sub-elements. These are not used in "voxel" plots:
*Default*: 0 0 0 (black)
*Default*: None
+
+``<plot>`` elements of ``type`` "wireframe_raytrace" or "solid_raytrace" can contain the
+following attributes or sub-elements.
+
+ :camera_position:
+ Location in 3D Cartesian space the camera is at.
+
+
+ *Default*: None - Required for all ``wireframe_raytrace`` or
+ ``solid_raytrace`` plots
+
+ :look_at:
+ Location in 3D Cartesian space the camera is looking at.
+
+
+ *Default*: None - Required for all ``wireframe_raytrace`` or
+ ``solid_raytrace`` 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
+
+``<plot>`` elements of ``type`` "solid_raytrace" 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
+
+``<plot>`` elements of ``type`` "wireframe_raytrace" 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/pythonapi/base.rst b/docs/source/pythonapi/base.rst
index 23df02f2ee7..d6d04444f2a 100644
--- a/docs/source/pythonapi/base.rst
+++ b/docs/source/pythonapi/base.rst
@@ -165,7 +165,8 @@ Geometry Plotting
:template: myclass.rst
openmc.Plot
- openmc.ProjectionPlot
+ openmc.WireframeRayTracePlot
+ openmc.SolidRayTracePlot
openmc.Plots
Running OpenMC
diff --git a/docs/source/releasenotes/0.14.0.rst b/docs/source/releasenotes/0.14.0.rst
index 445edfa077a..f82aa30850d 100644
--- a/docs/source/releasenotes/0.14.0.rst
+++ b/docs/source/releasenotes/0.14.0.rst
@@ -52,7 +52,7 @@ Compatibility Notes and Deprecations
New Features
------------
-- A new :class:`openmc.ProjectionPlot` class enables the generation of orthographic or
+- A new :class:`openmc.WireframeRayTracePlot` class enables the generation of orthographic or
perspective projection plots. (`#1926
<https://github.com/openmc-dev/openmc/pull/1926>`_)
- The :class:`openmc.model.RightCircularCylinder` class now supports optional
diff --git a/docs/source/usersguide/plots.rst b/docs/source/usersguide/plots.rst
index d453ea4537a..da0c69bdd8d 100644
--- a/docs/source/usersguide/plots.rst
+++ b/docs/source/usersguide/plots.rst
@@ -122,27 +122,82 @@ for doing this 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.
-----------------
-Projection Plots
-----------------
+----------------------
+Solid Ray-traced Plots
+----------------------
+
+.. image:: ../_images/phong_triso.png
+ :width: 300px
+
+The :class:`openmc.SolidRayTracePlot` 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
+<https://en.wikipedia.org/wiki/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.
+
+Solid ray-traced 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 solid ray-traced plot. That being said, these plots are
+not useful for detecting overlap or undefined regions, so it is recommended to
+use the slice plot approach for geometry debugging.
+
+Only a few inputs are required for a solid ray-traced plot. The camera location,
+where the camera is 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 that are opaque in the plot must correspond to
+materials if coloring by material, or cells if coloring by cell.
+
+A minimal solid ray-traced plot input could be::
+
+ plot = openmc.SolidRayTracePlot()
+ 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.
+
+---------------
+Wireframe Plots
+---------------
.. only:: html
.. image:: ../_images/hexlat_anim.gif
:width: 200px
-The :class:`openmc.ProjectionPlot` class presents an alternative method of
-producing 3D visualizations of OpenMC geometries. It was developed to overcome
-the primary shortcoming of voxel plots, that an enormous number of voxels must
-be employed to capture detailed geometric features. Projection plots perform
-volume rendering on material or cell volumes, with colors specified in the same
-manner as slice plots. This is done using the native ray tracing capabilities
+The :class:`openmc.WireframeRayTracePlot` 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 the 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 solid ray-traced 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.
+overlaps or leaks will work with wireframe plots.
-One drawback of projection plots is that particle tracks cannot be overlaid on
+One drawback of wireframe plots is that particle tracks cannot be overlaid on
them at present. Moreover, checking for overlap regions is not currently
-possible with projection plots. The image heading this section can be created by
+possible with wireframe plots. The image heading this section can be created by
adding the following code to the hexagonal lattice example packaged with OpenMC,
before exporting to plots.xml.
@@ -152,7 +207,7 @@ before exporting to plots.xml.
import numpy as np
for i in range(100):
phi = 2 * np.pi * i/100
- thisp = openmc.ProjectionPlot(plot_id = 4 + i)
+ thisp = openmc.WireframeRayTracePlot(plot_id = 4 + i)
thisp.filename = 'frame%s'%(str(i).zfill(3))
thisp.look_at = [0, 0, 0]
thisp.camera_position = [r * np.cos(phi), r * np.sin(phi), 6 * np.sin(phi)]
@@ -167,42 +222,45 @@ before exporting to plots.xml.
plot_file.append(thisp)
-This generates a sequence of png files which can be joined to form a gif. Each
+This generates a sequence of png files that can be joined to form a gif. Each
image specifies a different camera position using some simple periodic functions
-to create a perfectly looped gif. :attr:`ProjectionPlot.look_at` defines where
-the camera's centerline should point at. :attr:`ProjectionPlot.camera_position`
-similarly defines where the camera is situated in the universe level we seek to
-plot. The other settings resemble those employed by :class:`openmc.Plot`, with
-the exception of the :class:`ProjectionPlot.set_transparent` method and
-:attr:`ProjectionPlot.xs` dictionary. These are used to control volume rendering
-of material volumes. "xs" here stands for cross section, and it defines material
-opacities in units of inverse centimeters. Setting this value to a large number
-would make a material or cell opaque, and setting it to zero makes a material
-transparent. Thus, the :class:`ProjectionPlot.set_transparent` can be used to
-make all materials in the geometry transparent. From there, individual material
-or cell opacities can be tuned to produce the desired result.
+to create a perfectly looped gif. :attr:`~WireframeRayTracePlot.look_at` defines
+where the camera's centerline should point at.
+:attr:`~WireframeRayTracePlot.camera_position` similarly defines where the
+camera is situated in the universe level we seek to plot. The other settings
+resemble those employed by :class:`openmc.Plot`, with the exception of the
+:meth:`~WireframeRayTracePlot.set_transparent` method and
+:attr:`~WireframeRayTracePlot.xs` dictionary. These are used to control volume
+rendering of material volumes. "xs" here stands for cross section, and it
+defines material opacities in units of inverse centimeters. Setting this value
+to a large number would make a material or cell opaque, and setting it to zero
+makes a material transparent. Thus, the
+:meth:`~WireframeRayTracePlot.set_transparent` method can be used to make all
+materials in the geometry transparent. From there, individual material or cell
+opacities can be tuned to produce the desired result.
Two camera projections are available when using these plots, perspective and
orthographic. The default, perspective projection, is a cone of rays passing
through each pixel which radiate from the camera position and span the field of
view in the x and y positions. The horizontal field of view can be set with the
-:attr: `ProjectionPlot.horizontal_field_of_view` attribute, which is to be
-specified in units of degrees. The field of view only influences behavior in
+:attr:`~WireframeRayTracePlot.horizontal_field_of_view` attribute, which is to
+be specified in units of degrees. The field of view only influences behavior in
perspective projection mode.
In the orthographic projection, rays follow the same angle but originate from
different points. The horizontal width of this plane of ray starting points may
-be set with the :attr: `ProjectionPlot.orthographic_width` element. If this
-element is nonzero, the orthographic projection is employed. Left to its default
-value of zero, the perspective projection is employed.
-
-Lastly, projection plots come packaged with wireframe generation that can target
-either all surface/cell/material boundaries in the geometry, or only wireframing
-around specific regions. In the above example, we have set only the fuel region
-from the hexagonal lattice example to have a wireframe drawn around it. This is
-accomplished by setting the :attr: `ProjectionPlot.wireframe_domains`, which may
-be set to either material IDs or cell IDs. The
-:attr:`ProjectionPlot.wireframe_thickness` attribute sets the wireframe
+be set with the :attr:`~WireframeRayTracePlot.orthographic_width` attribute. If
+this element is nonzero, the orthographic projection is employed. Left to its
+default value of zero, the perspective projection is employed.
+
+Most importantly, wireframe plots come packaged with wireframe generation that
+can target either all surface/cell/material boundaries in the geometry, or only
+wireframing around specific regions. In the above example, we have set only the
+fuel region from the hexagonal lattice example to have a wireframe drawn around
+it. This is accomplished by setting the
+:attr:`~WireframeRayTracePlot.wireframe_domains` attribute, which may be set to
+either material IDs or cell IDs. The
+:attr:`~WireframeRayTracePlot.wireframe_thickness` attribute sets the wireframe
thickness in units of pixels.
.. note:: When setting specific material or cell regions to have wireframes
diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index d01020f8e7a..d7a3cad182b 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/dagmc.h b/include/openmc/dagmc.h
index 47fcfe237e3..9c4e47fdfd3 100644
--- a/include/openmc/dagmc.h
+++ b/include/openmc/dagmc.h
@@ -49,7 +49,8 @@ class DAGSurface : public Surface {
double evaluate(Position r) const override;
double distance(Position r, Direction u, bool coincident) const override;
Direction normal(Position r) const override;
- Direction reflect(Position r, Direction u, GeometryState* p) const override;
+ Direction reflect(
+ Position r, Direction u, GeometryState* p = nullptr) const override;
inline void to_hdf5_inner(hid_t group_id) const override {};
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index b1f9fabd110..aaac864e11e 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -39,10 +39,6 @@ class Particle : public ParticleData {
double speed() const;
- //! moves the particle by the distance length to its next location
- //! \param length the distance the particle is moved
- void move_distance(double length);
-
//! create a secondary particle
//
//! stores the current phase space attributes of the particle in the
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 2255d1a7d88..e2fac79e372 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -191,6 +191,13 @@ struct BoundaryInfo {
array<int, 3>
lattice_translation {}; //!< which way lattice indices will change
+ void reset()
+ {
+ distance = INFINITY;
+ surface = SURFACE_NONE;
+ coord_level = 0;
+ lattice_translation = {0, 0, 0};
+ }
// TODO: off-by-one
int surface_index() const { return std::abs(surface) - 1; }
};
@@ -226,6 +233,12 @@ class GeometryState {
n_coord_last_ = 1;
}
+ //! moves the particle by the specified distance to its next location
+ //! \param distance the distance the particle is moved
+ void move_distance(double distance);
+
+ void advance_to_boundary_from_void();
+
// Initialize all internal state from position and direction
void init_from_r_u(Position r_a, Direction u_a)
{
@@ -565,7 +578,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 ee6c3fbec58..69e801c5fc3 100644
--- a/include/openmc/plot.h
+++ b/include/openmc/plot.h
@@ -4,6 +4,7 @@
#include <cmath>
#include <sstream>
#include <unordered_map>
+#include <unordered_set>
#include "pugixml.hpp"
#include "xtensor/xarray.hpp"
@@ -61,6 +62,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;
};
@@ -70,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:
@@ -232,8 +246,8 @@ T SlicePlotBase::get_map() const
data.set_overlap(y, x);
}
} // inner for
- } // outer for
- } // omp parallel
+ }
+ }
return data;
}
@@ -268,32 +282,110 @@ class Plot : public PlottableInterface, public SlicePlotBase {
RGBColor meshlines_color_; //!< Color of meshlines on the plot
};
-class ProjectionPlot : public PlottableInterface {
-
+/**
+ * \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. WireframeRayTracePlot
+ * and SolidRayTracePlot provide concrete implementations of this class.
+ */
+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:
+protected:
+ Direction camera_x_axis() const
+ {
+ return {camera_to_model_[0], camera_to_model_[3], camera_to_model_[6]};
+ }
+
+ Direction camera_y_axis() const
+ {
+ return {camera_to_model_[1], camera_to_model_[4], camera_to_model_[7]};
+ }
+
+ Direction camera_z_axis() const
+ {
+ return {camera_to_model_[2], camera_to_model_[5], camera_to_model_[8]};
+ }
+
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<Position, Direction> get_pixel_ray(int horiz, int vert) const;
+
+ std::array<int, 2> 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
+
+ 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_ {C_NONE};
+
+ /*
+ * Cached camera-to-model matrix with column vectors of axes. The x-axis is
+ * the vector between the camera_position_ and look_at_; the y-axis is the
+ * cross product of the x-axis with the up_ vector, and the z-axis is the
+ * cross product of the x and y axes.
+ */
+ std::array<double, 9> camera_to_model_;
+};
+
+class ProjectionRay;
+
+/**
+ * \class WireframeRayTracePlot
+ * \brief Creates plots that are like colorful x-ray imaging
+ *
+ * WireframeRayTracePlot 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 WireframeRayTracePlot : public RayTracePlot {
+
+ friend class ProjectionRay;
+
+public:
+ WireframeRayTracePlot(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.
@@ -314,35 +406,138 @@ class ProjectionPlot : public PlottableInterface {
* if two surfaces bound a single cell, it allows drawing that sharp edge
* where the surfaces intersect.
*/
- int surface; // last surface ID intersected in this segment
+ int surface_index {-1}; // last surface index intersected in this segment
TrackSegment(int id_a, double length_a, int surface_a)
- : id(id_a), length(length_a), surface(surface_a)
+ : id(id_a), length(length_a), surface_index(surface_a)
{}
};
+ // which color IDs should be wireframed. If empty, all cells are wireframed.
+ vector<int> wireframe_ids_;
+
+ // Thickness of the wireframe lines. Can set to zero for no wireframe.
+ int wireframe_thickness_ {1};
+
+ RGBColor wireframe_color_ {BLACK}; // wireframe color
+ vector<double> xs_; // macro cross section values for cell volume rendering
+};
+
+/**
+ * \class SolidRayTracePlot
+ * \brief Plots 3D objects as the eye might see them.
+ *
+ * Plots a geometry with single-scattered Phong lighting plus a diffuse lighting
+ * contribution. The result is a physically reasonable, aesthetic 3D view of a
+ * geometry.
+ */
+class SolidRayTracePlot : public RayTracePlot {
+ friend class PhongRay;
+
+public:
+ SolidRayTracePlot(pugi::xml_node plot);
+
+ virtual void create_output() const;
+ virtual void print_info() 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);
+
+ std::unordered_set<int> 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();
+
+ // Stops the ray and exits tracing when called from on_intersection
+ void stop() { stop_ = true; }
+
+ // Sets the dist_ variable
+ void compute_distance();
+
+protected:
+ // Records how far the ray has traveled
+ double traversal_distance_ {0.0};
+
+private:
// Max intersections before we assume ray tracing is caught in an infinite
// loop:
static const int MAX_INTERSECTIONS = 1000000;
- std::array<int, 2> 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
+ bool hit_something_ {false};
+ bool stop_ {false};
- // which color IDs should be wireframed. If empty, all cells are wireframed.
- vector<int> wireframe_ids_;
+ unsigned event_counter_ {0};
+};
- /* The horizontal thickness, if using an orthographic projection.
- * If set to zero, we assume using a perspective projection.
+class ProjectionRay : public Ray {
+public:
+ ProjectionRay(Position r, Direction u, const WireframeRayTracePlot& plot,
+ vector<WireframeRayTracePlot::TrackSegment>& 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.
*/
- double orthographic_width_ {0.0};
+ const WireframeRayTracePlot& plot_;
- // Thickness of the wireframe lines. Can set to zero for no wireframe.
- int wireframe_thickness_ {1};
+ /* The ray runs through the geometry, and records the lengths of ray segments
+ * and cells they lie in along the way.
+ */
+ vector<WireframeRayTracePlot::TrackSegment>& line_segments_;
+};
- RGBColor wireframe_color_ {BLACK}; // wireframe color
- vector<double> xs_; // macro cross section values for cell volume rendering
+class PhongRay : public Ray {
+public:
+ PhongRay(Position r, Direction u, const SolidRayTracePlot& 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 SolidRayTracePlot& 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_;
};
//===============================================================================
diff --git a/include/openmc/position.h b/include/openmc/position.h
index dc60ab35a04..5d291d26b95 100644
--- a/include/openmc/position.h
+++ b/include/openmc/position.h
@@ -94,8 +94,24 @@ struct Position {
//! \result Reflected vector
Position reflect(Position n) const;
- //! Rotate the position based on a rotation matrix
- Position rotate(const vector<double>& rotation) const;
+ //! Rotate the position by applying a rotation matrix
+ template<typename T>
+ Position rotate(const T& rotation) const
+ {
+ return {x * rotation[0] + y * rotation[1] + z * rotation[2],
+ x * rotation[3] + y * rotation[4] + z * rotation[5],
+ x * rotation[6] + y * rotation[7] + z * rotation[8]};
+ }
+
+ //! Rotate the position by applying the inverse of a rotation matrix
+ //! using the fact that rotation matrices are orthonormal.
+ template<typename T>
+ Position inverse_rotate(const T& 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]};
+ }
// Data members
double x = 0.;
diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 498f71d4f9b..2397a64cc4c 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -62,7 +62,7 @@ class Surface {
Position r, Direction u, GeometryState* p = nullptr) const;
virtual Direction diffuse_reflect(
- Position r, Direction u, uint64_t* seed, GeometryState* p = nullptr) const;
+ Position r, Direction u, uint64_t* seed) const;
//! Evaluate the equation describing the surface.
//!
diff --git a/openmc/_xml.py b/openmc/_xml.py
index b40ecb258a9..17389a82f0c 100644
--- a/openmc/_xml.py
+++ b/openmc/_xml.py
@@ -82,7 +82,7 @@ def reorder_attributes(root):
def get_elem_tuple(elem, name, dtype=int):
- '''Helper function to get a tuple of values from an elem
+ """Helper function to get a tuple of values from an elem
Parameters
----------
@@ -97,7 +97,7 @@ def get_elem_tuple(elem, name, dtype=int):
-------
tuple of dtype
Data read from the tuple
- '''
+ """
subelem = elem.find(name)
if subelem is not None:
return tuple([dtype(x) for x in subelem.text.split()])
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 8dd13ef6b3c..5e93d98040c 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -124,7 +124,7 @@ def plots(self) -> openmc.Plots | None:
@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:
@@ -220,7 +220,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)
@@ -242,12 +243,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'))
@@ -538,11 +543,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("</model>\n")
diff --git a/openmc/plots.py b/openmc/plots.py
index a59a7f0eb4a..24607659844 100644
--- a/openmc/plots.py
+++ b/openmc/plots.py
@@ -555,13 +555,21 @@ def colorize(self, geometry, seed=1):
else:
domains = geometry.get_all_cells().values()
- # Set the seed for the random number generator
rng = np.random.RandomState(seed)
# Generate random colors for each feature
for domain in domains:
self.colors[domain] = rng.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
@@ -887,13 +895,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")
@@ -1051,7 +1053,8 @@ def to_vtk(self, output: PathLike | None = 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)
@@ -1072,20 +1075,17 @@ def to_vtk(self, output: PathLike | None = 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.
+ projections are more similar to a pinhole camera, and orthographic
+ projections preserve parallel lines and distances.
- .. versionadded:: 0.14.0
+ This is an abstract base class that :class:`WireframeRayTracePlot` and
+ :class:`SolidRayTracePlot` finish the implementation of.
+
+ .. versionadded:: 0.15.1
Parameters
----------
@@ -1113,21 +1113,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=''):
@@ -1138,10 +1123,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):
@@ -1195,6 +1176,161 @@ 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 cells 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)
+
+ if (tmp := elem.find("orthographic_width")) 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)
+
+ if elem.find("background") is not None:
+ self.background = get_elem_tuple(elem, "background")
+
+ # Set masking information
+ if (mask_elem := elem.find("mask")) 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 WireframeRayTracePlot(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.
+
+ .. versionchanged:: 0.15.1
+ Renamed from ProjectionPlot to WireframeRayTracePlot
+
+ Parameters
+ ----------
+ plot_id : int
+ Unique identifier for the plot
+ name : str
+ Name of the plot
+
+ Attributes
+ ----------
+ 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 :meth:`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
@@ -1221,15 +1357,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
@@ -1266,6 +1393,18 @@ def set_transparent(self, geometry):
for domain in domains:
self.xs[domain] = 0.0
+ def __repr__(self):
+ string = 'Wireframe Ray-traced 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
@@ -1275,15 +1414,8 @@ 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))
+ element.set("type", "wireframe_raytrace")
subelement = ET.SubElement(element, "wireframe_thickness")
subelement.text = str(self._wireframe_thickness)
@@ -1294,6 +1426,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")
@@ -1311,43 +1445,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
@@ -1359,60 +1458,193 @@ def from_xml_element(cls, elem):
Returns
-------
- openmc.ProjectionPlot
- ProjectionPlot object
+ openmc.WireframeRayTracePlot
+ WireframeRayTracePlot 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_id = int(elem.get("id"))
+ plot_name = get_text(elem, 'name', '')
+ plot = cls(plot_id, plot_name)
+ plot.type = "wireframe_raytrace"
- 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")
- if wireframe_thickness:
- plot.wireframe_thickness = int(wireframe_thickness)
+ # Attempt to get wireframe thickness.May not be present
+ wireframe_thickness = elem.find("wireframe_thickness")
+ if wireframe_thickness is not None:
+ plot.wireframe_thickness = int(wireframe_thickness.text)
wireframe_color = elem.get("wireframe_color")
if wireframe_color:
plot.wireframe_color = [int(item) for item in wireframe_color]
# Set plot colors
- colors = {}
- xs = {}
for color_elem in elem.findall("color"):
- uid = color_elem.get("id")
- colors[uid] = get_elem_tuple(color_elem, "rgb")
- xs[uid] = float(color_elem.get("xs"))
+ uid = int(color_elem.get("id"))
+ plot.colors[uid] = tuple(int(i)
+ for i in get_text(color_elem, 'rgb').split())
+ 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)
+
+class SolidRayTracePlot(RayTracePlot):
+ """Phong shading-based rendering of an OpenMC geometry
+
+ This class defines a plot that uses Phong shading to enhance the
+ visualization of an OpenMC geometry by incorporating diffuse lighting and
+ configurable opacity for certain regions. It extends :class:`RayTracePlot`
+ by adding parameters related to lighting and transparency.
+
+ .. versionadded:: 0.15.1
+
+ Parameters
+ ----------
+ plot_id : int, optional
+ Unique identifier for the plot
+ name : str, optional
+ Name of the plot
+
+ Attributes
+ ----------
+ light_position : tuple or list of float
+ Position of the light source in 3D space. Defaults to None, which places
+ the light at the camera position.
+ diffuse_fraction : float
+ Fraction of lighting that is diffuse (non-directional). Defaults to 0.1.
+ Must be between 0 and 1.
+ opaque_domains : list
+ List of domains (e.g., cells or materials) that should be rendered as
+ opaque rather than allowing transparency.
+ """
+
+ 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):
+ # Note that _check_domains_consistent_with_color_by checks
+ # the types within later. This is because we don't necessarily
+ # know what types are acceptable until the user has set the
+ # color_by attribute, too.
+ cv.check_type('opaque domains', x, Iterable)
+ self._opaque_domains = x
+
+ def __repr__(self):
+ string = 'Solid Ray-traced 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 solid ray-traced plot
+
+ Returns
+ -------
+ element : lxml.etree._Element
+ XML element containing plot data
+
+ """
+ element = super().to_xml_element()
+ element.set("type", "solid_raytrace")
+
+ # 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
+
+ def _read_phong_attributes(self, elem):
+ """Read attributes specific to the Phong plot from an XML element"""
+ if elem.find('light_position') is not None:
+ self.light_position = get_elem_tuple(elem, 'light_position', float)
+
+ diffuse_fraction = elem.find('diffuse_fraction')
+ if diffuse_fraction is not None:
+ self.diffuse_fraction = float(diffuse_fraction.text)
+
+ if elem.find('opaque_ids') is not None:
+ self.opaque_domains = list(get_elem_tuple(elem, 'opaque_ids', int))
+
+ @classmethod
+ def from_xml_element(cls, elem):
+ """Generate plot object from an XML element
+
+ Parameters
+ ----------
+ elem : lxml.etree._Element
+ XML element
+
+ Returns
+ -------
+ openmc.WireframeRayTracePlot
+ WireframeRayTracePlot object
+
+ """
+
+ plot_id = int(elem.get("id"))
+ plot_name = get_text(elem, 'name', '')
+ plot = cls(plot_id, plot_name)
+ plot.type = "solid_raytrace"
+
+ plot._read_xml_attributes(elem)
+ plot._read_phong_attributes(elem)
+
+ # Set plot colors
+ for color_elem in elem.findall("color"):
+ uid = color_elem.get("id")
+ plot.colors[uid] = get_elem_tuple(color_elem, "rgb")
return plot
@@ -1434,13 +1666,13 @@ class Plots(cv.CheckedList):
Parameters
----------
- plots : Iterable of openmc.Plot or openmc.ProjectionPlot
+ plots : Iterable of openmc.PlotBase
plots to add to the collection
"""
def __init__(self, plots=None):
- super().__init__((Plot, ProjectionPlot), 'plots collection')
+ super().__init__(PlotBase, 'plots collection')
self._plots_file = ET.Element("plots")
if plots is not None:
self += plots
@@ -1450,7 +1682,7 @@ def append(self, plot):
Parameters
----------
- plot : openmc.Plot or openmc.ProjectionPlot
+ plot : openmc.PlotBase
Plot to append
"""
@@ -1582,10 +1814,14 @@ def from_xml_element(cls, elem):
plots = cls()
for e in elem.findall('plot'):
plot_type = e.get('type')
- if plot_type == 'projection':
- plots.append(ProjectionPlot.from_xml_element(e))
- else:
+ if plot_type == 'wireframe_raytrace':
+ plots.append(WireframeRayTracePlot.from_xml_element(e))
+ elif plot_type == 'solid_raytrace':
+ plots.append(SolidRayTracePlot.from_xml_element(e))
+ elif plot_type in ('slice', 'voxel'):
plots.append(Plot.from_xml_element(e))
+ else:
+ raise ValueError("Unknown plot type: {}".format(plot_type))
return plots
@classmethod
diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index 52c45f21c28..e24b3fbf629 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -808,15 +808,12 @@ Direction DAGSurface::normal(Position r) const
Direction DAGSurface::reflect(Position r, Direction u, GeometryState* p) const
{
Expects(p);
- p->history().reset_to_last_intersection();
- moab::ErrorCode rval;
- moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
double pnt[3] = {r.x, r.y, r.z};
double dir[3];
- rval = dagmc_ptr_->get_angle(surf, pnt, dir, &p->history());
+ moab::ErrorCode rval =
+ dagmc_ptr_->get_angle(mesh_handle(), pnt, dir, &p->history());
MB_CHK_ERR_CONT(rval);
- p->last_dir() = u.reflect(dir);
- return p->last_dir();
+ return u.reflect(dir);
}
//==============================================================================
diff --git a/src/particle.cpp b/src/particle.cpp
index 4dbab3213ac..5b46e2e639b 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -75,13 +75,6 @@ double Particle::speed() const
}
}
-void Particle::move_distance(double length)
-{
- for (int j = 0; j < n_coord(); ++j) {
- coord(j).r += length * coord(j).u;
- }
-}
-
void Particle::create_secondary(
double wgt, Direction u, double E, ParticleType type)
{
diff --git a/src/particle_data.cpp b/src/particle_data.cpp
index d8a5bb9d99c..fef8359a356 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<double>& rotation)
{
r = r.rotate(rotation);
@@ -56,6 +56,40 @@ GeometryState::GeometryState()
clear();
}
+void GeometryState::advance_to_boundary_from_void()
+{
+ auto root_coord = this->coord(0);
+ const auto& root_universe = model::universes[model::root_universe];
+ boundary().reset();
+
+ for (auto c_i : root_universe->cells_) {
+ auto dist =
+ model::cells.at(c_i)->distance(root_coord.r, root_coord.u, 0, this);
+ if (dist.first < boundary().distance) {
+ boundary().distance = dist.first;
+ boundary().surface = dist.second;
+ }
+ }
+
+ // if no intersection or near-infinite intersection, reset
+ // boundary information
+ if (boundary().distance > 1e300) {
+ boundary().distance = INFTY;
+ boundary().surface = SURFACE_NONE;
+ return;
+ }
+
+ // move the particle up to (and just past) the boundary
+ move_distance(boundary().distance + TINY_BIT);
+}
+
+void GeometryState::move_distance(double length)
+{
+ for (int j = 0; j < n_coord(); ++j) {
+ coord(j).r += length * coord(j).u;
+ }
+}
+
ParticleData::ParticleData()
{
zero_delayed_bank();
diff --git a/src/plot.cpp b/src/plot.cpp
index b36ed6f5d93..093f89f8cec 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -204,18 +204,21 @@ void read_plots_xml(pugi::xml_node root)
int id = std::stoi(id_string);
if (check_for_node(node, "type")) {
std::string type_str = get_node_value(node, "type", true);
- if (type_str == "slice")
+ if (type_str == "slice") {
model::plots.emplace_back(
std::make_unique<Plot>(node, Plot::PlotType::slice));
- else if (type_str == "voxel")
+ } else if (type_str == "voxel") {
model::plots.emplace_back(
std::make_unique<Plot>(node, Plot::PlotType::voxel));
- else if (type_str == "projection")
- model::plots.emplace_back(std::make_unique<ProjectionPlot>(node));
- else
+ } else if (type_str == "wireframe_raytrace") {
+ model::plots.emplace_back(
+ std::make_unique<WireframeRayTracePlot>(node));
+ } else if (type_str == "solid_raytrace") {
+ model::plots.emplace_back(std::make_unique<SolidRayTracePlot>(node));
+ } else {
fatal_error(
fmt::format("Unsupported plot type '{}' in plot {}", type_str, id));
-
+ }
model::plot_map[model::plots.back()->id()] = model::plots.size() - 1;
} else {
fatal_error(fmt::format("Must specify plot type in plot {}", id));
@@ -264,8 +267,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) {
@@ -1036,26 +1039,49 @@ 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};
+}
+
+WireframeRayTracePlot::WireframeRayTracePlot(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)
+void WireframeRayTracePlot::set_wireframe_color(pugi::xml_node plot_node)
{
// Copy plot background color
if (check_for_node(plot_node, "wireframe_color")) {
@@ -1068,7 +1094,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;
@@ -1089,33 +1115,7 @@ void ProjectionPlot::set_output_path(pugi::xml_node node)
path_plot_ = filename;
}
-// Advances to the next boundary from outside the geometry
-// Returns -1 if no intersection found, and the surface index
-// if an intersection was found.
-int ProjectionPlot::advance_to_boundary_from_void(GeometryState& p)
-{
- constexpr double scoot = 1e-5;
- double min_dist = {INFINITY};
- auto coord = p.coord(0);
- Universe* uni = model::universes[model::root_universe].get();
- int intersected_surface = -1;
- for (auto c_i : uni->cells_) {
- auto dist = model::cells.at(c_i)->distance(coord.r, coord.u, 0, &p);
- if (dist.first < min_dist) {
- min_dist = dist.first;
- intersected_surface = dist.second;
- }
- }
- if (min_dist > 1e300)
- return -1;
- else { // advance the particle
- for (int j = 0; j < p.n_coord(); ++j)
- p.coord(j).r += (min_dist + scoot) * p.coord(j).u;
- return std::abs(intersected_surface);
- }
-}
-
-bool ProjectionPlot::trackstack_equivalent(
+bool WireframeRayTracePlot::trackstack_equivalent(
const std::vector<TrackSegment>& track1,
const std::vector<TrackSegment>& track2) const
{
@@ -1125,7 +1125,7 @@ bool ProjectionPlot::trackstack_equivalent(
return false;
for (int i = 0; i < track1.size(); ++i) {
if (track1[i].id != track2[i].id ||
- track1[i].surface != track2[i].surface) {
+ track1[i].surface_index != track2[i].surface_index) {
return false;
}
}
@@ -1152,7 +1152,7 @@ bool ProjectionPlot::trackstack_equivalent(
if (t1_i == track1.size() && t2_i == track2.size())
break;
// Check if surface different
- if (track1[t1_i].surface != track2[t2_i].surface)
+ if (track1[t1_i].surface_index != track2[t2_i].surface_index)
return false;
// Pretty sure this should not be used:
@@ -1160,7 +1160,7 @@ bool ProjectionPlot::trackstack_equivalent(
// t1_i != track1.size() - 1 &&
// track1[t1_i+1].id != track2[t2_i+1].id) return false;
if (t2_i != 0 && t1_i != 0 &&
- track1[t1_i - 1].surface != track2[t2_i - 1].surface)
+ track1[t1_i - 1].surface_index != track2[t2_i - 1].surface_index)
return false;
// Check if neighboring cells are different
@@ -1176,44 +1176,58 @@ bool ProjectionPlot::trackstack_equivalent(
}
}
-void ProjectionPlot::create_output() const
+std::pair<Position, Direction> 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<double> 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.
+ // 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<double>(pixels_[0]);
double p1 = static_cast<double>(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<Position, Direction> result;
+
+ // Generate the starting position/direction of the ray
+ if (orthographic_width_ == C_NONE) { // 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<double>(horiz) - p0 / 2.0) / p0;
+ double y_pix_coord = (static_cast<double>(vert) - p1 / 2.0) / p1;
+
+ result.first = camera_position_ +
+ camera_y_axis() * x_pix_coord * orthographic_width_ +
+ camera_z_axis() * y_pix_coord * orthographic_width_;
+ result.second = camera_x_axis();
+ }
+
+ return result;
+}
+
+void WireframeRayTracePlot::create_output() const
+{
size_t width = pixels_[0];
size_t height = pixels_[1];
ImageData data({width, height}, not_found_);
- // This array marks where the initial wireframe was drawn.
- // We convolve it with a filter that gets adjusted with the
- // wireframe thickness in order to thicken the lines.
+ // This array marks where the initial wireframe was drawn. We convolve it with
+ // a filter that gets adjusted with the wireframe thickness in order to
+ // thicken the lines.
xt::xtensor<int, 2> wireframe_initial({width, height}, 0);
/* Holds all of the track segments for the current rendered line of pixels.
@@ -1242,16 +1256,13 @@ 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++) {
- // Save bottom line of current work chunk to compare against later
- // I used to have this inside the below if block, but it causes a
- // spurious line to be drawn at the bottom of the image. Not sure
- // why, but moving it here fixes things.
+ // Save bottom line of current work chunk to compare against later. This
+ // used to be inside the below if block, but it causes a spurious line to
+ // be drawn at the bottom of the image. Not sure why, but moving it here
+ // fixes things.
if (tid == n_threads - 1)
old_segments = this_line_segments[n_threads - 1];
@@ -1259,126 +1270,48 @@ 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<double>(horiz) - p0 / 2.0) / p0;
- double y_pix_coord = (static_cast<double>(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<Position, Direction> 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 = p.surface_index();
- 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));
-
- // 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;
- 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.
Position current_color(
not_found_.red, not_found_.green, not_found_.blue);
const auto& segments = this_line_segments[tid][horiz];
- for (unsigned i = segments.size(); i-- > 0;) {
+
+ // There must be at least two cell intersections to color, front and
+ // back of the cell. Maybe an infinitely thick cell could be present
+ // with no back, but why would you want to color that? It's easier to
+ // just skip that edge case and not even color it.
+ if (segments.size() <= 1)
+ continue;
+
+ for (int i = segments.size() - 2; i >= 0; --i) {
int colormap_idx = segments[i].id;
RGBColor seg_color = colors_[colormap_idx];
Position seg_color_vec(
seg_color.red, seg_color.green, seg_color.blue);
- double mixing = std::exp(-xs_[colormap_idx] * segments[i].length);
+ double mixing =
+ std::exp(-xs_[colormap_idx] *
+ (segments[i + 1].length - segments[i].length));
current_color =
current_color * mixing + (1.0 - mixing) * seg_color_vec;
- RGBColor result;
- result.red = static_cast<uint8_t>(current_color.x);
- result.green = static_cast<uint8_t>(current_color.y);
- result.blue = static_cast<uint8_t>(current_color.z);
- data(horiz, vert) = result;
}
+ // save result converting from double-precision color coordinates to
+ // byte-sized
+ RGBColor result;
+ result.red = static_cast<uint8_t>(current_color.x);
+ result.green = static_cast<uint8_t>(current_color.y);
+ result.blue = static_cast<uint8_t>(current_color.z);
+ data(horiz, vert) = result;
+
// Check to draw wireframe in horizontal direction. No inter-thread
// comm.
if (horiz > 0) {
@@ -1451,9 +1384,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);
@@ -1462,7 +1394,13 @@ void ProjectionPlot::print_info() const
fmt::print("Pixels: {} {}\n", pixels_[0], pixels_[1]);
}
-void ProjectionPlot::set_opacities(pugi::xml_node node)
+void WireframeRayTracePlot::print_info() const
+{
+ fmt::print("Plot Type: Wireframe ray-traced\n");
+ RayTracePlot::print_info();
+}
+
+void WireframeRayTracePlot::set_opacities(pugi::xml_node node)
{
xs_.resize(colors_.size(), 1e6); // set to large value for opaque by default
@@ -1492,7 +1430,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 =
@@ -1503,7 +1441,7 @@ void ProjectionPlot::set_orthographic_width(pugi::xml_node node)
}
}
-void ProjectionPlot::set_wireframe_thickness(pugi::xml_node node)
+void WireframeRayTracePlot::set_wireframe_thickness(pugi::xml_node node)
{
if (check_for_node(node, "wireframe_thickness")) {
int wireframe_thickness =
@@ -1514,7 +1452,7 @@ void ProjectionPlot::set_wireframe_thickness(pugi::xml_node node)
}
}
-void ProjectionPlot::set_wireframe_ids(pugi::xml_node node)
+void WireframeRayTracePlot::set_wireframe_ids(pugi::xml_node node)
{
if (check_for_node(node, "wireframe_ids")) {
wireframe_ids_ = get_node_array<int>(node, "wireframe_ids");
@@ -1529,7 +1467,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<int> pxls = get_node_array<int>(node, "pixels");
if (pxls.size() != 2)
@@ -1539,19 +1477,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<double> camera_pos = get_node_array<double>(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<double> look_at = get_node_array<double>(node, "look_at");
if (look_at.size() != 3) {
@@ -1562,7 +1500,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")) {
@@ -1576,6 +1514,352 @@ void ProjectionPlot::set_field_of_view(pugi::xml_node node)
}
}
+SolidRayTracePlot::SolidRayTracePlot(pugi::xml_node node) : RayTracePlot(node)
+{
+ set_opaque_ids(node);
+ set_diffuse_fraction(node);
+ set_light_position(node);
+}
+
+void SolidRayTracePlot::print_info() const
+{
+ fmt::print("Plot Type: Solid ray-traced\n");
+ RayTracePlot::print_info();
+}
+
+void SolidRayTracePlot::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<Position, Direction> 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 SolidRayTracePlot::set_opaque_ids(pugi::xml_node node)
+{
+ if (check_for_node(node, "opaque_ids")) {
+ auto opaque_ids_tmp = get_node_array<int>(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_tmp)
+ x = color_by_ == PlotColorBy::mats ? model::material_map[x]
+ : model::cell_map[x];
+
+ opaque_ids_.insert(opaque_ids_tmp.begin(), opaque_ids_tmp.end());
+ }
+}
+
+void SolidRayTracePlot::set_light_position(pugi::xml_node node)
+{
+ if (check_for_node(node, "light_position")) {
+ auto light_pos_tmp = get_node_array<double>(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 SolidRayTracePlot::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()
+{
+ boundary() = distance_to_boundary(*this);
+}
+
+void Ray::trace()
+{
+ // To trace the ray from its origin all the way through the model, we have
+ // to proceed in two phases. In the first, the ray may or may not be found
+ // inside the model. If the ray is already in the model, phase one can be
+ // skipped. Otherwise, the ray has to be advanced to the boundary of the
+ // model where all the cells are defined. Importantly, this is assuming that
+ // the model is convex, which is a very reasonable assumption for any
+ // radiation transport model.
+ //
+ // After phase one is done, we can starting tracing from cell to cell within
+ // the model. This step can use neighbor lists to accelerate the ray tracing.
+
+ // Attempt to initialize the particle. We may have to enter a loop to move
+ // it up to the edge of the model.
+ bool inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
+
+ // Advance to the boundary of the model
+ while (!inside_cell) {
+ advance_to_boundary_from_void();
+ inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
+
+ // If true 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<int>::max()) {
+ warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u()));
+ return;
+ }
+
+ // Exit this loop and enter into cell-to-cell ray tracing (which uses
+ // neighbor lists)
+ if (inside_cell)
+ break;
+
+ // if there is no intersection with the model, we're done
+ if (boundary().surface == SURFACE_NONE)
+ return;
+
+ event_counter_++;
+ if (event_counter_ > MAX_INTERSECTIONS) {
+ warning("Likely infinite loop in ray traced plot");
+ return;
+ }
+ }
+
+ // Call the specialized logic for this type of ray. This is for the
+ // intersection for the first intersection if we had one.
+ if (boundary().surface != SURFACE_NONE) {
+ // set the geometry state's surface attribute to be used for
+ // surface normal computation
+ surface() = boundary().surface;
+ on_intersection();
+ if (stop_)
+ return;
+ }
+
+ // reset surface attribute to zero after the first intersection so that it
+ // doesn't perturb surface crossing logic from here on out
+ surface() = 0;
+
+ // This is the ray tracing loop within the model. It exits after exiting
+ // the model, which is equivalent to assuming that the model is convex.
+ // It would be nice to factor out the on_intersection at the end of this
+ // loop and then do "while (inside_cell)", but we can't guarantee it's
+ // on a surface in that case. There might be some other way to set it
+ // up that is perhaps a little more elegant, but this is what works just
+ // fine.
+ while (true) {
+
+ compute_distance();
+
+ // There are no more intersections to process
+ // if we hit the edge of the model, so stop
+ // the particle in that case. Also, just exit
+ // if a negative distance was somehow computed.
+ if (boundary().distance == INFTY || boundary().distance == INFINITY ||
+ boundary().distance < 0) {
+ return;
+ }
+
+ // See below comment where call_on_intersection is checked in an
+ // if statement for an explanation of this.
+ bool call_on_intersection {true};
+ if (boundary().distance < 10 * TINY_BIT) {
+ call_on_intersection = false;
+ }
+
+ // DAGMC surfaces expect us to go a little bit further than the advance
+ // distance to properly check cell inclusion.
+ boundary().distance += TINY_BIT;
+
+ // Advance particle, prepare for next intersection
+ for (int lev = 0; lev < n_coord(); ++lev) {
+ coord(lev).r += boundary().distance * coord(lev).u;
+ }
+ surface() = boundary().surface;
+ n_coord_last() = n_coord();
+ n_coord() = boundary().coord_level;
+ if (boundary().lattice_translation[0] != 0 ||
+ boundary().lattice_translation[1] != 0 ||
+ boundary().lattice_translation[2] != 0) {
+ cross_lattice(*this, boundary(), settings::verbosity >= 10);
+ }
+
+ // Record how far the ray has traveled
+ traversal_distance_ += boundary().distance;
+ inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
+
+ // Call the specialized logic for this type of ray. Note that we do not
+ // call this if the advance distance is very small. Unfortunately, it seems
+ // darn near impossible to get the particle advanced to the model boundary
+ // and through it without sometimes accidentally calling on_intersection
+ // twice. This incorrectly shades the region as occluded when it might not
+ // actually be. By screening out intersection distances smaller than a
+ // threshold 10x larger than the scoot distance used to advance up to the
+ // model boundary, we can avoid that situation.
+ if (call_on_intersection) {
+ on_intersection();
+ if (stop_)
+ return;
+ }
+
+ if (!inside_cell)
+ return;
+
+ event_counter_++;
+ if (event_counter_ > MAX_INTERSECTIONS) {
+ warning("Likely infinite loop in ray traced plot");
+ return;
+ }
+ }
+}
+
+void ProjectionRay::on_intersection()
+{
+ // 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,
+ traversal_distance_, boundary().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_.opaque_ids_.find(hit_id) == plot_.opaque_ids_.end())
+ 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();
+
+ // TODO
+ // Not sure what can cause a surface token to be invalid here, 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
+ // ray.
+ if (surface() == 0) {
+ result_color_ = plot_.overlap_color_;
+ stop();
+ return;
+ }
+
+ // Get surface pointer
+ const auto& surf = model::surfaces.at(surface_index());
+
+ Direction normal = surf->normal(r_local());
+ normal /= normal.norm();
+
+ // Need to apply translations to find the normal vector in
+ // 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 = normal.inverse_rotate(c.rotation_);
+ }
+ }
+
+ // use the normal opposed to the ray direction
+ if (normal.dot(u()) > 0.0) {
+ normal *= -1.0;
+ }
+
+ // Facing away from the light means no lighting
+ double dotprod = normal.dot(to_light);
+ dotprod = std::max(0.0, dotprod);
+
+ 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;
+
+ // OpenMC native CSG and DAGMC surfaces have some slight differences
+ // in how they interpret particles that are sitting on a surface.
+ // I don't know exactly why, but this makes everything work beautifully.
+ if (surf->geom_type() == GeometryType::DAG) {
+ surface() = 0;
+ } else {
+ surface() = -surface(); // go to other side
+ }
+
+ // Must fully restart coordinate search. Why? Not sure.
+ clear();
+
+ // Note this could likely be faster if we cached the previous
+ // cell we were in before the reflection. This is the easiest
+ // way to fully initialize all the sub-universe coordinates and
+ // directions though.
+ 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/src/position.cpp b/src/position.cpp
index 5b3613b2897..0361e99a0a9 100644
--- a/src/position.cpp
+++ b/src/position.cpp
@@ -75,13 +75,6 @@ Position Position::operator-() const
return {-x, -y, -z};
}
-Position Position::rotate(const vector<double>& rotation) const
-{
- return {x * rotation[0] + y * rotation[1] + z * rotation[2],
- x * rotation[3] + y * rotation[4] + z * rotation[5],
- x * rotation[6] + y * rotation[7] + z * rotation[8]};
-}
-
std::ostream& operator<<(std::ostream& os, Position r)
{
os << "(" << r.x << ", " << r.y << ", " << r.z << ")";
diff --git a/src/surface.cpp b/src/surface.cpp
index dbcaf849848..30a8d1f211f 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -141,7 +141,7 @@ Direction Surface::reflect(Position r, Direction u, GeometryState* p) const
}
Direction Surface::diffuse_reflect(
- Position r, Direction u, uint64_t* seed, GeometryState* p) const
+ Position r, Direction u, uint64_t* seed) const
{
// Diffuse reflect direction according to the normal.
// cosine distribution
diff --git a/tests/regression_tests/plot_projections/plots.xml b/tests/regression_tests/plot_projections/plots.xml
index 85689da1475..50d129cd97b 100644
--- a/tests/regression_tests/plot_projections/plots.xml
+++ b/tests/regression_tests/plot_projections/plots.xml
@@ -1,7 +1,7 @@
<?xml version="1.0"?>
<plots>
- <plot id="1" type="projection">
+ <plot id="1" type="wireframe_raytrace">
<look_at>0. 0. 0.</look_at>
<camera_position>20. 20. 20.</camera_position>
<pixels>200 200</pixels>
@@ -11,7 +11,7 @@
<field_of_view>70</field_of_view>
</plot>
- <plot id="2" type="projection">
+ <plot id="2" type="wireframe_raytrace">
<look_at>0. 0. 0.</look_at>
<camera_position>10. 10. 0.</camera_position>
<width>25 25</width>
@@ -22,7 +22,7 @@
<filename>example1</filename>
</plot>
- <plot id="3" color_by="material" type="projection">
+ <plot id="3" color_by="material" type="wireframe_raytrace">
<look_at>0. 0. 0.</look_at>
<camera_position>20. 20. 20.</camera_position>
<pixels>200 200</pixels>
@@ -31,7 +31,7 @@
<wireframe_ids>2</wireframe_ids>
</plot>
- <plot id="4" color_by="material" type="projection">
+ <plot id="4" color_by="material" type="wireframe_raytrace">
<look_at>0. 0. 0.</look_at>
<camera_position>0. 10.0 20.</camera_position>
<pixels>200 200</pixels>
@@ -39,7 +39,7 @@
<filename>example3.png</filename>
</plot>
- <plot id="5" type="projection">
+ <plot id="5" type="wireframe_raytrace">
<look_at>0. 0. 0.</look_at>
<camera_position>10. 10. 10.</camera_position>
<width>25 25</width>
@@ -52,4 +52,39 @@
<color id="3" rgb="255 0 0" xs="1.0"/>
</plot>
+ <plot id="6" color_by="material" type="solid_raytrace">
+ <look_at>0. 0. 0.</look_at>
+ <camera_position>10. 10. 10.</camera_position>
+ <pixels>200 200</pixels>
+ <filename>phong.png</filename>
+ <opaque_ids>1 3</opaque_ids>
+ <color id="1" rgb="0 0 255"/>
+ <color id="2" rgb="0 255 0"/>
+ <color id="3" rgb="255 0 100"/>
+ </plot>
+
+ <plot id="7" color_by="material" type="solid_raytrace">
+ <look_at>0. 0. 0.</look_at>
+ <camera_position>10. 10. 10.</camera_position>
+ <diffuse_fraction>0.5</diffuse_fraction>
+ <pixels>200 200</pixels>
+ <filename>phong_diffuse.png</filename>
+ <opaque_ids>1 3</opaque_ids>
+ <color id="1" rgb="0 0 255"/>
+ <color id="2" rgb="0 255 0"/>
+ <color id="3" rgb="255 0 100"/>
+ </plot>
+
+ <plot id="8" color_by="material" type="solid_raytrace">
+ <look_at>0. 0. 0.</look_at>
+ <camera_position>10. 10. 10.</camera_position>
+ <light_position>0. 10. 10.</light_position>
+ <pixels>200 200</pixels>
+ <filename>phong_move_light.png</filename>
+ <opaque_ids>1 3</opaque_ids>
+ <color id="1" rgb="0 0 255"/>
+ <color id="2" rgb="0 255 0"/>
+ <color id="3" rgb="255 0 100"/>
+ </plot>
+
</plots>
diff --git a/tests/regression_tests/plot_projections/results_true.dat b/tests/regression_tests/plot_projections/results_true.dat
index 1a67da3002a..672280fdd1d 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
+025804f1522eafd6e0e9566ce6b9b5603962f278de222c842fe3e06471290bb575676255bcd55e4f084bdcca4ee56d3c219827cb1ef2b5c3a90f7666986b55e9
\ 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()
diff --git a/tests/unit_tests/test_plots.py b/tests/unit_tests/test_plots.py
index 4efc12c9326..fad574ee697 100644
--- a/tests/unit_tests/test_plots.py
+++ b/tests/unit_tests/test_plots.py
@@ -4,6 +4,8 @@
import openmc.examples
import pytest
+from openmc.plots import _SVG_COLORS
+
@pytest.fixture(scope='module')
def myplot():
@@ -42,7 +44,7 @@ def myplot():
@pytest.fixture(scope='module')
def myprojectionplot():
- plot = openmc.ProjectionPlot(name='myprojectionplot')
+ plot = openmc.WireframeRayTracePlot(name='myprojectionplot')
plot.look_at = (0.0, 0.0, 0.0)
plot.camera_position = (4.0, 3.0, 0.0)
plot.pixels = (500, 500)
@@ -118,6 +120,31 @@ def test_repr_proj(myprojectionplot):
assert isinstance(r, str)
+def test_projection_plot_roundtrip(myprojectionplot):
+
+ elem = myprojectionplot.to_xml_element()
+
+ xml_plot = openmc.WireframeRayTracePlot.from_xml_element(elem)
+
+ svg_colors = _SVG_COLORS
+
+ assert xml_plot.name == myprojectionplot.name
+ assert xml_plot.look_at == myprojectionplot.look_at
+ assert xml_plot.camera_position == myprojectionplot.camera_position
+ assert xml_plot.pixels == myprojectionplot.pixels
+ assert xml_plot.filename == myprojectionplot.filename
+ assert xml_plot.background == svg_colors[myprojectionplot.background]
+ assert xml_plot.color_by == myprojectionplot.color_by
+ expected_colors = {m.id: svg_colors[c] for m, c in myprojectionplot.colors.items()}
+ assert xml_plot.colors == expected_colors
+ # TODO: needs geometry information
+ # assert xml_plot.mask_components == myprojectionplot.mask_components
+ assert xml_plot.mask_background == svg_colors[myprojectionplot.mask_background]
+ # assert xml_plot.overlap_color == svg_colors[myprojectionplot.overlap_color]
+ assert xml_plot.wireframe_thickness == myprojectionplot.wireframe_thickness
+ assert xml_plot.level == myprojectionplot.level
+
+
def test_from_geometry():
width = 25.
s = openmc.Sphere(r=width/2, boundary_type='vacuum')
@@ -182,7 +209,7 @@ def test_plots(run_in_tmpdir):
plots = openmc.Plots([p1, p2])
assert len(plots) == 2
- p3 = openmc.ProjectionPlot(name='plot3')
+ p3 = openmc.WireframeRayTracePlot(name='plot3')
plots = openmc.Plots([p1, p2, p3])
assert len(plots) == 3
@@ -223,6 +250,41 @@ def test_voxel_plot_roundtrip():
assert new_plot.color_by == plot.color_by
+def test_phong_plot_roundtrip():
+ plot = openmc.SolidRayTracePlot(name='my phong plot')
+ plot.id = 2300
+ plot.filename = 'phong1'
+ plot.pixels = (50, 50)
+ plot.look_at = (11., 12., 13.)
+ plot.camera_position = (22., 23., 24.)
+ plot.diffuse_fraction = 0.5
+ plot.horizontal_field_of_view = 90.0
+ plot.color_by = 'material'
+ plot.light_position = (8., 9., 10.)
+ plot.opaque_domains = [6, 7, 8]
+
+ elem = plot.to_xml_element()
+
+ repr(plot)
+
+ new_plot = openmc.SolidRayTracePlot.from_xml_element(elem)
+
+ assert new_plot.name == plot.name
+ assert new_plot.id == plot.id
+ assert new_plot.filename == plot.filename
+ assert new_plot.pixels == plot.pixels
+ assert new_plot.look_at == plot.look_at
+ assert new_plot.camera_position == plot.camera_position
+ assert new_plot.diffuse_fraction == plot.diffuse_fraction
+ assert new_plot.horizontal_field_of_view == plot.horizontal_field_of_view
+ assert new_plot.color_by == plot.color_by
+ assert new_plot.light_position == plot.light_position
+ assert new_plot.opaque_domains == plot.opaque_domains
+
+ # ensure the new object is valid to re-write to XML
+ new_elem = new_plot.to_xml_element()
+
+
def test_plot_directory(run_in_tmpdir):
pwr_pin = openmc.examples.pwr_pin_cell()