From 2d40130095b337552bf47249ce3292fa75f6cf7c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 17 Feb 2020 11:43:36 -0600 Subject: [PATCH] Factoring out a new SructuredMesh class. Plot method placement still needs to be resolved. Updating unstructured mesh to allow use of indices only. Going back to a single unstructure mesh class until we need the added abstraction. --- include/openmc/mesh.h | 94 +++++++++++++++++++++++++------------ src/mesh.cpp | 68 ++++++++++++++++++++------- src/tallies/filter_mesh.cpp | 13 +---- 3 files changed, 115 insertions(+), 60 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 8a66e06ab..8494e8578 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -69,6 +69,43 @@ public: //! \return Mesh bin virtual int get_bin(Position r) const = 0; + //! Get the number of mesh cells. + virtual int n_bins() const = 0; + + //! Get the number of mesh cell surfaces. + virtual int n_surface_bins() const = 0; + + //! Write mesh data to an HDF5 group + // + //! \param[in] group HDF5 group + virtual void to_hdf5(hid_t group) const = 0; + + //! Find the mesh lines that intersect an axis-aligned slice plot + // + //! \param[in] plot_ll The lower-left coordinates of the slice plot. + //! \param[in] plot_ur The upper-right coordinates of the slice plot. + //! \return A pair of vectors indicating where the mesh lines lie along each + //! of the plot's axes. For example an xy-slice plot will get back a vector + //! of x-coordinates and another of y-coordinates. These vectors may be + //! empty for low-dimensional meshes. + virtual std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const = 0; + + //! Get a label for the mesh bin + virtual std::string bin_label(int bin) const = 0; + + // Data members + int id_ {-1}; //!< User-specified ID + int n_dimension_; //!< Number of dimensions +}; + +class StructuredMesh : public Mesh +{ +public: + StructuredMesh() = default; + StructuredMesh(pugi::xml_node node) : Mesh {node} {}; + virtual ~StructuredMesh() = default; + //! Get bin given mesh indices // //! \param[in] Array of mesh indices @@ -88,32 +125,21 @@ public: //! \param[out] ijk Mesh indices virtual void get_indices_from_bin(int bin, int* ijk) const = 0; - //! Get the number of mesh cells. - virtual int n_bins() const = 0; + //! Get a label for the mesh bin + std::string bin_label(int bin) const override { + std::vector ijk(n_dimension_); + get_indices_from_bin(bin, ijk.data()); - //! Get the number of mesh cell surfaces. - virtual int n_surface_bins() const = 0; - - //! Find the mesh lines that intersect an axis-aligned slice plot - // - //! \param[in] plot_ll The lower-left coordinates of the slice plot. - //! \param[in] plot_ur The upper-right coordinates of the slice plot. - //! \return A pair of vectors indicating where the mesh lines lie along each - //! of the plot's axes. For example an xy-slice plot will get back a vector - //! of x-coordinates and another of y-coordinates. These vectors may be - //! empty for low-dimensional meshes. - virtual std::pair, std::vector> - plot(Position plot_ll, Position plot_ur) const = 0; - - //! Write mesh data to an HDF5 group - // - //! \param[in] group HDF5 group - virtual void to_hdf5(hid_t group) const = 0; + if (n_dim > 2) { + return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]); + } else if (n_dim > 1) { + return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]); + } else { + return fmt::format("Mesh Index ({})", ijk[0]) ; + } + } // Data members - - int id_ {-1}; //!< User-specified ID - int n_dimension_; //!< Number of dimensions xt::xtensor lower_left_; //!< Lower-left coordinates of mesh xt::xtensor upper_right_; //!< Upper-right coordinates of mesh }; @@ -122,7 +148,7 @@ public: //! Tessellation of n-dimensional Euclidean space by congruent squares or cubes //============================================================================== -class RegularMesh : public Mesh +class RegularMesh : public StructuredMesh { public: // Constructors @@ -188,7 +214,7 @@ private: }; -class RectilinearMesh : public Mesh +class RectilinearMesh : public StructuredMesh { public: // Constructors @@ -245,6 +271,7 @@ class UnstructuredMesh : public Mesh { public: UnstructuredMesh() { }; UnstructuredMesh(pugi::xml_node); + ~UnstructuredMesh() = default; //! Determine which bins were crossed by a particle. // @@ -316,14 +343,11 @@ 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; + int get_bin_from_index(int idx) const; - void get_indices(Position r, int* ijk, bool* in_mesh) const override; + int get_index(Position r, bool* in_mesh) const; - void get_indices_from_bin(int bin, int* ijk) const override; - - std::pair, std::vector> - plot(Position plot_ll, Position plot_ur) const override; + int get_index_from_bin(int bin) const; //! Builds a KDTree for all tetrahedra in the mesh. All //! triangles representing 2D faces of the mesh are @@ -333,6 +357,10 @@ intersect_track(const moab::CartVect& start, void build_kdtree(const moab::Range& all_tets); public: + + std::pair, std::vector> + plot(Position plot_ll, Position plot_ur) const override; + //! Determine which surface bins were crossed by a particle. // //! \param[in] p Particle to check @@ -356,6 +384,10 @@ public: int n_surface_bins() const override; + Position centroid(moab::EntityHandle tet) const; + + std::string bin_label(int bin) const override; + std::string filename_; //(new moab::Core()); + mbi_ = std::unique_ptr(new moab::Core()); // create meshset to load mesh into moab::ErrorCode rval = mbi_->create_meshset(moab::MESHSET_SET, meshset_); if (rval != moab::MB_SUCCESS) { @@ -1855,7 +1855,6 @@ UnstructuredMesh::compute_barycentric_data(const moab::Range& all_tets) { a = a.transpose().inverse(); baryc_data_.at(get_bin_from_ent_handle(tet)) = a; } - } // TODO: write this function @@ -1911,25 +1910,24 @@ UnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet) } int -UnstructuredMesh::get_bin_from_indices(const int* ijk) const { - if (ijk[0] >= n_bins()) { +UnstructuredMesh::get_bin_from_index(int idx) const { + if (idx >= n_bins()) { std::stringstream s; - s << "Invalid bin: " << ijk[0]; + s << "Invalid bin index: " << idx; fatal_error(s); } - int bin = ehs_[ijk[0]] - ehs_[0]; + return ehs_[idx] - ehs_[0]; +} + +int +UnstructuredMesh::get_index(Position r, bool* in_mesh) const { + int bin = get_bin(r); + *in_mesh = bin != -1; 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; +int UnstructuredMesh::get_index_from_bin(int bin) const { + return bin; } std::pair, std::vector> @@ -1972,6 +1970,42 @@ int UnstructuredMesh::n_surface_bins() const { return 2 * tris.size(); } +Position +UnstructuredMesh::centroid(moab::EntityHandle tet) const { + moab::ErrorCode rval; + + // look up the tet connectivity + std::vector conn; + rval = mbi_->get_connectivity(&tet, 1, conn); + if (rval != moab::MB_SUCCESS) { + warning("Failed to get connectivity of a mesh element."); + return {}; + } + + // get the coordinates + std::vector coords(conn.size()); + rval = mbi_->get_coords(&conn.front(), conn.size(), coords[0].array()); + if (rval != moab::MB_SUCCESS) { + warning("Failed to get the coordinates of a mesh element."); + return {}; + } + + // compute the centroid of the elements + moab::CartVect centroid(0.0); + for(const auto& coord : coords) { + centroid += coord; + } + centroid /= double(coords.size()); + + return {centroid[0], centroid[1], centroid[2]}; +} + +std::string +UnstructuredMesh::bin_label(int bin) const { + std::stringstream out; + out << "Mesh Index (" << bin << ")"; +}; + #endif diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 6c2cc8bd1..d2a240bcf 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -55,18 +55,7 @@ std::string MeshFilter::text_label(int bin) const { auto& mesh = *model::meshes[mesh_]; - int n_dim = mesh.n_dimension_; - - std::vector ijk(n_dim); - mesh.get_indices_from_bin(bin, ijk.data()); - - if (n_dim > 2) { - return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]); - } else if (n_dim > 1) { - return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]); - } else { - return fmt::format("Mesh Index ({})", ijk[0]) ; - } + return mesh.bin_label(bin); } void