diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3287b0d2a..b6e7c73a9 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -37,62 +37,9 @@ extern std::unordered_map mesh_map; } // namespace model -class Mesh { -public: - // Constructor - Mesh() {}; // empty constructor - Mesh(pugi::xml_node node); - - //! Determine which bins were crossed by a particle - // - //! \param[in] p Particle to check - //! \param[out] bins Bins that were crossed - //! \param[out] lengths Fraction of tracklength in each bin - virtual void bins_crossed(const Particle* p, std::vector& bins, - std::vector& lengths) const = 0; - - //! Write mesh data to an HDF5 group - // - //! \param[in] group HDF5 group - virtual void to_hdf5(hid_t group) const = 0; - - //! Determine which surface bins were crossed by a particle - // - //! \param[in] p Particle to check - //! \param[out] bins Surface bins that were crossed - virtual void - surface_bins_crossed(const Particle* p, std::vector& bins) const = 0; - - //! Get bin at a given position in space - // - //! \param[in] r Position to get bin for - //! \return Mesh bin - virtual int get_bin(Position r) const = 0; - - virtual std::string get_label_for_bin(int bin) const = 0; - - //! 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 - virtual xt::xarray - count_sites(const std::vector& bank, bool* outside) const; - - virtual double get_volume_frac(int bin = -1) const = 0; - - virtual int num_bins() const = 0; - - int id_ {-1}; //!< User-specified ID - int n_dimension_; //!< Number of dimensions -}; - -//============================================================================== -//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes -//============================================================================== - -class RegularMesh : public Mesh { +class Mesh +{ public: // Constructors and destructor Mesh() = default; @@ -171,6 +118,60 @@ public: xt::xtensor upper_right_; //!< Upper-right coordinates of mesh }; +// class Mesh { + +// public: +// // Constructors +// Mesh() = default; +// Mesh(pugi::xml_node node); + +// // Destructor +// virtual ~Mesh() = default; + +// //! Determine which bins were crossed by a particle +// // +// //! \param[in] p Particle to check +// //! \param[out] bins Bins that were crossed +// //! \param[out] lengths Fraction of tracklength in each bin +// virtual void bins_crossed(const Particle* p, std::vector& bins, +// std::vector& lengths) const = 0; + +// //! Write mesh data to an HDF5 group +// // +// //! \param[in] group HDF5 group +// virtual void to_hdf5(hid_t group) const = 0; + +// //! Determine which surface bins were crossed by a particle +// // +// //! \param[in] p Particle to check +// //! \param[out] bins Surface bins that were crossed +// virtual void +// surface_bins_crossed(const Particle* p, std::vector& bins) const = 0; + +// //! Get bin at a given position in space +// // +// //! \param[in] r Position to get bin for +// //! \return Mesh bin +// virtual int get_bin(Position r) const = 0; + +// virtual std::string get_label_for_bin(int bin) const = 0; + +// //! 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 +// virtual xt::xarray +// count_sites(const std::vector& bank, bool* outside) const; + +// virtual double get_volume_frac(int bin = -1) const = 0; + +// virtual int num_bins() const = 0; + +// int id_ {-1}; //!< User-specified ID +// int n_dimension_; //!< Number of dimensions +// }; + //============================================================================== //! Tessellation of n-dimensional Euclidean space by congruent squares or cubes //============================================================================== @@ -393,6 +394,15 @@ intersect_track(const moab::CartVect& start, //! \return MOAB EntityHandle of tet moab::EntityHandle get_ent_handle_from_bin(int bin) const; + int get_bin_from_indices(const int* ijk) const override; + + void get_indices(Position r, int* ijk, bool* in_mesh) const override; + + void get_indices_from_bin(int bin, int* ijk) const override; + + std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const override; + //! Builds a KDTree for all tetrahedra in the mesh. All //! triangles representing 2D faces of the mesh are //! added to the tree as well. @@ -420,9 +430,11 @@ public: std::string get_label_for_bin(int bin) const; - double get_volume_frac(int bin = -1) const; + // double get_volume_frac(int bin = -1) const; - int num_bins() const; + int n_bins() const override; + + int n_surface_bins() const override; std::string filename_; // shape = {m}; // Create array of zeros @@ -850,28 +848,22 @@ RegularMesh::count_sites(const Particle::Bank* bank, int64_t length, return counts; } -std::string RegularMesh::get_label_for_bin(int bin) const { - int ijk[n_dimension_]; - get_indices_from_bin(bin, ijk); +// std::string RegularMesh::get_label_for_bin(int bin) const { +// int ijk[n_dimension_]; +// get_indices_from_bin(bin, ijk); - std::stringstream out; - out << "Mesh Index (" << ijk[0]; - if (n_dimension_ > 1) out << ", " << ijk[1]; - if (n_dimension_ > 2) out << ", " << ijk[2]; - out << ")"; +// std::stringstream out; +// out << "Mesh Index (" << ijk[0]; +// if (n_dimension_ > 1) out << ", " << ijk[1]; +// if (n_dimension_ > 2) out << ", " << ijk[2]; +// out << ")"; - return out.str(); -} +// return out.str(); +// } -double RegularMesh::get_volume_frac(int bin) const { - return volume_frac_; -} - -int RegularMesh::num_bins() const { - int n_bins = 1; - for (auto v : shape_) n_bins *= v; - return n_bins; -} +// double RegularMesh::get_volume_frac(int bin) const { +// return volume_frac_; +// } //============================================================================== // RectilinearMesh implementation @@ -1936,10 +1928,35 @@ UnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) return in_tet; } +int +UnstructuredMesh::get_bin_from_indices(const int* ijk) const { + if (ijk[0] >= n_bins()) { + std::stringstream s; + s << "Invalid bin: " << ijk[0]; + fatal_error(s); + } + int bin = ehs_[ijk[0]] - ehs_[0]; + return bin; +} + +void +UnstructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const { + int bin = get_bin(r); + ijk[0]= bin; + *in_mesh = bin != -1; +} + +void UnstructuredMesh::get_indices_from_bin(int bin, int* ijk) const { + ijk[0] = bin; +} + +std::pair, std::vector> +UnstructuredMesh::plot(Position plot_ll, Position plot_ur) const { return {}; } + int UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { int bin = eh - ehs_[0]; - if (bin >= num_bins()) { + if (bin >= n_bins()) { std::stringstream s; s << "Invalid bin: " << bin; fatal_error(s); @@ -1949,7 +1966,7 @@ UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { moab::EntityHandle UnstructuredMesh::get_ent_handle_from_bin(int bin) const { - if (bin >= num_bins()) { + if (bin >= n_bins()) { std::stringstream s; s << "Invalid bin index: " << bin; fatal_error(s); @@ -1957,19 +1974,31 @@ UnstructuredMesh::get_ent_handle_from_bin(int bin) const { return ehs_[bin]; } -double UnstructuredMesh::get_volume_frac(int bin) const { - if (bin == -1) { return 0.0; } - if (bin > ehs_.size() || bin < -1) { - std::stringstream msg; - msg << "Invalid bin " << bin << " for umesh with id " << id_; - fatal_error(msg); - } - moab::EntityHandle tet = get_ent_handle_from_bin(bin); - return tet_volume(tet); +// double UnstructuredMesh::get_volume_frac(int bin) const { +// if (bin == -1) { return 0.0; } +// if (bin > ehs_.size() || bin < -1) { +// std::stringstream msg; +// msg << "Invalid bin " << bin << " for umesh with id " << id_; +// fatal_error(msg); +// } +// moab::EntityHandle tet = get_ent_handle_from_bin(bin); +// return tet_volume(tet); +// } + +int UnstructuredMesh::n_bins() const { + return ehs_.size(); } -int UnstructuredMesh::num_bins() const { - return ehs_.size(); +int UnstructuredMesh::n_surface_bins() const { + // collect all triangles in the set of tets for this mesh + moab::Range tris; + moab::ErrorCode rval; + rval = mbi_->get_entities_by_type(0, moab::MBTRI, tris); + if (rval != moab::MB_SUCCESS) { + warning("Failed to get all triangles in the mesh instance"); + return -1; + } + return 2 * tris.size(); } #endif