mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Address comments on #1246
This commit is contained in:
Sterling Harper
parent
7cf1822410
commit
b804ddc952
5 changed files with 83 additions and 76 deletions
|
|
@ -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<double, 1> lower_left_; //!< Lower-left coordinates of mesh
|
||||
xt::xtensor<double, 1> 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<double> count_sites(const std::vector<Particle::Bank>& bank,
|
||||
xt::xtensor<double, 1> count_sites(const std::vector<Particle::Bank>& bank,
|
||||
bool* outside) const;
|
||||
|
||||
// Data members
|
||||
|
||||
double volume_frac_; //!< Volume fraction of each mesh element
|
||||
xt::xarray<int> shape_; //!< Number of mesh elements in each dimension
|
||||
xt::xarray<double> lower_left_; //!< Lower-left coordinates of mesh
|
||||
xt::xarray<double> upper_right_; //!< Upper-right coordinates of mesh
|
||||
xt::xarray<double> width_; //!< Width of each mesh element
|
||||
xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
|
||||
xt::xtensor<double, 1> 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<int> shape_; //!< Number of mesh elements in each dimension
|
||||
xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
|
||||
|
||||
private:
|
||||
std::vector<std::vector<double>> grid_;
|
||||
|
|
|
|||
|
|
@ -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), ...]
|
||||
|
||||
"""
|
||||
|
||||
|
|
|
|||
16
src/mesh.cpp
16
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<double>& 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<double>
|
||||
xt::xtensor<double, 1>
|
||||
RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
|
||||
bool* outside) const
|
||||
{
|
||||
|
|
@ -846,6 +846,18 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node)
|
|||
shape_ = {static_cast<int>(grid_[0].size()) - 1,
|
||||
static_cast<int>(grid_[1].size()) - 1,
|
||||
static_cast<int>(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<int>& bins,
|
||||
|
|
|
|||
102
src/plot.cpp
102
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<double>, std::vector<double>> 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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 <entropy> element "
|
||||
|
|
@ -498,15 +498,16 @@ void read_settings_xml()
|
|||
model::meshes.push_back(std::make_unique<RegularMesh>(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<RegularMesh*>(model::meshes[i_entropy_mesh].get());
|
||||
if (index_entropy_mesh >= 0) {
|
||||
auto* m = dynamic_cast<RegularMesh*>(
|
||||
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) {
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue