Addressing comments from @paulromano's review

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<double> mesh_CDF_;

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<Position, 4> coords, uint64_t* seed) const;
 
 private:

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<Mesh*>(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<UnstructuredMesh*>(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<double> 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);
 }
 
 //==============================================================================

src/mesh.cpp

@@ -195,14 +195,17 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
   }
 }
 
-Position UnstructuredMesh::sample_tet(std::array<Position, 4> coords, uint64_t* seed) const {
+Position UnstructuredMesh::sample_tet(
+  std::array<Position, 4> 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<Position, 4> 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<moab::EntityHandle> 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<moab::EntityHandle, 4> 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());