#ifndef OPENMC_PLOT_H #define OPENMC_PLOT_H #include #include #include #include #include #include "openmc/tensor.h" #include "pugixml.hpp" #include "hdf5.h" #include "openmc/cell.h" #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/geometry.h" #include "openmc/particle.h" #include "openmc/position.h" #include "openmc/random_lcg.h" #include "openmc/ray.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_match.h" #include "openmc/xml_interface.h" namespace openmc { //=============================================================================== // Global variables //=============================================================================== class PlottableInterface; namespace model { extern std::unordered_map plot_map; //!< map of plot ids to index extern vector> plots; //!< Plot instance container extern uint64_t plotter_seed; // Stream index used by the plotter } // namespace model //=============================================================================== // RGBColor holds color information for plotted objects //=============================================================================== struct RGBColor { // Constructors RGBColor() : red(0), green(0), blue(0) {}; RGBColor(const int v[3]) : red(v[0]), green(v[1]), blue(v[2]) {}; RGBColor(int r, int g, int b) : red(r), green(g), blue(b) {}; RGBColor(const vector& v) { if (v.size() != 3) { throw std::out_of_range("Incorrect vector size for RGBColor."); } red = v[0]; green = v[1]; blue = v[2]; } bool operator==(const RGBColor& other) { 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; }; // some default colors const RGBColor WHITE {255, 255, 255}; const RGBColor RED {255, 0, 0}; const RGBColor BLACK {0, 0, 0}; /** * \class PlottableInterface * \brief Interface for plottable objects. * * PlottableInterface classes must have unique IDs. If no ID (or -1) is * provided, the next available ID is assigned automatically. 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. */ typedef tensor::Tensor ImageData; class PlottableInterface { public: PlottableInterface() = default; void set_default_colors(); private: void set_id(pugi::xml_node plot_node); int id_ {C_NONE}; // unique plot ID void set_bg_color(pugi::xml_node plot_node); void set_universe(pugi::xml_node plot_node); void set_color_by(pugi::xml_node plot_node); void set_user_colors(pugi::xml_node plot_node); void set_overlap_color(pugi::xml_node plot_node); void set_mask(pugi::xml_node plot_node); protected: // Plot output filename, derived classes have logic to set it std::string path_plot_; public: enum class PlotColorBy { cells = 0, mats = 1 }; // Generates image data based on plot parameters and returns it virtual ImageData create_image() const = 0; // Creates the output image named path_plot_ virtual void create_output() const = 0; // Write populated image data to file void write_image(const ImageData& data) const; // Print useful info to the terminal virtual void print_info() const = 0; const std::string& path_plot() const { return path_plot_; } std::string& path_plot() { return path_plot_; } int id() const { return id_; } void set_id(int id = C_NONE); int level() const { return level_; } PlotColorBy color_by() const { return color_by_; } // Public color-related data PlottableInterface(pugi::xml_node plot_node); virtual ~PlottableInterface() = default; int level_ {-1}; // Universe level to plot bool color_overlaps_ {false}; // Show overlapping cells? PlotColorBy color_by_ {PlotColorBy::mats}; // Plot coloring (cell/material) RGBColor not_found_ {WHITE}; // Plot background color RGBColor overlap_color_ {RED}; // Plot overlap color vector colors_; // Plot colors }; struct IdData { // Constructor IdData(size_t h_res, size_t v_res, bool include_filter = false); // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); void set_overlap(size_t y, size_t x, int overlap_idx); // Members tensor::Tensor data_; //!< 2D array of cell & material ids }; struct PropertyData { // Constructor PropertyData(size_t h_res, size_t v_res, bool include_filter = false); // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); void set_overlap(size_t y, size_t x, int overlap_idx); // Members tensor::Tensor data_; //!< 2D array of temperature & density data }; struct RasterData { // Constructor RasterData(size_t h_res, size_t v_res, bool include_filter = false); // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); void set_overlap(size_t y, size_t x, int overlap_idx); // Members tensor::Tensor id_data_; //!< [v_res, h_res, 3 or 4]: cell, instance, mat, [filter_bin] tensor::Tensor property_data_; //!< [v_res, h_res, 2]: temperature, density bool include_filter_; //!< Whether filter bin index is included }; //=============================================================================== // Plot class //=============================================================================== class SlicePlotBase { public: template T get_map(int32_t filter_index = -1) const; enum class PlotBasis { xy = 1, xz = 2, yz = 3 }; // Accessors const std::array& pixels() const { return pixels_; } std::array& pixels() { return pixels_; } // Members public: Position origin_; //!< Plot origin in geometry Direction u_span_; //!< Full-width span vector in geometry Direction v_span_; //!< Full-height span vector in geometry array pixels_; //!< Plot size in pixels bool show_overlaps_; //!< Show overlapping cells? int slice_level_ {-1}; //!< Plot universe level private: }; template T SlicePlotBase::get_map(int32_t filter_index) const { size_t width = pixels_[0]; size_t height = pixels_[1]; // Determine if filter is being used bool include_filter = (filter_index >= 0); Filter* filter = nullptr; if (include_filter) { filter = model::tally_filters[filter_index].get(); } // size data array T data(width, height, include_filter); // compute pixel steps and top-left pixel center Direction u_step = u_span_ / static_cast(width); Direction v_step = v_span_ / static_cast(height); Position start = origin_ - 0.5 * u_span_ + 0.5 * v_span_ + 0.5 * u_step - 0.5 * v_step; // Validate that span vectors define a valid plane Position cross = u_span_.cross(v_span_); if (cross.norm() == 0.0) { fatal_error("Slice span vectors are invalid (zero area)."); } // Use an arbitrary direction that is not aligned with any coordinate axis. // The direction has no physical meaning for plotting but is used by // Surface::sense() to break ties when a pixel is coincident with a surface. Direction dir = {1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 0.0}; #pragma omp parallel { Particle p; p.r() = start; p.u() = dir; p.coord(0).universe() = model::root_universe; int level = slice_level_; int j {}; FilterMatch match; #pragma omp for for (int y = 0; y < height; y++) { Position row = start - v_step * static_cast(y); for (int x = 0; x < width; x++) { p.r() = row + u_step * static_cast(x); p.n_coord() = 1; // local variables bool found_cell = exhaustive_find_cell(p); j = p.n_coord() - 1; if (level >= 0) { j = level; } if (found_cell) { data.set_value(y, x, p, j, filter, &match); } if (show_overlaps_) { int overlap_idx = check_cell_overlap(p, false); if (overlap_idx >= 0) { data.set_overlap(y, x, overlap_idx); } } } // inner for } } return data; } // Represents either a voxel or pixel plot class Plot : public PlottableInterface, public SlicePlotBase { public: enum class PlotType { slice = 1, voxel = 2 }; Plot(pugi::xml_node plot, PlotType type); private: void set_output_path(pugi::xml_node plot_node); void set_basis(pugi::xml_node plot_node); void set_origin(pugi::xml_node plot_node); void set_width(pugi::xml_node plot_node); void set_meshlines(pugi::xml_node plot_node); public: // Add mesh lines to ImageData void draw_mesh_lines(ImageData& data) const; ImageData create_image() const override; void create_voxel() const; void create_output() const override; void print_info() const override; PlotType type_; //!< Plot type (Slice/Voxel) Position width_; //!< Axis-aligned width from plot.xml PlotBasis basis_; //!< Basis from plot.xml for slice plots int meshlines_width_; //!< Width of lines added to the plot int index_meshlines_mesh_ {-1}; //!< Index of the mesh to draw on the plot RGBColor meshlines_color_; //!< Color of meshlines on the plot }; /** * \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: RayTracePlot() = default; RayTracePlot(pugi::xml_node plot); // Standard getters. No setting since it's done from XML. const Position& camera_position() const { return camera_position_; } Position& camera_position() { return camera_position_; } const Position& look_at() const { return look_at_; } Position& look_at() { return look_at_; } const double& horizontal_field_of_view() const { return horizontal_field_of_view_; } double& horizontal_field_of_view() { return horizontal_field_of_view_; } void print_info() const override; const std::array& pixels() const { return pixels_; } std::array& pixels() { return pixels_; } const Direction& up() const { return up_; } Direction& up() { return up_; } //! brief Updates the cached camera-to-model matrix after changes to //! camera parameters. void update_view(); 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 get_pixel_ray(int horiz, int vert) const; 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_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 std::array pixels_ {100, 100}; // pixel dimension of resulting image 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 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); ImageData create_image() const override; void create_output() const override; void print_info() const override; 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); /* Checks if a vector of two TrackSegments is equivalent. We define this * to mean not having matching intersection lengths, but rather having * a matching sequence of surface/cell/material intersections. */ struct TrackSegment; bool trackstack_equivalent(const vector& track1, const vector& track2) const; /* Used for drawing wireframe and colors. We record the list of * surface/cell/material intersections and the corresponding lengths as a ray * traverses the geometry, then color by iterating in reverse. */ struct TrackSegment { int id; // material or cell ID (which is being colored) double length; // length of this track intersection /* Recording this allows us to draw edges on the wireframe. For instance * if two surfaces bound a single cell, it allows drawing that sharp edge * where the surfaces intersect. */ int surface_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_index(surface_a) {} }; // which color IDs should be wireframed. If empty, all cells are wireframed. vector 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 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() = default; SolidRayTracePlot(pugi::xml_node plot); ImageData create_image() const override; void create_output() const override; void print_info() const override; const std::unordered_set& opaque_ids() const { return opaque_ids_; } std::unordered_set& opaque_ids() { return opaque_ids_; } const Position& light_location() const { return light_location_; } Position& light_location() { return light_location_; } const double& diffuse_fraction() const { return diffuse_fraction_; } double& diffuse_fraction() { return diffuse_fraction_; } 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 opaque_ids_; double diffuse_fraction_ {0.1}; // By default, the light is at the camera unless otherwise specified. Position light_location_; }; class ProjectionRay : public Ray { public: ProjectionRay(Position r, Direction u, const WireframeRayTracePlot& plot, vector& line_segments) : Ray(r, u), plot_(plot), line_segments_(line_segments) {} void on_intersection() override; private: /* Store a reference to the plot object which is running this ray, in order * to access some of the plot settings which influence the behavior where * intersections are. */ const WireframeRayTracePlot& plot_; /* The ray runs through the geometry, and records the lengths of ray segments * and cells they lie in along the way. */ vector& line_segments_; }; class PhongRay : public Ray { public: PhongRay(Position r, Direction u, const SolidRayTracePlot& plot) : Ray(r, u), plot_(plot) { result_color_ = plot_.not_found_; } 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_; }; //=============================================================================== // Non-member functions //=============================================================================== /* Write a PPM image * filename - name of output file * data - image data to write */ void output_ppm(const std::string& filename, const ImageData& data); #ifdef USE_LIBPNG /* Write a PNG image * filename - name of output file * data - image data to write */ void output_png(const std::string& filename, const ImageData& data); #endif //! Initialize a voxel file //! \param[in] id of an open hdf5 file //! \param[in] dimensions of the voxel file (dx, dy, dz) //! \param[out] dataspace pointer to voxel data //! \param[out] dataset pointer to voxesl data //! \param[out] pointer to memory space of voxel data void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset, hid_t* memspace); //! Write a section of the voxel data to hdf5 //! \param[in] voxel slice //! \param[out] dataspace pointer to voxel data //! \param[out] dataset pointer to voxesl data //! \param[out] pointer to data to write void voxel_write_slice( int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf); //! Close voxel file entities //! \param[in] data space to close //! \param[in] dataset to close //! \param[in] memory space to close void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace); //=============================================================================== // External functions //=============================================================================== //! Read plot specifications from a plots.xml file void read_plots_xml(); //! Read plot specifications from an XML Node //! \param[in] XML node containing plot info void read_plots_xml(pugi::xml_node root); //! Clear memory void free_memory_plot(); //! Create a randomly generated RGB color //! \return RGBColor with random value RGBColor random_color(); } // namespace openmc #endif // OPENMC_PLOT_H