crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-23 19:45:34 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5
OpenMC/src/plot.cpp
Olaf Schumann 4905acc20a fix memory leak in output_png()
2022-01-13 14:29:20 +01:00

1024 lines
30 KiB
C++
Raw Blame History

#include "openmc/plot.h"
#include <algorithm>
#include <cstdio>
#include <fstream>
#include <sstream>
#include "xtensor/xview.hpp"
#include <fmt/core.h>
#include <fmt/ostream.h>
#ifdef USE_LIBPNG
#include <png.h>
#endif
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
#include "openmc/material.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/output.h"
#include "openmc/particle.h"
#include "openmc/progress_bar.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level
constexpr int32_t NOT_FOUND {-2};
constexpr int32_t OVERLAP {-3};
IdData::IdData(size_t h_res, size_t v_res) : data_({v_res, h_res, 3}, NOT_FOUND)
{}
void IdData::set_value(size_t y, size_t x, const Particle& p, int level)
{
// set cell data
if (p.n_coord() <= level) {
data_(y, x, 0) = NOT_FOUND;
data_(y, x, 1) = NOT_FOUND;
} else {
data_(y, x, 0) = model::cells.at(p.coord(level).cell)->id_;
data_(y, x, 1) = level == p.n_coord() - 1
? p.cell_instance()
: cell_instance_at_level(p, level);
}
// set material data
Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get();
if (p.material() == MATERIAL_VOID) {
data_(y, x, 2) = MATERIAL_VOID;
return;
} else if (c->type_ == Fill::MATERIAL) {
Material* m = model::materials.at(p.material()).get();
data_(y, x, 2) = m->id_;
}
}
void IdData::set_overlap(size_t y, size_t x)
{
xt::view(data_, y, x, xt::all()) = OVERLAP;
}
PropertyData::PropertyData(size_t h_res, size_t v_res)
: data_({v_res, h_res, 2}, NOT_FOUND)
{}
void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level)
{
Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get();
data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
Material* m = model::materials.at(p.material()).get();
data_(y, x, 1) = m->density_gpcc_;
}
}
void PropertyData::set_overlap(size_t y, size_t x)
{
data_(y, x) = OVERLAP;
}
//==============================================================================
// Global variables
//==============================================================================
namespace model {
std::unordered_map<int, int> plot_map;
vector<Plot> plots;
uint64_t plotter_seed = 1;
} // namespace model
//==============================================================================
// RUN_PLOT controls the logic for making one or many plots
//==============================================================================
extern "C" int openmc_plot_geometry()
{
for (auto& pl : model::plots) {
write_message(5, "Processing plot {}: {}...", pl.id_, pl.path_plot_);
if (PlotType::slice == pl.type_) {
// create 2D image
create_image(pl);
} else if (PlotType::voxel == pl.type_) {
// create voxel file for 3D viewing
create_voxel(pl);
}
}
return 0;
}
void read_plots_xml()
{
// Check if plots.xml exists; this is only necessary when the plot runmode is
// initiated. Otherwise, we want to read plots.xml because it may be called
// later via the API. In that case, its ok for a plots.xml to not exist
std::string filename = settings::path_input + "plots.xml";
if (!file_exists(filename) && settings::run_mode == RunMode::PLOTTING) {
fatal_error(fmt::format("Plots XML file '{}' does not exist!", filename));
}
write_message("Reading plot XML file...", 5);
// Parse plots.xml file
pugi::xml_document doc;
doc.load_file(filename.c_str());
pugi::xml_node root = doc.document_element();
for (auto node : root.children("plot")) {
model::plots.emplace_back(node);
model::plot_map[model::plots.back().id_] = model::plots.size() - 1;
}
}
void free_memory_plot()
{
model::plots.clear();
model::plot_map.clear();
}
//==============================================================================
// CREATE_IMAGE creates an image based on user input from a plots.xml <plot>
// specification in the PNG/PPM format
//==============================================================================
void create_image(Plot const& pl)
{
size_t width = pl.pixels_[0];
size_t height = pl.pixels_[1];
ImageData data({width, height}, pl.not_found_);
// generate ids for the plot
auto ids = pl.get_map<IdData>();
// assign colors
for (size_t y = 0; y < height; y++) {
for (size_t x = 0; x < width; x++) {
int idx = pl.color_by_ == PlotColorBy::cells ? 0 : 2;
auto id = ids.data_(y, x, idx);
// no setting needed if not found
if (id == NOT_FOUND) {
continue;
}
if (id == OVERLAP) {
data(x, y) = pl.overlap_color_;
continue;
}
if (PlotColorBy::cells == pl.color_by_) {
data(x, y) = pl.colors_[model::cell_map[id]];
} else if (PlotColorBy::mats == pl.color_by_) {
if (id == MATERIAL_VOID) {
data(x, y) = WHITE;
continue;
}
data(x, y) = pl.colors_[model::material_map[id]];
} // color_by if-else
} // x for loop
} // y for loop
// draw mesh lines if present
if (pl.index_meshlines_mesh_ >= 0) {
draw_mesh_lines(pl, data);
}
// create image file
#ifdef USE_LIBPNG
output_png(pl, data);
#else
output_ppm(pl, data);
#endif
}
void Plot::set_id(pugi::xml_node plot_node)
{
// Copy data into plots
if (check_for_node(plot_node, "id")) {
id_ = std::stoi(get_node_value(plot_node, "id"));
} else {
fatal_error("Must specify plot id in plots XML file.");
}
// Check to make sure 'id' hasn't been used
if (model::plot_map.find(id_) != model::plot_map.end()) {
fatal_error(
fmt::format("Two or more plots use the same unique ID: {}", id_));
}
}
void Plot::set_type(pugi::xml_node plot_node)
{
// Copy plot type
// Default is slice
type_ = PlotType::slice;
// check type specified on plot node
if (check_for_node(plot_node, "type")) {
std::string type_str = get_node_value(plot_node, "type", true);
// set type using node value
if (type_str == "slice") {
type_ = PlotType::slice;
} else if (type_str == "voxel") {
type_ = PlotType::voxel;
} else {
// if we're here, something is wrong
fatal_error(
fmt::format("Unsupported plot type '{}' in plot {}", type_str, id_));
}
}
}
void Plot::set_output_path(pugi::xml_node plot_node)
{
// Set output file path
std::string filename;
if (check_for_node(plot_node, "filename")) {
filename = get_node_value(plot_node, "filename");
} else {
filename = fmt::format("plot_{}", id_);
}
// add appropriate file extension to name
switch (type_) {
case PlotType::slice:
#ifdef USE_LIBPNG
filename.append(".png");
#else
filename.append(".ppm");
#endif
break;
case PlotType::voxel:
filename.append(".h5");
break;
}
path_plot_ = filename;
// Copy plot pixel size
vector<int> pxls = get_node_array<int>(plot_node, "pixels");
if (PlotType::slice == type_) {
if (pxls.size() == 2) {
pixels_[0] = pxls[0];
pixels_[1] = pxls[1];
} else {
fatal_error(
fmt::format("<pixels> must be length 2 in slice plot {}", id_));
}
} else if (PlotType::voxel == type_) {
if (pxls.size() == 3) {
pixels_[0] = pxls[0];
pixels_[1] = pxls[1];
pixels_[2] = pxls[2];
} else {
fatal_error(
fmt::format("<pixels> must be length 3 in voxel plot {}", id_));
}
}
}
void Plot::set_bg_color(pugi::xml_node plot_node)
{
// Copy plot background color
if (check_for_node(plot_node, "background")) {
vector<int> bg_rgb = get_node_array<int>(plot_node, "background");
if (PlotType::voxel == type_) {
if (mpi::master) {
warning(fmt::format("Background color ignored in voxel plot {}", id_));
}
}
if (bg_rgb.size() == 3) {
not_found_ = bg_rgb;
} else {
fatal_error(fmt::format("Bad background RGB in plot {}", id_));
}
}
}
void Plot::set_basis(pugi::xml_node plot_node)
{
// Copy plot basis
if (PlotType::slice == type_) {
std::string pl_basis = "xy";
if (check_for_node(plot_node, "basis")) {
pl_basis = get_node_value(plot_node, "basis", true);
}
if ("xy" == pl_basis) {
basis_ = PlotBasis::xy;
} else if ("xz" == pl_basis) {
basis_ = PlotBasis::xz;
} else if ("yz" == pl_basis) {
basis_ = PlotBasis::yz;
} else {
fatal_error(
fmt::format("Unsupported plot basis '{}' in plot {}", pl_basis, id_));
}
}
}
void Plot::set_origin(pugi::xml_node plot_node)
{
// Copy plotting origin
auto pl_origin = get_node_array<double>(plot_node, "origin");
if (pl_origin.size() == 3) {
origin_ = pl_origin;
} else {
fatal_error(fmt::format("Origin must be length 3 in plot {}", id_));
}
}
void Plot::set_width(pugi::xml_node plot_node)
{
// Copy plotting width
vector<double> pl_width = get_node_array<double>(plot_node, "width");
if (PlotType::slice == type_) {
if (pl_width.size() == 2) {
width_.x = pl_width[0];
width_.y = pl_width[1];
} else {
fatal_error(
fmt::format("<width> must be length 2 in slice plot {}", id_));
}
} else if (PlotType::voxel == type_) {
if (pl_width.size() == 3) {
pl_width = get_node_array<double>(plot_node, "width");
width_ = pl_width;
} else {
fatal_error(
fmt::format("<width> must be length 3 in voxel plot {}", id_));
}
}
}
void Plot::set_universe(pugi::xml_node plot_node)
{
// Copy plot universe level
if (check_for_node(plot_node, "level")) {
level_ = std::stoi(get_node_value(plot_node, "level"));
if (level_ < 0) {
fatal_error(fmt::format("Bad universe level in plot {}", id_));
}
} else {
level_ = PLOT_LEVEL_LOWEST;
}
}
void Plot::set_default_colors(pugi::xml_node plot_node)
{
// Copy plot color type and initialize all colors randomly
std::string pl_color_by = "cell";
if (check_for_node(plot_node, "color_by")) {
pl_color_by = get_node_value(plot_node, "color_by", true);
}
if ("cell" == pl_color_by) {
color_by_ = PlotColorBy::cells;
colors_.resize(model::cells.size());
} else if ("material" == pl_color_by) {
color_by_ = PlotColorBy::mats;
colors_.resize(model::materials.size());
} else {
fatal_error(fmt::format(
"Unsupported plot color type '{}' in plot {}", pl_color_by, id_));
}
for (auto& c : colors_) {
c = random_color();
// make sure we don't interfere with some default colors
while (c == RED || c == WHITE) {
c = random_color();
}
}
}
void Plot::set_user_colors(pugi::xml_node plot_node)
{
if (!plot_node.select_nodes("color").empty() && PlotType::voxel == type_) {
if (mpi::master) {
warning(
fmt::format("Color specifications ignored in voxel plot {}", id_));
}
}
for (auto cn : plot_node.children("color")) {
// Make sure 3 values are specified for RGB
vector<int> user_rgb = get_node_array<int>(cn, "rgb");
if (user_rgb.size() != 3) {
fatal_error(fmt::format("Bad RGB in plot {}", id_));
}
// Ensure that there is an id for this color specification
int col_id;
if (check_for_node(cn, "id")) {
col_id = std::stoi(get_node_value(cn, "id"));
} else {
fatal_error(
fmt::format("Must specify id for color specification in plot {}", id_));
}
// Add RGB
if (PlotColorBy::cells == color_by_) {
if (model::cell_map.find(col_id) != model::cell_map.end()) {
col_id = model::cell_map[col_id];
colors_[col_id] = user_rgb;
} else {
warning(fmt::format(
"Could not find cell {} specified in plot {}", col_id, id_));
}
} else if (PlotColorBy::mats == color_by_) {
if (model::material_map.find(col_id) != model::material_map.end()) {
col_id = model::material_map[col_id];
colors_[col_id] = user_rgb;
} else {
warning(fmt::format(
"Could not find material {} specified in plot {}", col_id, id_));
}
}
} // color node loop
}
void Plot::set_meshlines(pugi::xml_node plot_node)
{
// Deal with meshlines
pugi::xpath_node_set mesh_line_nodes = plot_node.select_nodes("meshlines");
if (!mesh_line_nodes.empty()) {
if (PlotType::voxel == type_) {
warning(fmt::format("Meshlines ignored in voxel plot {}", id_));
}
if (mesh_line_nodes.size() == 1) {
// Get first meshline node
pugi::xml_node meshlines_node = mesh_line_nodes[0].node();
// Check mesh type
std::string meshtype;
if (check_for_node(meshlines_node, "meshtype")) {
meshtype = get_node_value(meshlines_node, "meshtype");
} else {
fatal_error(fmt::format(
"Must specify a meshtype for meshlines specification in plot {}",
id_));
}
// Ensure that there is a linewidth for this meshlines specification
std::string meshline_width;
if (check_for_node(meshlines_node, "linewidth")) {
meshline_width = get_node_value(meshlines_node, "linewidth");
meshlines_width_ = std::stoi(meshline_width);
} else {
fatal_error(fmt::format(
"Must specify a linewidth for meshlines specification in plot {}",
id_));
}
// Check for color
if (check_for_node(meshlines_node, "color")) {
// Check and make sure 3 values are specified for RGB
vector<int> ml_rgb = get_node_array<int>(meshlines_node, "color");
if (ml_rgb.size() != 3) {
fatal_error(
fmt::format("Bad RGB for meshlines color in plot {}", id_));
}
meshlines_color_ = ml_rgb;
}
// Set mesh based on type
if ("ufs" == meshtype) {
if (!simulation::ufs_mesh) {
fatal_error(fmt::format("No UFS mesh for meshlines on plot {}", id_));
} else {
for (int i = 0; i < model::meshes.size(); ++i) {
if (const auto* m =
dynamic_cast<const RegularMesh*>(model::meshes[i].get())) {
if (m == simulation::ufs_mesh) {
index_meshlines_mesh_ = i;
}
}
}
if (index_meshlines_mesh_ == -1)
fatal_error("Could not find the UFS mesh for meshlines plot");
}
} else if ("entropy" == meshtype) {
if (!simulation::entropy_mesh) {
fatal_error(
fmt::format("No entropy mesh for meshlines on plot {}", id_));
} else {
for (int i = 0; i < model::meshes.size(); ++i) {
if (const auto* m =
dynamic_cast<const RegularMesh*>(model::meshes[i].get())) {
if (m == simulation::entropy_mesh) {
index_meshlines_mesh_ = i;
}
}
}
if (index_meshlines_mesh_ == -1)
fatal_error("Could not find the entropy mesh for meshlines plot");
}
} else if ("tally" == meshtype) {
// Ensure that there is a mesh id if the type is tally
int tally_mesh_id;
if (check_for_node(meshlines_node, "id")) {
tally_mesh_id = std::stoi(get_node_value(meshlines_node, "id"));
} else {
std::stringstream err_msg;
fatal_error(fmt::format("Must specify a mesh id for meshlines tally "
"mesh specification in plot {}",
id_));
}
// find the tally index
int idx;
int err = openmc_get_mesh_index(tally_mesh_id, &idx);
if (err != 0) {
fatal_error(fmt::format("Could not find mesh {} specified in "
"meshlines for plot {}",
tally_mesh_id, id_));
}
index_meshlines_mesh_ = idx;
} else {
fatal_error(fmt::format("Invalid type for meshlines on plot {}", id_));
}
} else {
fatal_error(fmt::format("Mutliple meshlines specified in plot {}", id_));
}
}
}
void Plot::set_mask(pugi::xml_node plot_node)
{
// Deal with masks
pugi::xpath_node_set mask_nodes = plot_node.select_nodes("mask");
if (!mask_nodes.empty()) {
if (PlotType::voxel == type_) {
if (mpi::master) {
warning(fmt::format("Mask ignored in voxel plot {}", id_));
}
}
if (mask_nodes.size() == 1) {
// Get pointer to mask
pugi::xml_node mask_node = mask_nodes[0].node();
// Determine how many components there are and allocate
vector<int> iarray = get_node_array<int>(mask_node, "components");
if (iarray.size() == 0) {
fatal_error(
fmt::format("Missing <components> in mask of plot {}", id_));
}
// First we need to change the user-specified identifiers to indices
// in the cell and material arrays
for (auto& col_id : iarray) {
if (PlotColorBy::cells == color_by_) {
if (model::cell_map.find(col_id) != model::cell_map.end()) {
col_id = model::cell_map[col_id];
} else {
fatal_error(fmt::format("Could not find cell {} specified in the "
"mask in plot {}",
col_id, id_));
}
} else if (PlotColorBy::mats == color_by_) {
if (model::material_map.find(col_id) != model::material_map.end()) {
col_id = model::material_map[col_id];
} else {
fatal_error(fmt::format("Could not find material {} specified in "
"the mask in plot {}",
col_id, id_));
}
}
}
// Alter colors based on mask information
for (int j = 0; j < colors_.size(); j++) {
if (std::find(iarray.begin(), iarray.end(), j) == iarray.end()) {
if (check_for_node(mask_node, "background")) {
vector<int> bg_rgb = get_node_array<int>(mask_node, "background");
colors_[j] = bg_rgb;
} else {
colors_[j] = WHITE;
}
}
}
} else {
fatal_error(fmt::format("Mutliple masks specified in plot {}", id_));
}
}
}
void Plot::set_overlap_color(pugi::xml_node plot_node)
{
color_overlaps_ = false;
if (check_for_node(plot_node, "show_overlaps")) {
color_overlaps_ = get_node_value_bool(plot_node, "show_overlaps");
// check for custom overlap color
if (check_for_node(plot_node, "overlap_color")) {
if (!color_overlaps_) {
warning(fmt::format(
"Overlap color specified in plot {} but overlaps won't be shown.",
id_));
}
vector<int> olap_clr = get_node_array<int>(plot_node, "overlap_color");
if (olap_clr.size() == 3) {
overlap_color_ = olap_clr;
} else {
fatal_error(fmt::format("Bad overlap RGB in plot {}", id_));
}
}
}
// make sure we allocate the vector for counting overlap checks if
// they're going to be plotted
if (color_overlaps_ && settings::run_mode == RunMode::PLOTTING) {
settings::check_overlaps = true;
model::overlap_check_count.resize(model::cells.size(), 0);
}
}
Plot::Plot(pugi::xml_node plot_node)
: index_meshlines_mesh_ {-1}, overlap_color_ {RED}
{
set_id(plot_node);
set_type(plot_node);
set_output_path(plot_node);
set_bg_color(plot_node);
set_basis(plot_node);
set_origin(plot_node);
set_width(plot_node);
set_universe(plot_node);
set_default_colors(plot_node);
set_user_colors(plot_node);
set_meshlines(plot_node);
set_mask(plot_node);
set_overlap_color(plot_node);
} // End Plot constructor
//==============================================================================
// OUTPUT_PPM writes out a previously generated image to a PPM file
//==============================================================================
void output_ppm(Plot const& pl, const ImageData& data)
{
// Open PPM file for writing
std::string fname = pl.path_plot_;
fname = strtrim(fname);
std::ofstream of;
of.open(fname);
// Write header
of << "P6\n";
of << pl.pixels_[0] << " " << pl.pixels_[1] << "\n";
of << "255\n";
of.close();
of.open(fname, std::ios::binary | std::ios::app);
// Write color for each pixel
for (int y = 0; y < pl.pixels_[1]; y++) {
for (int x = 0; x < pl.pixels_[0]; x++) {
RGBColor rgb = data(x, y);
of << rgb.red << rgb.green << rgb.blue;
}
}
of << "\n";
}
//==============================================================================
// OUTPUT_PNG writes out a previously generated image to a PNG file
//==============================================================================
#ifdef USE_LIBPNG
void output_png(Plot const& pl, const ImageData& data)
{
// Open PNG file for writing
std::string fname = pl.path_plot_;
fname = strtrim(fname);
auto fp = std::fopen(fname.c_str(), "wb");
// Initialize write and info structures
auto png_ptr =
png_create_write_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr);
auto info_ptr = png_create_info_struct(png_ptr);
// Setup exception handling
if (setjmp(png_jmpbuf(png_ptr)))
fatal_error("Error during png creation");
png_init_io(png_ptr, fp);
// Write header (8 bit colour depth)
int width = pl.pixels_[0];
int height = pl.pixels_[1];
png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGB,
PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
png_write_info(png_ptr, info_ptr);
// Allocate memory for one row (3 bytes per pixel - RGB)
std::vector<png_byte> row(3 * width);
// Write color for each pixel
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
RGBColor rgb = data(x, y);
row[3 * x] = rgb.red;
row[3 * x + 1] = rgb.green;
row[3 * x + 2] = rgb.blue;
}
png_write_row(png_ptr, row.data());
}
// End write
png_write_end(png_ptr, nullptr);
// Clean up data structures
std::fclose(fp);
png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
png_destroy_write_struct(&png_ptr, &info_ptr);
}
#endif
//==============================================================================
// DRAW_MESH_LINES draws mesh line boundaries on an image
//==============================================================================
void draw_mesh_lines(Plot const& pl, ImageData& data)
{
RGBColor rgb;
rgb = pl.meshlines_color_;
int ax1, ax2;
switch (pl.basis_) {
case PlotBasis::xy:
ax1 = 0;
ax2 = 1;
break;
case PlotBasis::xz:
ax1 = 0;
ax2 = 2;
break;
case PlotBasis::yz:
ax1 = 1;
ax2 = 2;
break;
default:
UNREACHABLE();
}
Position ll_plot {pl.origin_};
Position ur_plot {pl.origin_};
ll_plot[ax1] -= pl.width_[0] / 2.;
ll_plot[ax2] -= pl.width_[1] / 2.;
ur_plot[ax1] += pl.width_[0] / 2.;
ur_plot[ax2] += pl.width_[1] / 2.;
Position width = ur_plot - ll_plot;
// Find the (axis-aligned) lines of the mesh that intersect this plot.
auto axis_lines =
model::meshes[pl.index_meshlines_mesh_]->plot(ll_plot, ur_plot);
// Find the bounds along the second axis (accounting for low-D meshes).
int ax2_min, ax2_max;
if (axis_lines.second.size() > 0) {
double frac = (axis_lines.second.back() - ll_plot[ax2]) / width[ax2];
ax2_min = (1.0 - frac) * pl.pixels_[1];
if (ax2_min < 0)
ax2_min = 0;
frac = (axis_lines.second.front() - ll_plot[ax2]) / width[ax2];
ax2_max = (1.0 - frac) * pl.pixels_[1];
if (ax2_max > pl.pixels_[1])
ax2_max = pl.pixels_[1];
} else {
ax2_min = 0;
ax2_max = pl.pixels_[1];
}
// Iterate across the first axis and draw lines.
for (auto ax1_val : axis_lines.first) {
double frac = (ax1_val - ll_plot[ax1]) / width[ax1];
int ax1_ind = frac * pl.pixels_[0];
for (int ax2_ind = ax2_min; ax2_ind < ax2_max; ++ax2_ind) {
for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
if (ax1_ind + plus >= 0 && ax1_ind + plus < pl.pixels_[0])
data(ax1_ind + plus, ax2_ind) = rgb;
if (ax1_ind - plus >= 0 && ax1_ind - plus < pl.pixels_[0])
data(ax1_ind - plus, ax2_ind) = rgb;
}
}
}
// Find the bounds along the first axis.
int ax1_min, ax1_max;
if (axis_lines.first.size() > 0) {
double frac = (axis_lines.first.front() - ll_plot[ax1]) / width[ax1];
ax1_min = frac * pl.pixels_[0];
if (ax1_min < 0)
ax1_min = 0;
frac = (axis_lines.first.back() - ll_plot[ax1]) / width[ax1];
ax1_max = frac * pl.pixels_[0];
if (ax1_max > pl.pixels_[0])
ax1_max = pl.pixels_[0];
} else {
ax1_min = 0;
ax1_max = pl.pixels_[0];
}
// Iterate across the second axis and draw lines.
for (auto ax2_val : axis_lines.second) {
double frac = (ax2_val - ll_plot[ax2]) / width[ax2];
int ax2_ind = (1.0 - frac) * pl.pixels_[1];
for (int ax1_ind = ax1_min; ax1_ind < ax1_max; ++ax1_ind) {
for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
if (ax2_ind + plus >= 0 && ax2_ind + plus < pl.pixels_[1])
data(ax1_ind, ax2_ind + plus) = rgb;
if (ax2_ind - plus >= 0 && ax2_ind - plus < pl.pixels_[1])
data(ax1_ind, ax2_ind - plus) = rgb;
}
}
}
}
//==============================================================================
// CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D
// geometry visualization. It works the same way as create_image by dragging a
// particle across the geometry for the specified number of voxels. The first 3
// int's in the binary are the number of x, y, and z voxels. The next 3
// double's are the widths of the voxels in the x, y, and z directions. The
// next 3 double's are the x, y, and z coordinates of the lower left
// point. Finally the binary is filled with entries of four int's each. Each
// 'row' in the binary contains four int's: 3 for x,y,z position and 1 for
// cell or material id. For 1 million voxels this produces a file of
// approximately 15MB.
// =============================================================================
void create_voxel(Plot const& pl)
{
// compute voxel widths in each direction
array<double, 3> vox;
vox[0] = pl.width_[0] / (double)pl.pixels_[0];
vox[1] = pl.width_[1] / (double)pl.pixels_[1];
vox[2] = pl.width_[2] / (double)pl.pixels_[2];
// initial particle position
Position ll = pl.origin_ - pl.width_ / 2.;
// Open binary plot file for writing
std::ofstream of;
std::string fname = std::string(pl.path_plot_);
fname = strtrim(fname);
hid_t file_id = file_open(fname, 'w');
// write header info
write_attribute(file_id, "filetype", "voxel");
write_attribute(file_id, "version", VERSION_VOXEL);
write_attribute(file_id, "openmc_version", VERSION);
#ifdef GIT_SHA1
write_attribute(file_id, "git_sha1", GIT_SHA1);
#endif
// Write current date and time
write_attribute(file_id, "date_and_time", time_stamp().c_str());
array<int, 3> pixels;
std::copy(pl.pixels_.begin(), pl.pixels_.end(), pixels.begin());
write_attribute(file_id, "num_voxels", pixels);
write_attribute(file_id, "voxel_width", vox);
write_attribute(file_id, "lower_left", ll);
// Create dataset for voxel data -- note that the dimensions are reversed
// since we want the order in the file to be z, y, x
hsize_t dims[3];
dims[0] = pl.pixels_[2];
dims[1] = pl.pixels_[1];
dims[2] = pl.pixels_[0];
hid_t dspace, dset, memspace;
voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
PlotBase pltbase;
pltbase.width_ = pl.width_;
pltbase.origin_ = pl.origin_;
pltbase.basis_ = PlotBasis::xy;
pltbase.pixels_ = pl.pixels_;
pltbase.level_ = -1; // all universes for voxel files
pltbase.color_overlaps_ = pl.color_overlaps_;
ProgressBar pb;
for (int z = 0; z < pl.pixels_[2]; z++) {
// update progress bar
pb.set_value(100. * (double)z / (double)(pl.pixels_[2] - 1));
// update z coordinate
pltbase.origin_.z = ll.z + z * vox[2];
// generate ids using plotbase
IdData ids = pltbase.get_map<IdData>();
// select only cell/material ID data and flip the y-axis
int idx = pl.color_by_ == PlotColorBy::cells ? 0 : 2;
xt::xtensor<int32_t, 2> data_slice =
xt::view(ids.data_, xt::all(), xt::all(), idx);
xt::xtensor<int32_t, 2> data_flipped = xt::flip(data_slice, 0);
// Write to HDF5 dataset
voxel_write_slice(z, dspace, dset, memspace, data_flipped.data());
}
voxel_finalize(dspace, dset, memspace);
file_close(file_id);
}
void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
hid_t* memspace)
{
// Create dataspace/dataset for voxel data
*dspace = H5Screate_simple(3, dims, nullptr);
*dset = H5Dcreate(file_id, "data", H5T_NATIVE_INT, *dspace, H5P_DEFAULT,
H5P_DEFAULT, H5P_DEFAULT);
// Create dataspace for a slice of the voxel
hsize_t dims_slice[2] {dims[1], dims[2]};
*memspace = H5Screate_simple(2, dims_slice, nullptr);
// Select hyperslab in dataspace
hsize_t start[3] {0, 0, 0};
hsize_t count[3] {1, dims[1], dims[2]};
H5Sselect_hyperslab(*dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
}
void voxel_write_slice(
int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf)
{
hssize_t offset[3] {x, 0, 0};
H5Soffset_simple(dspace, offset);
H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf);
}
void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
{
H5Dclose(dset);
H5Sclose(dspace);
H5Sclose(memspace);
}
RGBColor random_color(void)
{
return {int(prn(&model::plotter_seed) * 255),
int(prn(&model::plotter_seed) * 255), int(prn(&model::plotter_seed) * 255)};
}
extern "C" int openmc_id_map(const void* plot, int32_t* data_out)
{
auto plt = reinterpret_cast<const PlotBase*>(plot);
if (!plt) {
set_errmsg("Invalid slice pointer passed to openmc_id_map");
return OPENMC_E_INVALID_ARGUMENT;
}
if (plt->color_overlaps_ && model::overlap_check_count.size() == 0) {
model::overlap_check_count.resize(model::cells.size());
}
auto ids = plt->get_map<IdData>();
// write id data to array
std::copy(ids.data_.begin(), ids.data_.end(), data_out);
return 0;
}
extern "C" int openmc_property_map(const void* plot, double* data_out)
{
auto plt = reinterpret_cast<const PlotBase*>(plot);
if (!plt) {
set_errmsg("Invalid slice pointer passed to openmc_id_map");
return OPENMC_E_INVALID_ARGUMENT;
}
if (plt->color_overlaps_ && model::overlap_check_count.size() == 0) {
model::overlap_check_count.resize(model::cells.size());
}
auto props = plt->get_map<PropertyData>();
// write id data to array
std::copy(props.data_.begin(), props.data_.end(), data_out);
return 0;
}
} // namespace openmc
Reference in a new issue View git blame Copy permalink