From b804ddc952696ea00b75ef4748f36650f4ca7370 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 6 Jun 2019 15:19:15 -0400 Subject: [PATCH] Address comments on #1246 --- include/openmc/mesh.h | 14 +++--- openmc/mesh.py | 12 ++--- src/mesh.cpp | 16 ++++++- src/plot.cpp | 102 ++++++++++++++++++++---------------------- src/settings.cpp | 15 ++++--- 5 files changed, 83 insertions(+), 76 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 0e8d75ea3..6d59826aa 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -10,7 +10,7 @@ #include "hdf5.h" #include "pugixml.hpp" -#include "xtensor/xarray.hpp" +#include "xtensor/xtensor.hpp" #include "openmc/particle.h" #include "openmc/position.h" @@ -105,6 +105,8 @@ public: int id_ {-1}; //!< User-specified ID int n_dimension_; //!< Number of dimensions + xt::xtensor lower_left_; //!< Lower-left coordinates of mesh + xt::xtensor upper_right_; //!< Upper-right coordinates of mesh }; //============================================================================== @@ -158,16 +160,14 @@ public: //! \param[in] bank Array of bank sites //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin - xt::xarray count_sites(const std::vector& bank, + xt::xtensor count_sites(const std::vector& bank, bool* outside) const; // Data members double volume_frac_; //!< Volume fraction of each mesh element - xt::xarray shape_; //!< Number of mesh elements in each dimension - xt::xarray lower_left_; //!< Lower-left coordinates of mesh - xt::xarray upper_right_; //!< Upper-right coordinates of mesh - xt::xarray width_; //!< Width of each mesh element + xt::xtensor shape_; //!< Number of mesh elements in each dimension + xt::xtensor width_; //!< Width of each mesh element private: bool intersects_1d(Position& r0, Position r1, int* ijk) const; @@ -218,7 +218,7 @@ public: // Data members - xt::xarray shape_; //!< Number of mesh elements in each dimension + xt::xtensor shape_; //!< Number of mesh elements in each dimension private: std::vector> grid_; diff --git a/openmc/mesh.py b/openmc/mesh.py index ef3212b90..5b6b8ad3e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -104,9 +104,9 @@ class Mesh(MeshBase): are given, it is assumed that the mesh is an x-y mesh. width : Iterable of float The width of mesh cells in each direction. - indices : list of tuple - A list of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, - 1), ...] + indices : Iterable of tuple + An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), + (2, 1, 1), ...] """ @@ -427,9 +427,9 @@ class RectilinearMesh(MeshBase): Mesh boundary points along the y-axis. z_grid : Iterable of float Mesh boundary points along the z-axis. - indices : list of tuple - A list of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, - 1), ...] + indices : Iterable of tuple + An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), + (2, 1, 1), ...] """ diff --git a/src/mesh.cpp b/src/mesh.cpp index 870a58e6c..315381706 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -752,7 +752,7 @@ RegularMesh::plot(Position plot_ll, Position plot_ur) const for (int i_ax = 0; i_ax < 2; ++i_ax) { int axis = axes[i_ax]; if (axis == -1) continue; - std::vector& lines {axis_lines[i_ax]}; + auto& lines {axis_lines[i_ax]}; double coord = lower_left_[axis]; for (int i = 0; i < shape_[axis] + 1; ++i) { @@ -778,7 +778,7 @@ void RegularMesh::to_hdf5(hid_t group) const close_group(mesh_group); } -xt::xarray +xt::xtensor RegularMesh::count_sites(const std::vector& bank, bool* outside) const { @@ -846,6 +846,18 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) shape_ = {static_cast(grid_[0].size()) - 1, static_cast(grid_[1].size()) - 1, static_cast(grid_[2].size()) - 1}; + + for (const auto& g : grid_) { + if (g.size() < 2) fatal_error("x-, y-, and z- grids for rectilinear meshes " + "must each have at least 2 points"); + for (int i = 1; i < g.size(); ++i) { + if (g[i] <= g[i-1]) fatal_error("Values in for x-, y-, and z- grids for " + "rectilinear meshes must be sorted and unique."); + } + } + + lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()}; + upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()}; } void RectilinearMesh::bins_crossed(const Particle* p, std::vector& bins, diff --git a/src/plot.cpp b/src/plot.cpp index 1ad4f8317..6bc82e164 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -727,67 +727,61 @@ void draw_mesh_lines(Plot pl, ImageData& data) Position width = ur_plot - ll_plot; // Find the (axis-aligned) lines of the mesh that intersect this plot. - std::pair, std::vector> axis_lines; - axis_lines = model::meshes[pl.index_meshlines_mesh_]->plot(ll_plot, ur_plot); + auto axis_lines = model::meshes[pl.index_meshlines_mesh_] + ->plot(ll_plot, ur_plot); - // Draw lines perpendicular to the first axis. - { - // 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]; - } + // 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; - } + // 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; } } } - // Draw lines perpendicular to the second axis. - { - // 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]; - } + // 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; - } + // 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; } } } diff --git a/src/settings.cpp b/src/settings.cpp index 63827c3f4..ac31d2670 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -478,7 +478,7 @@ void read_settings_xml() read_meshes(root); // Shannon Entropy mesh - int32_t i_entropy_mesh = -1; + int32_t index_entropy_mesh = -1; if (check_for_node(root, "entropy_mesh")) { int temp = std::stoi(get_node_value(root, "entropy_mesh")); if (model::mesh_map.find(temp) == model::mesh_map.end()) { @@ -486,7 +486,7 @@ void read_settings_xml() msg << "Mesh " << temp << " specified for Shannon entropy does not exist."; fatal_error(msg); } - i_entropy_mesh = model::mesh_map.at(temp); + index_entropy_mesh = model::mesh_map.at(temp); } else if (check_for_node(root, "entropy")) { warning("Specifying a Shannon entropy mesh via the element " @@ -498,15 +498,16 @@ void read_settings_xml() model::meshes.push_back(std::make_unique(node_entropy)); // Set entropy mesh index - i_entropy_mesh = model::meshes.size() - 1; + index_entropy_mesh = model::meshes.size() - 1; // Assign ID and set mapping model::meshes.back()->id_ = 10000; - model::mesh_map[10000] = i_entropy_mesh; + model::mesh_map[10000] = index_entropy_mesh; } - if (i_entropy_mesh >= 0) { - auto* m = dynamic_cast(model::meshes[i_entropy_mesh].get()); + if (index_entropy_mesh >= 0) { + auto* m = dynamic_cast( + model::meshes[index_entropy_mesh].get()); if (!m) fatal_error("Only regular meshes can be used as an entropy mesh"); simulation::entropy_mesh = m; @@ -552,7 +553,7 @@ void read_settings_xml() // Assign ID and set mapping model::meshes.back()->id_ = 10001; - model::mesh_map[10001] = i_entropy_mesh; + model::mesh_map[10001] = index_entropy_mesh; } if (i_ufs_mesh >= 0) {