diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index b83443cec6..293ae3a298 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -18,6 +18,7 @@ #ifdef DAGMC #include "moab/Core.hpp" #include "moab/AdaptiveKDTree.hpp" +#include "moab/Matrix3.hpp" #endif namespace openmc { @@ -26,11 +27,13 @@ namespace openmc { // Global variables //============================================================================== + class Mesh; namespace model { extern std::vector> meshes; + extern std::unordered_map mesh_map; } // namespace model @@ -50,14 +53,6 @@ public: virtual void bins_crossed(const Particle* p, std::vector& bins, std::vector& lengths) const = 0; - //! Check where a line segment intersects the mesh and if it intersects at all - // - //! \param[in,out] r0 In: starting position, out: intersection point - //! \param[in] r1 Ending position - //! \param[out] ijk Indices of the mesh bin containing the intersection point - //! \return Whether the line segment connecting r0 and r1 intersects mesh - virtual bool intersects(Position& r0, Position r1, int* ijk) const = 0; - //! Write mesh data to an HDF5 group // //! \param[in] group HDF5 group @@ -69,14 +64,6 @@ public: //! \return Mesh bin virtual int get_bin(Position r) 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 = 0; - int id_ {-1}; //!< User-specified ID int n_dimension_; //!< Number of dimensions }; @@ -281,6 +268,7 @@ private: #ifdef DAGMC class UnstructuredMesh : public Mesh { +public: UnstructuredMesh() { }; UnstructuredMesh(pugi::xml_node); @@ -294,24 +282,20 @@ class UnstructuredMesh : public Mesh { bool intersects(Position& r0, Position r1, int* ijk); + bool point_in_tet(const Position& r, moab::EntityHandle tet) const; + //! Write mesh data to an HDF5 group // //! \param[in] group HDF5 group - void to_hdf5(hid_t group); + void to_hdf5(hid_t group) 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); + int get_bin(Position r) 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); + void compute_barycentric_data(const moab::Range& all_tets); int get_bin_from_ent_handle(moab::EntityHandle eh) const; @@ -325,6 +309,7 @@ private: moab::EntityHandle meshset_; std::shared_ptr mbi_; std::unique_ptr kdtree_; + std::vector baryc_data_; }; #endif diff --git a/src/mesh.cpp b/src/mesh.cpp index 5ed8086764..1f480fbc3d 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -32,6 +32,7 @@ namespace openmc { namespace model { + std::vector> meshes; std::unordered_map mesh_map; @@ -1605,6 +1606,89 @@ UnstructuredMesh::bins_crossed(const Particle* p, std::vector& bins, }; + + +int +UnstructuredMesh::get_bin(Position r) const { + moab::CartVect pnt(r.x, r.y, r.z); + moab::AdaptiveKDTreeIter kdtree_iter; + moab::ErrorCode rval = kdtree_->point_search(pnt.array(), kdtree_iter); + if (rval != moab::MB_SUCCESS) { return -1; } + + moab::EntityHandle leaf = kdtree_iter.handle(); + moab::Range tets; + rval = mbi_->get_entities_by_dimension(leaf, 3, tets); + for (const auto& tet : tets) { + if (point_in_tet(r, tet)) { + return get_bin_from_ent_handle(tet); + } + } + + return -1; +} + +void +UnstructuredMesh::compute_barycentric_data(const moab::Range& all_tets) { + moab::ErrorCode rval; + for (auto& tet : all_tets) { + moab::Range verts; + rval = mbi_->get_connectivity(&tet, 1, verts); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get connectivity of tet on umesh: " + filename_); + } + + moab::CartVect p[4]; + rval = mbi_->get_coords(verts, p[0].array()); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get coordinates of a tet in umesh: " + filename_); + } + + moab::Matrix3 a(p[1] - p[0], p[2] - p[0], p[3] - p[0], true); + + baryc_data_.push_back(a.transpose().inverse()); + } +} + +// TODO: write this function +void +UnstructuredMesh::to_hdf5(hid_t group) const { } + +bool +UnstructuredMesh::point_in_tet(const Position& r, moab::EntityHandle tet) const { + + moab::ErrorCode rval; + + // get tet vertices + moab::Range verts; + rval = mbi_->get_connectivity(&tet, 1, verts); + if (rval != moab::MB_SUCCESS) { + warning("Failed to get vertices of tet in umesh: " + filename_); + return false; + } + + moab::EntityHandle v_zero = verts[0]; + moab::CartVect p_zero; + rval = mbi_->get_coords(&v_zero, 1, p_zero.array()); + if (rval != moab::MB_SUCCESS) { + warning("Failed to get coordinates of a vertex in umesh: " + filename_); + return false; + } + + moab::CartVect pos(r.x, r.y, r.z); + + // look up barycentric data + int idx = get_bin_from_ent_handle(tet); + const moab::Matrix3& a_inv = baryc_data_[idx]; + + moab::CartVect bary_coords = a_inv * (pos - p_zero); + + bool in_tet = (bary_coords[0] >= 0 && bary_coords[1] >= 0 && bary_coords[2] >= 0 && + bary_coords[0] + bary_coords[1] + bary_coords[2] <= 1.); + + return in_tet; + +} + int UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { auto pos = ehs_.find(eh); @@ -1638,6 +1722,9 @@ void read_meshes(pugi::xml_node root) model::meshes.push_back(std::make_unique(node)); } else if (mesh_type == "rectilinear") { model::meshes.push_back(std::make_unique(node)); + } + else if (mesh_type == "unstructured") { + model::meshes.push_back(std::make_unique(node)); } else { fatal_error("Invalid mesh type: " + mesh_type); }