From 7f8e993e93673b58ecfd6ce961a7c93534413add Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 29 May 2019 15:48:58 -0400 Subject: [PATCH] Implement rectilinear meshline plots --- include/openmc/mesh.h | 8 ++++++ src/mesh.cpp | 27 ++++++++++++++++++++ src/plot.cpp | 10 ++++---- tests/regression_tests/plot/plots.xml | 4 +-- tests/regression_tests/plot/results_true.dat | 2 +- tests/regression_tests/plot/tallies.xml | 20 +++++++++++++++ 6 files changed, 62 insertions(+), 9 deletions(-) create mode 100644 tests/regression_tests/plot/tallies.xml diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 8233d6779..0e8d75ea3 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -85,6 +85,14 @@ public: //! Get the number of mesh cell surfaces. virtual int n_surface_bins() const = 0; + //! Find the mesh lines that intersect an axis-aligned slice plot + // + //! \param[in] plot_ll The lower-left coordinates of the slice plot. + //! \param[in] plot_ur The upper-right coordinates of the slice plot. + //! \return A pair of vectors indicating where the mesh lines lie along each + //! of the plot's axes. For example an xy-slice plot will get back a vector + //! of x-coordinates and another of y-coordinates. These vectors may be + //! empty for low-dimensional meshes. virtual std::pair, std::vector> plot(Position plot_ll, Position plot_ur) const = 0; diff --git a/src/mesh.cpp b/src/mesh.cpp index f398d641d..870a58e6c 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -732,6 +732,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::pair, std::vector> RegularMesh::plot(Position plot_ll, Position plot_ur) const { + // Figure out which axes lie in the plane of the plot. std::array axes {-1, -1}; if (plot_ur.z == plot_ll.z) { axes[0] = 0; @@ -746,6 +747,7 @@ RegularMesh::plot(Position plot_ll, Position plot_ur) const fatal_error("Can only plot mesh lines on an axis-aligned plot"); } + // Get the coordinates of the mesh lines along both of the axes. std::array, 2> axis_lines; for (int i_ax = 0; i_ax < 2; ++i_ax) { int axis = axes[i_ax]; @@ -1162,6 +1164,31 @@ int RectilinearMesh::n_surface_bins() const std::pair, std::vector> RectilinearMesh::plot(Position plot_ll, Position plot_ur) const { + // Figure out which axes lie in the plane of the plot. + std::array axes {-1, -1}; + if (plot_ur.z == plot_ll.z) { + axes = {0, 1}; + } else if (plot_ur.y == plot_ll.y) { + axes = {0, 2}; + } else if (plot_ur.x == plot_ll.x) { + axes = {1, 2}; + } else { + fatal_error("Can only plot mesh lines on an axis-aligned plot"); + } + + // Get the coordinates of the mesh lines along both of the axes. + std::array, 2> axis_lines; + for (int i_ax = 0; i_ax < 2; ++i_ax) { + int axis = axes[i_ax]; + std::vector& lines {axis_lines[i_ax]}; + + for (auto coord : grid_[axis]) { + if (coord >= plot_ll[axis] && coord <= plot_ur[axis]) + lines.push_back(coord); + } + } + + return {axis_lines[0], axis_lines[1]}; } void RectilinearMesh::to_hdf5(hid_t group) const diff --git a/src/plot.cpp b/src/plot.cpp index bfb6b335f..1ad4f8317 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -735,11 +735,11 @@ void draw_mesh_lines(Plot pl, ImageData& data) // 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]; + 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.back() - ll_plot[ax2]) / width[ax2]; - ax2_max = frac * pl.pixels_[1]; + 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; @@ -780,7 +780,7 @@ void draw_mesh_lines(Plot pl, ImageData& data) // 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]; + 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]) diff --git a/tests/regression_tests/plot/plots.xml b/tests/regression_tests/plot/plots.xml index 075124a8f..ce63da144 100644 --- a/tests/regression_tests/plot/plots.xml +++ b/tests/regression_tests/plot/plots.xml @@ -6,10 +6,7 @@ 25 25 200 200 - - @@ -17,6 +14,7 @@ 25 25 200 200 + diff --git a/tests/regression_tests/plot/results_true.dat b/tests/regression_tests/plot/results_true.dat index 128aa9bbe..db6beb9cd 100644 --- a/tests/regression_tests/plot/results_true.dat +++ b/tests/regression_tests/plot/results_true.dat @@ -1 +1 @@ -ae51cbf147d69804b03f9d88c03d98533665a9822438d1d8b70b40928a1db35d20e18d12ee697b5e1e217799d6e8c352524e8d2227328e719e85c4e6d96f825c \ No newline at end of file +fd61f6a5de632f06a39e2a938480ba3dc314810655e684291e4b255bb8c142a16f10fec414c587ecd727fa559d3760d6e5b43f53fe86e90238a69d9c5698db8e \ No newline at end of file diff --git a/tests/regression_tests/plot/tallies.xml b/tests/regression_tests/plot/tallies.xml new file mode 100644 index 000000000..b7e678ca0 --- /dev/null +++ b/tests/regression_tests/plot/tallies.xml @@ -0,0 +1,20 @@ + + + + + -10 10 + -10 10 + -10 0 5 7.5 8.75 10 + + + + mesh + 2 + + + + 1 + total + + +