From b1af80697caca15e47cdb08463fd3ebde1526cce Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 8 Apr 2020 14:07:37 -0500 Subject: [PATCH] Restructuring unstructured mesh inheritance a bit. --- include/openmc/mesh.h | 44 +++++++++------ src/mesh.cpp | 128 ++++++++++++++++++++---------------------- 2 files changed, 86 insertions(+), 86 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index e9d10b6c8..150555680 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -284,6 +284,11 @@ public: std::string bin_label(int bin) const override; + //! Write mesh data to an HDF5 group. + // + //! \param[in] group HDF5 group + void to_hdf5(hid_t group) const; + //! Add a variable to the mesh instance virtual void add_score(const std::string& var_name) = 0; @@ -296,20 +301,24 @@ public: virtual Position centroid(int bin) const = 0; + virtual double volume(int bin) const = 0; + + virtual std::string library() const = 0; + // Data members - std::string filename_; + std::string filename_; //!< Path to unstructured mesh file }; #ifdef DAGMC -class UnstructuredMesh : public UnstructuredMeshBase { +class MOABUnstructuredMesh : public UnstructuredMeshBase { public: - UnstructuredMesh() = default; - UnstructuredMesh(pugi::xml_node); - UnstructuredMesh(const std::string& filename); - ~UnstructuredMesh() = default; + MOABUnstructuredMesh() = default; + MOABUnstructuredMesh(pugi::xml_node); + MOABUnstructuredMesh(const std::string& filename); + ~MOABUnstructuredMesh() = default; void bins_crossed(const Particle& p, std::vector& bins, @@ -324,27 +333,26 @@ public: //! \param[out] bins Surface bins that were crossed void surface_bins_crossed(const Particle& p, std::vector& bins) const; - //! Write mesh data to an HDF5 group. - // - //! \param[in] group HDF5 group - void to_hdf5(hid_t group) const; - //! Get bin at a given position. // //! \param[in] r Position to get bin for //! \return Mesh bin int get_bin(Position r) const; + std::string library() const override; + int n_bins() const override; int n_surface_bins() const override; //! Retrieve a centroid for the mesh cell // - // \param[in] tet MOAB EntityHandle of the tetrahedron + // \param[in] bin Bin to return the volume for // \return The centroid of the element Position centroid(int bin) const override; + double volume(int bin) const override; + //! Add a score to the mesh instance void add_score(const std::string& score); @@ -364,8 +372,6 @@ public: //! Write the mesh with any current tally data void write(std::string base_filename) const; - std::string filename_; //!< Path to unstructured mesh file - bool intersects(Position& r0, Position r1, int* ijk); private: @@ -492,6 +498,8 @@ public: int get_bin(Position r) const override; + std::string library() const override; + int n_bins() const override; int n_surface_bins() const override; @@ -504,8 +512,6 @@ public: void write(std::string filename) const override; - void to_hdf5(hid_t group) const override; - //! Add a variable to the libmesh mesh instance void add_score(const std::string& var_name) override; @@ -514,6 +520,10 @@ public: std::vector values, std::vector std_dev) override; + Position centroid(int bin) const override; + + double volume(int bin) const override; + private: //! Setup method for the mesh. Builds data structures, @@ -551,8 +561,6 @@ private: double first(const libMesh::Node& a, const libMesh::Node& b) const; - Position centroid(int bin) const override; - // Data members private: diff --git a/src/mesh.cpp b/src/mesh.cpp index 955c264ab..0f92032e3 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -136,6 +136,30 @@ UnstructuredMeshBase::bin_label(int bin) const { return fmt::format("Mesh Index ({})", bin); }; +void +UnstructuredMeshBase::to_hdf5(hid_t group) const +{ + hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_)); + + write_dataset(mesh_group, "type", "unstructured"); + write_dataset(mesh_group, "filename", filename_); + write_dataset(mesh_group, "library", this->library()); + // write volume of each tet + std::vector tet_vols; + xt::xtensor centroids({this->n_bins(), 3}); + for (int i = 0; i < this->n_bins(); i++) { + tet_vols.emplace_back(this->volume(i)); + auto c = this->centroid(i); + xt::view(centroids, i, xt::all()) = xt::xarray({c.x, c.y, c.z}); + } + + write_dataset(mesh_group, "volumes", tet_vols); + write_dataset(mesh_group, "centroids", centroids); + + close_group(mesh_group); +} + + //============================================================================== // RegularMesh implementation //============================================================================== @@ -1404,7 +1428,7 @@ extern "C" int openmc_add_unstructured_mesh(const char filename[], bool valid_lib = false; #ifdef DAGMC if (library == "moab") { - model::meshes.push_back(std::move(std::make_unique(filename))); + model::meshes.push_back(std::move(std::make_unique(filename))); valid_lib = true; } #endif @@ -1590,16 +1614,16 @@ openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x, #ifdef DAGMC -UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : UnstructuredMeshBase(node) { +MOABUnstructuredMesh::MOABUnstructuredMesh(pugi::xml_node node) : UnstructuredMeshBase(node) { initialize(); } -UnstructuredMesh::UnstructuredMesh(const std::string& filename) { +MOABUnstructuredMesh::MOABUnstructuredMesh(const std::string& filename) { filename_ = filename; initialize(); } -void UnstructuredMesh::initialize() { +void MOABUnstructuredMesh::initialize() { // create MOAB instance mbi_ = std::make_unique(); // load unstructured mesh file @@ -1637,7 +1661,7 @@ void UnstructuredMesh::initialize() { } void -UnstructuredMesh::build_kdtree(const moab::Range& all_tets) +MOABUnstructuredMesh::build_kdtree(const moab::Range& all_tets) { moab::Range all_tris; int adj_dim = 2; @@ -1672,7 +1696,7 @@ UnstructuredMesh::build_kdtree(const moab::Range& all_tets) } void -UnstructuredMesh::intersect_track(const moab::CartVect& start, +MOABUnstructuredMesh::intersect_track(const moab::CartVect& start, const moab::CartVect& dir, double track_len, std::vector& hits) const { @@ -1703,7 +1727,7 @@ UnstructuredMesh::intersect_track(const moab::CartVect& start, } void -UnstructuredMesh::bins_crossed(const Particle& p, +MOABUnstructuredMesh::bins_crossed(const Particle* p, std::vector& bins, std::vector& lengths) const { @@ -1780,7 +1804,7 @@ UnstructuredMesh::bins_crossed(const Particle& p, }; moab::EntityHandle -UnstructuredMesh::get_tet(const Position& r) const +MOABUnstructuredMesh::get_tet(const Position& r) const { moab::CartVect pos(r.x, r.y, r.z); // find the leaf of the kd-tree for this position @@ -1807,7 +1831,13 @@ UnstructuredMesh::get_tet(const Position& r) const return 0; } -double UnstructuredMesh::tet_volume(moab::EntityHandle tet) const { +double MOABUnstructuredMesh::volume(int bin) const { + return tet_volume(get_ent_handle_from_bin(bin)); +} + +std::string MOABUnstructuredMesh::library() const { return "moab"; } + +double MOABUnstructuredMesh::tet_volume(moab::EntityHandle tet) const { std::vector conn; moab::ErrorCode rval = mbi_->get_connectivity(&tet, 1, conn); if (rval != moab::MB_SUCCESS) { @@ -1823,13 +1853,13 @@ double UnstructuredMesh::tet_volume(moab::EntityHandle tet) const { return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0])); } -void UnstructuredMesh::surface_bins_crossed(const Particle& p, std::vector& bins) const { +void MOABUnstructuredMesh::surface_bins_crossed(const Particle* p, std::vector& bins) const { // TODO: Implement triangle crossings here throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; } int -UnstructuredMesh::get_bin(Position r) const { +MOABUnstructuredMesh::get_bin(Position r) const { moab::EntityHandle tet = get_tet(r); if (tet == 0) { return -1; @@ -1839,7 +1869,7 @@ UnstructuredMesh::get_bin(Position r) const { } void -UnstructuredMesh::compute_barycentric_data(const moab::Range& tets) { +MOABUnstructuredMesh::compute_barycentric_data(const moab::Range& tets) { moab::ErrorCode rval; baryc_data_.clear(); @@ -1868,32 +1898,8 @@ UnstructuredMesh::compute_barycentric_data(const moab::Range& tets) { } } -void -UnstructuredMesh::to_hdf5(hid_t group) const -{ - hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_)); - - write_dataset(mesh_group, "type", "unstructured"); - write_dataset(mesh_group, "filename", filename_); - write_dataset(mesh_group, "library", "moab"); - // write volume of each tet - std::vector tet_vols; - xt::xtensor centroids({ehs_.size(), 3}); - for (int i = 0; i < ehs_.size(); i++) { - const auto& eh = ehs_[i]; - tet_vols.emplace_back(this->tet_volume(eh)); - Position c = this->centroid(eh); - xt::view(centroids, i, xt::all()) = xt::xarray({c.x, c.y, c.z}); - } - - write_dataset(mesh_group, "volumes", tet_vols); - write_dataset(mesh_group, "centroids", centroids); - - close_group(mesh_group); -} - bool -UnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) const { +MOABUnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) const { moab::ErrorCode rval; @@ -1928,7 +1934,7 @@ UnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) } int -UnstructuredMesh::get_bin_from_index(int idx) const { +MOABUnstructuredMesh::get_bin_from_index(int idx) const { if (idx >= n_bins()) { fatal_error(fmt::format("Invalid bin index: {}", idx)); } @@ -1936,25 +1942,25 @@ UnstructuredMesh::get_bin_from_index(int idx) const { } int -UnstructuredMesh::get_index(const Position& r, +MOABUnstructuredMesh::get_index(const Position& r, bool* in_mesh) const { int bin = get_bin(r); *in_mesh = bin != -1; return bin; } -int UnstructuredMesh::get_index_from_bin(int bin) const { +int MOABUnstructuredMesh::get_index_from_bin(int bin) const { return bin; } std::pair, std::vector> -UnstructuredMesh::plot(Position plot_ll, Position plot_ur) const { +MOABUnstructuredMesh::plot(Position plot_ll, Position plot_ur) const { // TODO: Implement mesh lines return {}; } int -UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { +MOABUnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { int bin = eh - ehs_[0]; if (bin >= n_bins()) { fatal_error(fmt::format("Invalid bin: {}", bin)); @@ -1963,18 +1969,18 @@ UnstructuredMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { } moab::EntityHandle -UnstructuredMesh::get_ent_handle_from_bin(int bin) const { +MOABUnstructuredMesh::get_ent_handle_from_bin(int bin) const { if (bin >= n_bins()) { fatal_error(fmt::format("Invalid bin index: ", bin)); } return ehs_[0] + bin; } -int UnstructuredMesh::n_bins() const { +int MOABUnstructuredMesh::n_bins() const { return ehs_.size(); } -int UnstructuredMesh::n_surface_bins() const { +int MOABUnstructuredMesh::n_surface_bins() const { // collect all triangles in the set of tets for this mesh moab::Range tris; moab::ErrorCode rval; @@ -1987,7 +1993,7 @@ int UnstructuredMesh::n_surface_bins() const { } Position -UnstructuredMesh::centroid(int bin) const { +MOABUnstructuredMesh::centroid(int bin) const { moab::ErrorCode rval; auto tet = this->get_ent_handle_from_bin(bin); @@ -2019,7 +2025,7 @@ UnstructuredMesh::centroid(int bin) const { } std::pair -UnstructuredMesh::get_score_tags(std::string score) const { +MOABUnstructuredMesh::get_score_tags(std::string score) const { moab::ErrorCode rval; // add a tag to the mesh // all scores are treated as a single value @@ -2061,7 +2067,7 @@ UnstructuredMesh::get_score_tags(std::string score) const { } void -UnstructuredMesh::add_score(const std::string& score) { +MOABUnstructuredMesh::add_score(const std::string& score) { auto score_tags = get_score_tags(score); } @@ -2094,7 +2100,7 @@ void UnstructuredMesh::remove_score(std::string score) const { } void -UnstructuredMesh::set_score_data(const std::string& score, +MOABUnstructuredMesh::set_score_data(const std::string& score, std::vector values, std::vector std_dev) { auto score_tags = this->get_score_tags(score); @@ -2126,7 +2132,7 @@ UnstructuredMesh::set_score_data(const std::string& score, } void -UnstructuredMesh::write(std::string base_filename) const { +MOABUnstructuredMesh::write(std::string base_filename) const { // add extension to the base name auto filename = base_filename + ".vtk"; write_message(5, "Writing unstructured mesh {}...", filename); @@ -2172,6 +2178,8 @@ LibMesh::centroid(int bin) const { return {centroid(0), centroid(1), centroid(2)}; } +std::string LibMesh::library() const { return "libmesh"; } + void LibMesh::initialize() { // always 3 for unstructured meshes n_dimension_ = 3; @@ -2358,23 +2366,7 @@ LibMesh::inside_tet(const libMesh::Point& r, return e->contains_point(r, FP_COINCIDENT); } -void LibMesh::to_hdf5(hid_t group) const -{ - hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); - - write_dataset(mesh_group, "type", "unstructured"); - write_dataset(mesh_group, "filename", filename_); - write_dataset(mesh_group, "library", "libmesh"); - - // write volume of each tet - std::vector tet_vols; - for (int i = 0; i < m_->n_elem(); i++) { - tet_vols.emplace_back(m_->elem_ref(i).volume()); - } - write_dataset(mesh_group, "volumes", tet_vols); - - close_group(mesh_group); -} +double LibMesh::volume(int bin) const { return m_->elem_ref(bin).volume(); } #endif // LIBMESH @@ -2407,7 +2399,7 @@ void read_meshes(pugi::xml_node root) #ifdef DAGMC } else if (mesh_type == "unstructured" && mesh_lib == "moab") { - model::meshes.push_back(std::make_unique(node)); + model::meshes.push_back(std::make_unique(node)); #else } else if (mesh_type == "unstructured") { fatal_error("Unstructured mesh support is disabled.");