From ce2bf83254554f61ee66cde71af64e659c80853e Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 24 May 2019 10:50:57 -0400 Subject: [PATCH] Start implementing RectilinearMesh --- include/openmc/mesh.h | 82 ++++- src/eigenvalue.cpp | 23 +- src/mesh.cpp | 500 +++++++++++++++++++++++++++-- src/plot.cpp | 16 +- src/settings.cpp | 2 +- src/tallies/filter_mesh.cpp | 40 ++- src/tallies/filter_meshsurface.cpp | 7 +- 7 files changed, 588 insertions(+), 82 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index f53bba224..f7e41cf70 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -21,20 +21,48 @@ namespace openmc { // Global variables //============================================================================== -class RegularMesh; +//class RegularMesh; +class Mesh; namespace model { -extern std::vector> meshes; +//extern std::vector> meshes; +extern std::vector> meshes; extern std::unordered_map mesh_map; } // namespace model +class Mesh +{ +public: + // Constructors + Mesh() = default; + Mesh(pugi::xml_node node); + + virtual void bins_crossed(const Particle* p, std::vector& bins, + std::vector& lengths) const = 0; + + virtual void surface_bins_crossed(const Particle* p, std::vector& 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& bins, + virtual void bins_crossed(const Particle* p, std::vector& bins, std::vector& 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& bins) const; + virtual void surface_bins_crossed(const Particle* p, std::vector& 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 count_sites(const std::vector& bank, - bool* outside) const; + virtual xt::xarray + count_sites(const std::vector& 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 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 @@ -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& bins, + std::vector& lengths) const; + + virtual void surface_bins_crossed(const Particle* p, std::vector& 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 shape_; //!< Number of mesh elements in each dimension + std::vector> grid_; +}; + //============================================================================== // Non-member functions //============================================================================== diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 08a4d3f91..0f4b49f46 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -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( + model::meshes[settings::index_entropy_mesh].get()); // Get source weight in each mesh bin bool sites_outside; - xt::xtensor p = m->count_sites(simulation::fission_bank, + xt::xtensor 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( + 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( + 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; } diff --git a/src/mesh.cpp b/src/mesh.cpp index ddfa272a6..676c5888f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -31,7 +31,7 @@ namespace openmc { namespace model { -std::vector> meshes; +std::vector> meshes; std::unordered_map 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(node, "dimension"); @@ -783,6 +781,425 @@ RegularMesh::count_sites(const std::vector& bank, return counts; } +//============================================================================== +// RectilinearMesh implementation +//============================================================================== + +RectilinearMesh::RectilinearMesh(pugi::xml_node node) + : Mesh {node} +{ + grid_.resize(3); + grid_[0] = get_node_array(node, "x_grid"); + grid_[1] = get_node_array(node, "y_grid"); + grid_[2] = get_node_array(node, "z_grid"); + + shape_ = {static_cast(grid_[0].size()) - 1, + static_cast(grid_[1].size()) - 1, + static_cast(grid_[2].size()) - 1}; +} + +void RectilinearMesh::bins_crossed(const Particle* p, std::vector& bins, + std::vector& 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& 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(model::meshes[index].get())->shape_.data(); + *n = dynamic_cast(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 shape = {static_cast(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(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(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(model::meshes[index].get()); std::vector shape = {static_cast(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(node)); + if (mesh_type == "regular") { + model::meshes.push_back(std::make_unique(node)); + } else if (mesh_type == "rectilinear") { + model::meshes.push_back(std::make_unique(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; diff --git a/src/plot.cpp b/src/plot.cpp index 57466eca1..070726bcb 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -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(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]; diff --git a/src/settings.cpp b/src/settings.cpp index 4c041a080..27f8c31d2 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -508,7 +508,7 @@ void read_settings_xml() } if (index_entropy_mesh >= 0) { - auto& m = *model::meshes[index_entropy_mesh]; + auto& m = *dynamic_cast(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 diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index f69a96652..9658e552e 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -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(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(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(m_ptr)) { + n_bins_ = 1; + for (auto dim : m->shape_) n_bins_ *= dim; + } else if (auto* m = dynamic_cast(m_ptr)) { + n_bins_ = 1; + for (auto dim : m->shape_) n_bins_ *= dim; + } } //============================================================================== diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index 57f10b09c..8f4ae4e35 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -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(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(model::meshes[mesh_].get()); + n_bins_ = 4 * m.n_dimension_; + for (auto dim : m.shape_) n_bins_ *= dim; } //==============================================================================