From 7a49a29f38c6537a1346754ec66cf7d9a0aaa450 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 29 May 2019 15:29:54 -0400 Subject: [PATCH] Refactor and fix RegularMesh meshline plots --- include/openmc/mesh.h | 9 ++ src/mesh.cpp | 39 +++++ src/plot.cpp | 141 ++++++++++--------- tests/regression_tests/plot/plots.xml | 3 + tests/regression_tests/plot/results_true.dat | 2 +- 5 files changed, 127 insertions(+), 67 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 12158339d..8233d6779 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -85,6 +85,9 @@ public: //! Get the number of mesh cell surfaces. virtual int n_surface_bins() const = 0; + virtual std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const = 0; + //! Write mesh data to an HDF5 group // //! \param[in] group HDF5 group @@ -127,6 +130,9 @@ public: int n_surface_bins() const override; + std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const override; + void to_hdf5(hid_t group) const override; // New methods @@ -187,6 +193,9 @@ public: int n_surface_bins() const override; + std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const override; + void to_hdf5(hid_t group) const override; // New methods diff --git a/src/mesh.cpp b/src/mesh.cpp index 2cb91b98f..f398d641d 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -729,6 +729,40 @@ void RegularMesh::surface_bins_crossed(const Particle* p, } } +std::pair, std::vector> +RegularMesh::plot(Position plot_ll, Position plot_ur) const +{ + std::array axes {-1, -1}; + if (plot_ur.z == plot_ll.z) { + axes[0] = 0; + if (n_dimension_ > 1) axes[1] = 1; + } else if (plot_ur.y == plot_ll.y) { + axes[0] = 0; + if (n_dimension_ > 2) axes[1] = 2; + } else if (plot_ur.x == plot_ll.x) { + if (n_dimension_ > 1) axes[0] = 1; + if (n_dimension_ > 2) axes[1] = 2; + } else { + fatal_error("Can only plot mesh lines on an axis-aligned plot"); + } + + std::array, 2> axis_lines; + 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]}; + + double coord = lower_left_[axis]; + for (int i = 0; i < shape_[axis] + 1; ++i) { + if (coord >= plot_ll[axis] && coord <= plot_ur[axis]) + lines.push_back(coord); + coord += width_[axis]; + } + } + + return {axis_lines[0], axis_lines[1]}; +} + void RegularMesh::to_hdf5(hid_t group) const { hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); @@ -1125,6 +1159,11 @@ int RectilinearMesh::n_surface_bins() const return 4 * n_dimension_ * n_bins(); } +std::pair, std::vector> +RectilinearMesh::plot(Position plot_ll, Position plot_ur) const +{ +} + void RectilinearMesh::to_hdf5(hid_t group) const { hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); diff --git a/src/plot.cpp b/src/plot.cpp index 43e7ed69c..bfb6b335f 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -698,19 +698,19 @@ void draw_mesh_lines(Plot pl, ImageData& data) RGBColor rgb; rgb = pl.meshlines_color_; - int outer, inner; + int ax1, ax2; switch(pl.basis_) { case PlotBasis::xy : - outer = 0; - inner = 1; + ax1 = 0; + ax2 = 1; break; case PlotBasis::xz : - outer = 0; - inner = 2; + ax1 = 0; + ax2 = 2; break; case PlotBasis::yz : - outer = 1; - inner = 2; + ax1 = 1; + ax2 = 2; break; default: UNREACHABLE(); @@ -719,70 +719,79 @@ void draw_mesh_lines(Plot pl, ImageData& data) Position ll_plot {pl.origin_}; Position ur_plot {pl.origin_}; - ll_plot[outer] -= pl.width_[0] / 2.; - ll_plot[inner] -= pl.width_[1] / 2.; - ur_plot[outer] += pl.width_[0] / 2.; - ur_plot[inner] += pl.width_[1] / 2.; + 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; - auto& m = *dynamic_cast(model::meshes[pl.index_meshlines_mesh_].get()); + // 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); - int ijk_ll[3], ijk_ur[3]; - bool in_mesh; - m.get_indices(ll_plot, &(ijk_ll[0]), &in_mesh); - m.get_indices(ur_plot, &(ijk_ur[0]), &in_mesh); - - // Fortran/C++ index correction - ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++; - - Position r_ll, r_ur; - // sweep through all meshbins on this plane and draw borders - for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) { - for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) { - // check if we're in the mesh for this ijk - if (i > 0 && i <= m.shape_[outer] && j >0 && j <= m.shape_[inner] ) { - int outrange[3], inrange[3]; - // get xyz's of lower left and upper right of this mesh cell - r_ll[outer] = m.lower_left_[outer] + m.width_[outer] * (i - 1); - r_ll[inner] = m.lower_left_[inner] + m.width_[inner] * (j - 1); - r_ur[outer] = m.lower_left_[outer] + m.width_[outer] * i; - r_ur[inner] = m.lower_left_[inner] + m.width_[inner] * j; - - // map the xyz ranges to pixel ranges - double frac = (r_ll[outer] - ll_plot[outer]) / width[outer]; - outrange[0] = int(frac * double(pl.pixels_[0])); - frac = (r_ur[outer] - ll_plot[outer]) / width[outer]; - outrange[1] = int(frac * double(pl.pixels_[0])); - - frac = (r_ur[inner] - ll_plot[inner]) / width[inner]; - inrange[0] = int((1. - frac) * (double)pl.pixels_[1]); - frac = (r_ll[inner] - ll_plot[inner]) / width[inner]; - inrange[1] = int((1. - frac) * (double)pl.pixels_[1]); - - // draw lines - for (int out_ = outrange[0]; out_ <= outrange[1]; out_++) { - for (int plus = 0; plus <= pl.meshlines_width_; plus++) { - data(out_, inrange[0] + plus) = rgb; - data(out_, inrange[1] + plus) = rgb; - data(out_, inrange[0] - plus) = rgb; - data(out_, inrange[1] - plus) = rgb; - } - } - - for (int in_ = inrange[0]; in_ <= inrange[1]; in_++) { - for (int plus = 0; plus <= pl.meshlines_width_; plus++) { - data(outrange[0] + plus, in_) = rgb; - data(outrange[1] + plus, in_) = rgb; - data(outrange[0] - plus, in_) = rgb; - data(outrange[1] - plus, in_) = rgb; - } - } - - } // end if(in mesh) + // 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.front() - ll_plot[ax2]) / width[ax2]; + ax2_min = frac * pl.pixels_[1]; + if (ax2_min < 0) ax2_min = 0; + frac = (axis_lines.second.back() - ll_plot[ax2]) / width[ax2]; + ax2_max = frac * pl.pixels_[1]; + if (ax2_max > pl.pixels_[1]) ax2_max = pl.pixels_[1]; + } else { + ax2_min = 0; + ax2_max = pl.pixels_[1]; } - } // end outer loops -} // end draw_mesh_lines + + // 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]; + } + + // 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 = 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 diff --git a/tests/regression_tests/plot/plots.xml b/tests/regression_tests/plot/plots.xml index ecfe69125..075124a8f 100644 --- a/tests/regression_tests/plot/plots.xml +++ b/tests/regression_tests/plot/plots.xml @@ -6,7 +6,10 @@ 25 25 200 200 + + diff --git a/tests/regression_tests/plot/results_true.dat b/tests/regression_tests/plot/results_true.dat index 18251e8d0..128aa9bbe 100644 --- a/tests/regression_tests/plot/results_true.dat +++ b/tests/regression_tests/plot/results_true.dat @@ -1 +1 @@ -6b3eb36488b995b42233e6b7c0164cf6c92a1823b3c91985b97fd076f86c0aad93fbdb74e70026f5f12c61a6aee52a09588171a7fd839511ee7d9efc3952beb2 \ No newline at end of file +ae51cbf147d69804b03f9d88c03d98533665a9822438d1d8b70b40928a1db35d20e18d12ee697b5e1e217799d6e8c352524e8d2227328e719e85c4e6d96f825c \ No newline at end of file