diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index ece7067590..939d055bad 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -113,7 +113,6 @@ public: int32_t n_sources() const { return this->mesh()->n_bins(); } private: - Mesh* mesh_ptr_ {nullptr}; int32_t mesh_idx_ {C_NONE}; double total_strength_ {0.0}; std::vector mesh_CDF_; diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3103204865..73a4f0bb74 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -78,7 +78,8 @@ public: //! Sample a mesh volume using a certain seed // //! \param[in] seed Seed to use for random sampling - //! \param[out] r Position within tet + //! \param[in] bin Bin value of the tet sampled + //! \return sampled position within tet virtual Position sample(uint64_t* seed, int32_t bin) const=0; //! Get the volume of a mesh bin @@ -536,6 +537,8 @@ public: true}; //!< Write tallies onto the unstructured mesh at the end of a run std::string filename_; //!< Path to unstructured mesh file + ElementType element_type(int bin) const; + protected: //! Set the length multiplier to apply to each point in the mesh void set_length_multiplier(const double length_multiplier); @@ -545,7 +548,12 @@ protected: 1.0}; //!< Constant multiplication factor to apply to mesh coordinates bool specified_length_multiplier_ {false}; - //! Sample a tetrahedron for an unstructured mesh + //! Sample barycentric coordinates given a seed and the vertex positions and + //! return the sampled position + // + //! \param[in] coords Coordinates of the tetrahedron + //! \param[in] seed Random number generation seed + //! \return Sampled position within the tetrahedron Position sample_tet(std::array coords, uint64_t* seed) const; private: diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index c1e000fcde..fbe5ce6583 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -194,9 +194,19 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) // Get pointer to spatial distribution mesh_idx_ = model::mesh_map.at(mesh_id); - // Check whether mesh pointer points to a mesh - mesh_ptr_ = dynamic_cast(model::meshes[mesh_idx_].get()); - if (!mesh_ptr_) {fatal_error("Mesh passed to spatial distribution is not a mesh object"); } + auto mesh_ptr = + dynamic_cast(model::meshes.at(mesh_idx_).get()); + if (!mesh_ptr) { + fatal_error("Only unstructured mesh is supported for source sampling."); + } + + // ensure that the unstructured mesh contains only linear tets + for (int bin = 0; bin < mesh_ptr->n_bins(); bin++) { + if (mesh_ptr->element_type(bin) != ElementType::LINEAR_TET) { + fatal_error( + "Mesh specified for source must contain only linear tetrahedra."); + } + } int32_t n_bins = this->n_sources(); std::vector strengths(n_bins, 0.0); @@ -222,7 +232,7 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) if (get_node_value_bool(node, "volume_normalized")) { for (int i = 0; i < n_bins; i++) { - mesh_strengths_[i] *= mesh_ptr_->volume(i); + mesh_strengths_[i] *= mesh()->volume(i); } } @@ -244,7 +254,7 @@ Position MeshSpatial::sample(uint64_t* seed) const double eta = prn(seed); // Sample over the CDF defined in initialization above int32_t elem_idx = lower_bound_index(mesh_CDF_.begin(), mesh_CDF_.end(), eta); - return mesh_ptr_->sample(seed, elem_idx); + return mesh()->sample(seed, elem_idx); } //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 2102f953de..39b5a225e8 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -195,14 +195,17 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) } } -Position UnstructuredMesh::sample_tet(std::array coords, uint64_t* seed) const { +Position UnstructuredMesh::sample_tet( + std::array coords, uint64_t* seed) const +{ // Uniform distribution double s = prn(seed); double t = prn(seed); double u = prn(seed); - // From PyNE implementation of moab tet sampling - // C. Rocchini and P. Cignoni, “Generating Random Points in a Tetrahedron,” + // From PyNE implementation of moab tet sampling C. Rocchini & P. Cignoni + // (2000) Generating Random Points in a Tetrahedron, Journal of Graphics + // Tools, 5:4, 9-12, DOI: 10.1080/10867651.2000.10487528 if (s + t > 1) { s = 1.0 - s; t = 1.0 - t; @@ -222,7 +225,6 @@ Position UnstructuredMesh::sample_tet(std::array coords, uint64_t* + t*(coords[2]-coords[0]) + u*(coords[3]-coords[0]) + coords[0]; - } const std::string UnstructuredMesh::mesh_type = "unstructured"; @@ -313,6 +315,18 @@ void UnstructuredMesh::set_length_multiplier(double length_multiplier) specified_length_multiplier_ = true; } +ElementType UnstructuredMesh::element_type(int bin) const +{ + auto conn = connectivity(bin); + + if (conn.size() == 4) + return ElementType::LINEAR_TET; + else if (conn.size() == 8) + return ElementType::LINEAR_HEX; + else + return ElementType::UNSUPPORTED; +} + StructuredMesh::MeshIndex StructuredMesh::get_indices( Position r, bool& in_mesh) const { @@ -2095,10 +2109,14 @@ Position MOABMesh::sample(uint64_t* seed, int32_t bin) const { moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin); // Get vertex coordinates for MOAB tet - vector conn1; - moab::ErrorCode rval = mbi_->get_connectivity(&tet_ent, 1, conn1); - if (rval != moab::MB_SUCCESS) { - fatal_error("Failed to get tet connectivity"); + std::array conn1; + int conn1_size; + moab::ErrorCode rval = + mbi_->get_connectivity(&tet_ent, 1, conn1.data(), conn1_size); + if (rval != moab::MB_SUCCESS || conn1_size != 4) { + fatal_error(fmt::format( + "Failed to get tet connectivity or connectivity size () is invalid.", + conn1_size)); } moab::CartVect p[4]; rval = mbi_->get_coords(conn1.data(), conn1.size(), p[0].array());