crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Start implementing RectilinearMesh

Browse source
This commit is contained in:
Sterling Harper 2019-05-24 10:50:57 -04:00
parent 7b6ec80acf
commit ce2bf83254
7 changed files with 588 additions and 82 deletions
Download patch file Download diff file Expand all files Collapse all files

82
include/openmc/mesh.h
View file

@ -21,20 +21,48 @@ namespace openmc {
// Global variables
//==============================================================================
class RegularMesh;
//class RegularMesh;
class Mesh;
namespace model {
extern std::vector<std::unique_ptr<RegularMesh>> meshes;
//extern std::vector<std::unique_ptr<RegularMesh>> meshes;
extern std::vector<std::unique_ptr<Mesh>> meshes;
extern std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace model
class Mesh
{
public:
// Constructors
Mesh() = default;
Mesh(pugi::xml_node node);
virtual void bins_crossed(const Particle* p, std::vector<int>& bins,
std::vector<double>& lengths) const = 0;
virtual void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const = 0;
virtual int get_bin(Position r) const = 0;
virtual int get_bin_from_indices(const int* ijk) const = 0;
virtual void get_indices(Position r, int* ijk, bool* in_mesh) const = 0;
virtual void get_indices_from_bin(int bin, int* ijk) const = 0;
virtual void to_hdf5(hid_t group) const = 0;
int id_ {-1}; //!< User-specified ID
};
//==============================================================================
//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
//==============================================================================
class RegularMesh {
class RegularMesh : public Mesh
{
public:
// Constructors
RegularMesh() = default;
@ -47,39 +75,39 @@ public:
//! \param[in] p Particle to check
//! \param[out] bins Bins that were crossed
//! \param[out] lengths Fraction of tracklength in each bin
void bins_crossed(const Particle* p, std::vector<int>& bins,
virtual void bins_crossed(const Particle* p, std::vector<int>& bins,
std::vector<double>& lengths) const;
//! Determine which surface bins were crossed by a particle
//
//! \param[in] p Particle to check
//! \param[out] bins Surface bins that were crossed
void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
virtual void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
//! Get bin at a given position in space
//
//! \param[in] r Position to get bin for
//! \return Mesh bin
int get_bin(Position r) const;
virtual int get_bin(Position r) const;
//! Get bin given mesh indices
//
//! \param[in] Array of mesh indices
//! \return Mesh bin
int get_bin_from_indices(const int* ijk) const;
virtual int get_bin_from_indices(const int* ijk) const;
//! Get mesh indices given a position
//
//! \param[in] r Position to get indices for
//! \param[out] ijk Array of mesh indices
//! \param[out] in_mesh Whether position is in mesh
void get_indices(Position r, int* ijk, bool* in_mesh) const;
virtual void get_indices(Position r, int* ijk, bool* in_mesh) const;
//! Get mesh indices corresponding to a mesh bin
//
//! \param[in] bin Mesh bin
//! \param[out] ijk Mesh indices
void get_indices_from_bin(int bin, int* ijk) const;
virtual void get_indices_from_bin(int bin, int* ijk) const;
//! Check where a line segment intersects the mesh and if it intersects at all
//
@ -92,19 +120,18 @@ public:
//! Write mesh data to an HDF5 group
//
//! \param[in] group HDF5 group
void to_hdf5(hid_t group) const;
virtual void to_hdf5(hid_t group) const;
//! Count number of bank sites in each mesh bin / energy bin
//
//! \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,
bool* outside) const;
virtual xt::xarray<double>
count_sites(const std::vector<Particle::Bank>& bank, bool* outside) const;
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
double volume_frac_; //!< Volume fraction of each mesh element
int n_dimension_; //!< Number of dimensions
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
@ -116,6 +143,33 @@ private:
bool intersects_3d(Position& r0, Position r1, int* ijk) const;
};
class RectilinearMesh : public Mesh
{
public:
// Constructors
RectilinearMesh(pugi::xml_node node);
virtual void bins_crossed(const Particle* p, std::vector<int>& bins,
std::vector<double>& lengths) const;
virtual void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
virtual int get_bin(Position r) const;
virtual int get_bin_from_indices(const int* ijk) const;
virtual void get_indices(Position r, int* ijk, bool* in_mesh) const;
virtual void get_indices_from_bin(int bin, int* ijk) const;
virtual void to_hdf5(hid_t group) const;
bool intersects(Position& r0, Position r1, int* ijk) const;
xt::xarray<int> shape_; //!< Number of mesh elements in each dimension
std::vector<std::vector<double>> grid_;
};
//==============================================================================
// Non-member functions
//==============================================================================

23
src/eigenvalue.cpp
View file

@ -532,12 +532,13 @@ int openmc_get_keff(double* k_combined)
void shannon_entropy()
{
// Get pointer to entropy mesh
auto& m = model::meshes[settings::index_entropy_mesh];
// Get reference to entropy mesh
auto& m = *dynamic_cast<const RegularMesh*>(
model::meshes[settings::index_entropy_mesh].get());
// Get source weight in each mesh bin
bool sites_outside;
xt::xtensor<double, 1> p = m->count_sites(simulation::fission_bank,
xt::xtensor<double, 1> p = m.count_sites(simulation::fission_bank,
&sites_outside);
// display warning message if there were sites outside entropy box
@ -564,7 +565,8 @@ void shannon_entropy()
void ufs_count_sites()
{
auto &m = model::meshes[settings::index_ufs_mesh];
auto& m = *dynamic_cast<const RegularMesh*>(
model::meshes[settings::index_entropy_mesh].get());
if (simulation::current_batch == 1 && simulation::current_gen == 1) {
// On the first generation, just assume that the source is already evenly
@ -572,12 +574,12 @@ void ufs_count_sites()
// biased
auto s = xt::view(simulation::source_frac, xt::all());
s = m->volume_frac_;
s = m.volume_frac_;
} else {
// count number of source sites in each ufs mesh cell
bool sites_outside;
simulation::source_frac = m->count_sites(simulation::source_bank,
simulation::source_frac = m.count_sites(simulation::source_bank,
&sites_outside);
// Check for sites outside of the mesh
@ -587,7 +589,7 @@ void ufs_count_sites()
#ifdef OPENMC_MPI
// Send source fraction to all processors
int n_bins = xt::prod(m->shape_)();
int n_bins = xt::prod(m.shape_)();
MPI_Bcast(simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
#endif
@ -605,17 +607,18 @@ void ufs_count_sites()
double ufs_get_weight(const Particle* p)
{
auto& m = model::meshes[settings::index_ufs_mesh];
auto& m = *dynamic_cast<const RegularMesh*>(
model::meshes[settings::index_entropy_mesh].get());
// Determine indices on ufs mesh for current location
int mesh_bin = m->get_bin(p->r());
int mesh_bin = m.get_bin(p->r());
if (mesh_bin < 0) {
p->write_restart();
fatal_error("Source site outside UFS mesh!");
}
if (simulation::source_frac(mesh_bin) != 0.0) {
return m->volume_frac_ / simulation::source_frac(mesh_bin);
return m.volume_frac_ / simulation::source_frac(mesh_bin);
} else {
return 1.0;
}

500
src/mesh.cpp
View file

@ -31,7 +31,7 @@ namespace openmc {
namespace model {
std::vector<std::unique_ptr<RegularMesh>> meshes;
std::vector<std::unique_ptr<Mesh>> meshes;
std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace model
@ -63,10 +63,10 @@ inline bool check_intersection_point(double x1, double x0, double y1,
}
//==============================================================================
// RegularMesh implementation
// Mesh implementation
//==============================================================================
RegularMesh::RegularMesh(pugi::xml_node node)
Mesh::Mesh(pugi::xml_node node)
{
// Copy mesh id
if (check_for_node(node, "id")) {
@ -78,17 +78,15 @@ RegularMesh::RegularMesh(pugi::xml_node node)
std::to_string(id_));
}
}
}
// Read mesh type
if (check_for_node(node, "type")) {
auto temp = get_node_value(node, "type", true, true);
if (temp == "regular") {
// TODO: move elsewhere
} else {
fatal_error("Invalid mesh type: " + temp);
}
}
//==============================================================================
// RegularMesh implementation
//==============================================================================
RegularMesh::RegularMesh(pugi::xml_node node)
: Mesh {node}
{
// Determine number of dimensions for mesh
if (check_for_node(node, "dimension")) {
shape_ = get_node_xarray<int>(node, "dimension");
@ -783,6 +781,425 @@ RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
return counts;
}
//==============================================================================
// RectilinearMesh implementation
//==============================================================================
RectilinearMesh::RectilinearMesh(pugi::xml_node node)
: Mesh {node}
{
grid_.resize(3);
grid_[0] = get_node_array<double>(node, "x_grid");
grid_[1] = get_node_array<double>(node, "y_grid");
grid_[2] = get_node_array<double>(node, "z_grid");
shape_ = {static_cast<int>(grid_[0].size()) - 1,
static_cast<int>(grid_[1].size()) - 1,
static_cast<int>(grid_[2].size()) - 1};
}
void RectilinearMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
std::vector<double>& lengths) const
{
// ========================================================================
// Determine where the track intersects the mesh and if it intersects at all.
// Copy the starting and ending coordinates of the particle.
Position last_r {p->r_last_};
Position r {p->r()};
Direction u {p->u()};
// Compute the length of the entire track.
double total_distance = (r - last_r).norm();
// While determining if this track intersects the mesh, offset the starting
// and ending coords by a bit. This avoid finite-precision errors that can
// occur when the mesh surfaces coincide with lattice or geometric surfaces.
Position r0 = last_r + TINY_BIT*u;
Position r1 = r - TINY_BIT*u;
// Determine the mesh indices for the starting and ending coords.
int ijk0[3], ijk1[3];
bool start_in_mesh;
get_indices(r0, ijk0, &start_in_mesh);
bool end_in_mesh;
get_indices(r1, ijk1, &end_in_mesh);
// Reset coordinates and check for a mesh intersection if necessary.
if (start_in_mesh) {
// The initial coords lie in the mesh, use those coords for tallying.
r0 = last_r;
} else {
// The initial coords do not lie in the mesh. Check to see if the particle
// eventually intersects the mesh and compute the relevant coords and
// indices.
if (!intersects(r0, r1, ijk0)) return;
}
r1 = r;
// ========================================================================
// Find which mesh cells are traversed and the length of each traversal.
while (true) {
if (std::equal(ijk0, ijk0+3, ijk1)) {
// The track ends in this cell. Use the particle end location rather
// than the mesh surface and stop iterating.
double distance = (r1 - r0).norm();
bins.push_back(get_bin_from_indices(ijk0));
lengths.push_back(distance / total_distance);
break;
}
// The track exits this cell. Determine the distance to each mesh surface.
double d[3];
for (int k = 0; k < 3; ++k) {
if (std::fabs(u[k]) < FP_PRECISION) {
d[k] = INFTY;
} else if (u[k] > 0) {
double xyz_cross = grid_[k][ijk0[k]];
d[k] = (xyz_cross - r0[k]) / u[k];
} else {
double xyz_cross = grid_[k][ijk0[k] - 1];
d[k] = (xyz_cross - r0[k]) / u[k];
}
}
// Pick the closest mesh surface and append this traversal to the output.
auto j = std::min_element(d, d+3) - d;
double distance = d[j];
bins.push_back(get_bin_from_indices(ijk0));
lengths.push_back(distance / total_distance);
// Translate to the oncoming mesh surface.
r0 += distance * u;
// Increment the indices into the next mesh cell.
if (u[j] > 0.0) {
++ijk0[j];
} else {
--ijk0[j];
}
// If the next indices are invalid, then the track has left the mesh and
// we are done.
bool in_mesh = true;
for (int i = 0; i < 3; ++i) {
if (ijk0[i] < 1 || ijk0[i] > shape_[i]) {
in_mesh = false;
break;
}
}
if (!in_mesh) break;
}
}
void RectilinearMesh::surface_bins_crossed(const Particle* p,
std::vector<int>& bins) const
{
// ========================================================================
// Determine if the track intersects the tally mesh.
// Copy the starting and ending coordinates of the particle.
Position r0 {p->r_last_current_};
Position r1 {p->r()};
Direction u {p->u()};
// Determine indices for starting and ending location.
int ijk0[3], ijk1[3];
bool start_in_mesh;
get_indices(r0, ijk0, &start_in_mesh);
bool end_in_mesh;
get_indices(r1, ijk1, &end_in_mesh);
// Check if the track intersects any part of the mesh.
if (!start_in_mesh) {
Position r0_copy = r0;
int ijk0_copy[3];
for (int i = 0; i < 3; ++i) ijk0_copy[i] = ijk0[i];
if (!intersects(r0_copy, r1, ijk0_copy)) return;
}
// ========================================================================
// Find which mesh surfaces are crossed.
// Calculate number of surface crossings
int n_cross = 0;
for (int i = 0; i < 3; ++i) n_cross += std::abs(ijk1[i] - ijk0[i]);
if (n_cross == 0) return;
// Bounding coordinates
Position xyz_cross;
for (int i = 0; i < 3; ++i) {
if (u[i] > 0.0) {
xyz_cross[i] = grid_[i][ijk0[i]];
} else {
xyz_cross[i] = grid_[i][ijk0[i] - 1];
}
}
for (int j = 0; j < n_cross; ++j) {
// Set the distances to infinity
Position d {INFTY, INFTY, INFTY};
// Determine closest bounding surface. We need to treat
// special case where the cosine of the angle is zero since this would
// result in a divide-by-zero.
double distance = INFTY;
for (int i = 0; i < 3; ++i) {
if (u[i] == 0) {
d[i] = INFTY;
} else {
d[i] = (xyz_cross[i] - r0[i])/u[i];
}
distance = std::min(distance, d[i]);
}
// Loop over the dimensions
for (int i = 0; i < 3; ++i) {
// Check whether distance is the shortest distance
if (distance == d[i]) {
// Check whether the current indices are within the mesh bounds
bool in_mesh = true;
for (int j = 0; j < 3; ++j) {
if (ijk0[j] < 1 || ijk0[j] > shape_[j]) {
in_mesh = false;
break;
}
}
// Check whether particle is moving in positive i direction
if (u[i] > 0) {
// Outward current on i max surface
if (in_mesh) {
int i_surf = 4*i + 3;
int i_mesh = get_bin_from_indices(ijk0);
int i_bin = 4*3*i_mesh + i_surf - 1;
bins.push_back(i_bin);
}
// Advance position
++ijk0[i];
xyz_cross[i] = grid_[i][ijk0[i]];
in_mesh = true;
for (int j = 0; j < 3; ++j) {
if (ijk0[j] < 1 || ijk0[j] > shape_[j]) {
in_mesh = false;
break;
}
}
// If the particle crossed the surface, tally the inward current on
// i min surface
if (in_mesh) {
int i_surf = 4*i + 2;
int i_mesh = get_bin_from_indices(ijk0);
int i_bin = 4*3*i_mesh + i_surf - 1;
bins.push_back(i_bin);
}
} else {
// The particle is moving in the negative i direction
// Outward current on i min surface
if (in_mesh) {
int i_surf = 4*i + 1;
int i_mesh = get_bin_from_indices(ijk0);
int i_bin = 4*3*i_mesh + i_surf - 1;
bins.push_back(i_bin);
}
// Advance position
--ijk0[i];
xyz_cross[i] = grid_[i][ijk0[i] - 1];
in_mesh = true;
for (int j = 0; j < 3; ++j) {
if (ijk0[j] < 1 || ijk0[j] > shape_[j]) {
in_mesh = false;
break;
}
}
// If the particle crossed the surface, tally the inward current on
// i max surface
if (in_mesh) {
int i_surf = 4*i + 4;
int i_mesh = get_bin_from_indices(ijk0);
int i_bin = 4*3*i_mesh + i_surf - 1;
bins.push_back(i_bin);
}
}
}
}
// Calculate new coordinates
r0 += distance * u;
}
}
int RectilinearMesh::get_bin(Position r) const
{
// Determine indices
int ijk[3];
bool in_mesh;
get_indices(r, ijk, &in_mesh);
if (!in_mesh) return -1;
// Convert indices to bin
return get_bin_from_indices(ijk);
}
int RectilinearMesh::get_bin_from_indices(const int* ijk) const
{
return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1;
}
void RectilinearMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
{
*in_mesh = true;
for (int i = 0; i < 3; ++i) {
if (r[i] < grid_[i].front() || r[i] > grid_[i].back()) {
ijk[i] = -1;
*in_mesh = false;
} else {
ijk[i] = lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1;
}
}
}
void RectilinearMesh::get_indices_from_bin(int bin, int* ijk) const
{
ijk[0] = bin % shape_[0] + 1;
ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1;
ijk[2] = bin / (shape_[0] * shape_[1]) + 1;
}
void RectilinearMesh::to_hdf5(hid_t group) const
{
hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
write_dataset(mesh_group, "type", "rectilinear");
close_group(mesh_group);
}
bool RectilinearMesh::intersects(Position& r0, Position r1, int* ijk) const
{
// Copy coordinates of starting point
double x0 = r0.x;
double y0 = r0.y;
double z0 = r0.z;
// Copy coordinates of ending point
double x1 = r1.x;
double y1 = r1.y;
double z1 = r1.z;
// Copy coordinates of mesh lower_left
double xm0 = grid_[0].front();
double ym0 = grid_[1].front();
double zm0 = grid_[2].front();
// Copy coordinates of mesh upper_right
double xm1 = grid_[0].back();
double ym1 = grid_[1].back();
double zm1 = grid_[2].back();
double min_dist = INFTY;
// Check if line intersects left surface -- calculate the intersection point
// (y,z)
if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
ijk[0] = 1;
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
}
}
}
// Check if line intersects back surface -- calculate the intersection point
// (x,z)
if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
ijk[1] = 1;
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
}
}
}
// Check if line intersects bottom surface -- calculate the intersection
// point (x,y)
if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
ijk[2] = 1;
}
}
}
// Check if line intersects right surface -- calculate the intersection point
// (y,z)
if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
ijk[0] = shape_[0];
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
}
}
}
// Check if line intersects front surface -- calculate the intersection point
// (x,z)
if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
ijk[1] = shape_[1];
ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
}
}
}
// Check if line intersects top surface -- calculate the intersection point
// (x,y)
if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
ijk[2] = shape_[2];
}
}
}
return min_dist < INFTY;
}
//==============================================================================
// C API functions
//==============================================================================
@ -846,8 +1263,8 @@ openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
set_errmsg("Index in meshes array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*dims = model::meshes[index]->shape_.data();
*n = model::meshes[index]->n_dimension_;
*dims = dynamic_cast<RegularMesh*>(model::meshes[index].get())->shape_.data();
*n = dynamic_cast<RegularMesh*>(model::meshes[index].get())->n_dimension_;
return 0;
}
@ -862,9 +1279,9 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
// Copy dimension
std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
auto& m = model::meshes[index];
m->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
m->n_dimension_ = m->shape_.size();
auto& m = *dynamic_cast<RegularMesh*>(model::meshes[index].get());
m.shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
m.n_dimension_ = m.shape_.size();
return 0;
}
@ -878,16 +1295,17 @@ openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width,
return OPENMC_E_OUT_OF_BOUNDS;
}
auto& m = model::meshes[index];
if (m->lower_left_.dimension() == 0) {
//auto& m = model::meshes[index];
auto& m = *dynamic_cast<RegularMesh*>(model::meshes[index].get());
if (m.lower_left_.dimension() == 0) {
set_errmsg("Mesh parameters have not been set.");
return OPENMC_E_ALLOCATE;
}
*ll = m->lower_left_.data();
*ur = m->upper_right_.data();
*width = m->width_.data();
*n = m->n_dimension_;
*ll = m.lower_left_.data();
*ur = m.upper_right_.data();
*width = m.width_.data();
*n = m.n_dimension_;
return 0;
}
@ -901,20 +1319,21 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
return OPENMC_E_OUT_OF_BOUNDS;
}
auto& m = model::meshes[index];
//auto& m = model::meshes[index];
auto& m = *dynamic_cast<RegularMesh*>(model::meshes[index].get());
std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
if (ll && ur) {
m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_;
m.lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
m.upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_;
} else if (ll && width) {
m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
m->upper_right_ = m->lower_left_ + m->shape_ * m->width_;
m.lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
m.width_ = xt::adapt(width, n, xt::no_ownership(), shape);
m.upper_right_ = m.lower_left_ + m.shape_ * m.width_;
} else if (ur && width) {
m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
m->lower_left_ = m->upper_right_ - m->shape_ * m->width_;
m.upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
m.width_ = xt::adapt(width, n, xt::no_ownership(), shape);
m.lower_left_ = m.upper_right_ - m.shape_ * m.width_;
} else {
set_errmsg("At least two parameters must be specified.");
return OPENMC_E_INVALID_ARGUMENT;
@ -930,8 +1349,21 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
void read_meshes(pugi::xml_node root)
{
for (auto node : root.children("mesh")) {
std::string mesh_type;
if (check_for_node(node, "type")) {
mesh_type = get_node_value(node, "type", true, true);
} else {
mesh_type = "regular";
}
// Read mesh and add to vector
model::meshes.push_back(std::make_unique<RegularMesh>(node));
if (mesh_type == "regular") {
model::meshes.push_back(std::make_unique<RegularMesh>(node));
} else if (mesh_type == "rectilinear") {
model::meshes.push_back(std::make_unique<RectilinearMesh>(node));
} else {
fatal_error("Invalid mesh type: " + mesh_type);
}
// Map ID to position in vector
model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;

16
src/plot.cpp
View file

@ -707,12 +707,12 @@ void draw_mesh_lines(Plot pl, ImageData& data)
Position width = ur_plot - ll_plot;
auto& m = model::meshes[pl.index_meshlines_mesh_];
auto& m = *dynamic_cast<RegularMesh*>(model::meshes[pl.index_meshlines_mesh_].get());
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);
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]++;
@ -722,13 +722,13 @@ void draw_mesh_lines(Plot pl, ImageData& data)
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] ) {
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;
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];

2
src/settings.cpp
View file

@ -508,7 +508,7 @@ void read_settings_xml()
}
if (index_entropy_mesh >= 0) {
auto& m = *model::meshes[index_entropy_mesh];
auto& m = *dynamic_cast<RegularMesh*>(model::meshes[index_entropy_mesh].get());
if (m.shape_.dimension() == 0) {
// If the user did not specify how many mesh cells are to be used in
// each direction, we automatically determine an appropriate number of

40
src/tallies/filter_mesh.cpp
View file

@ -55,27 +55,43 @@ MeshFilter::to_statepoint(hid_t filter_group) const
std::string
MeshFilter::text_label(int bin) const
{
auto& mesh = *model::meshes[mesh_];
int n_dim = mesh.n_dimension_;
if (auto* mesh = dynamic_cast<RegularMesh*>(model::meshes[mesh_].get())) {
int n_dim = mesh->n_dimension_;
int ijk[n_dim];
mesh.get_indices_from_bin(bin, ijk);
int ijk[n_dim];
mesh->get_indices_from_bin(bin, ijk);
std::stringstream out;
out << "Mesh Index (" << ijk[0];
if (n_dim > 1) out << ", " << ijk[1];
if (n_dim > 2) out << ", " << ijk[2];
out << ")";
std::stringstream out;
out << "Mesh Index (" << ijk[0];
if (n_dim > 1) out << ", " << ijk[1];
if (n_dim > 2) out << ", " << ijk[2];
out << ")";
return out.str();
return out.str();
} else if (auto* mesh = dynamic_cast<RectilinearMesh*>(model::meshes[mesh_].get())) {
int ijk[3];
mesh->get_indices_from_bin(bin, ijk);
std::stringstream out;
out << "Mesh Index (" << ijk[0] << ", " << ijk[1] << ", " << ijk[2] << ")";
return out.str();
}
}
void
MeshFilter::set_mesh(int32_t mesh)
{
mesh_ = mesh;
n_bins_ = 1;
for (auto dim : model::meshes[mesh_]->shape_) n_bins_ *= dim;
auto* m_ptr = model::meshes[mesh_].get();
if (auto* m = dynamic_cast<RegularMesh*>(m_ptr)) {
n_bins_ = 1;
for (auto dim : m->shape_) n_bins_ *= dim;
} else if (auto* m = dynamic_cast<RectilinearMesh*>(m_ptr)) {
n_bins_ = 1;
for (auto dim : m->shape_) n_bins_ *= dim;
}
}
//==============================================================================

7
src/tallies/filter_meshsurface.cpp
View file

@ -18,7 +18,7 @@ MeshSurfaceFilter::get_all_bins(const Particle* p, int estimator,
std::string
MeshSurfaceFilter::text_label(int bin) const
{
auto& mesh = *model::meshes[mesh_];
auto& mesh = *dynamic_cast<RegularMesh*>(model::meshes[mesh_].get());
int n_dim = mesh.n_dimension_;
// Get flattend mesh index and surface index.
@ -75,8 +75,9 @@ void
MeshSurfaceFilter::set_mesh(int32_t mesh)
{
mesh_ = mesh;
n_bins_ = 4 * model::meshes[mesh_]->n_dimension_;
for (auto dim : model::meshes[mesh_]->shape_) n_bins_ *= dim;
auto &m = *dynamic_cast<RegularMesh*>(model::meshes[mesh_].get());
n_bins_ = 4 * m.n_dimension_;
for (auto dim : m.shape_) n_bins_ *= dim;
}
//==============================================================================