diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index eece0f1ff1..6bf49b2211 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -75,6 +75,12 @@ public:
 
   // Methods
 
+  //! Sample a mesh volume using a certain seed
+  //
+  //! \param[in] seed Seed to use for random sampling
+  //! \param[out] r Position within tet
+  virtual Position sample(uint64_t* seed) const=0;
+
   //! Determine which bins were crossed by a particle
   //
   //! \param[in] r0 Previous position of the particle
@@ -160,6 +166,8 @@ public:
     }
   };
 
+  Position sample(uint64_t* seed) const override;
+
   int get_bin(Position r) const override;
 
   int n_bins() const override;
@@ -460,7 +468,12 @@ public:
   static const std::string mesh_type;
   virtual std::string get_mesh_type() const override;
 
+  //! Sample a tetrahedron for an unstructured mesh
+  virtual Position sample_tet(std::array<Position, 4> coords, uint64_t* seed) const;
+
   // Overridden Methods
+  // TODO Position sample(uint64_t* seed) const=0;
+
   void surface_bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins) const override;
 
@@ -552,6 +565,8 @@ public:
 
   // Overridden Methods
 
+  Position sample(uint64_t* seed) const override;
+
   void bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins, vector<double>& lengths) const override;
 
@@ -714,6 +729,8 @@ public:
   void bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins, vector<double>& lengths) const override;
 
+  Position sample(uint64_t* seed) const override;
+
   int get_bin(Position r) const override;
 
   int n_bins() const override;
diff --git a/src/mesh.cpp b/src/mesh.cpp
index d9ad0ce8ad..191f0af5d6 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -195,6 +195,37 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
   }
 }
 
+// Sample barycentric coordinates given a seed and the vertex positions and return the sampled position
+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,”
+  if (s + t > 1) {
+      s = 1.0 - s;
+      t = 1.0 - t;
+  }
+  if (s + t + u > 1) {
+    if (t + u > 1) {
+      double old_t = t;
+      t = 1.0 - u;
+      u = 1.0 - s - old_t;
+    }else if (t + u <= 1) {
+      double old_s = s;
+      s = 1.0 - t - u;
+      u = old_s + t + u - 1;
+    }
+  }
+  return s*(coords[1]-coords[0]) 
+    + t*(coords[2]-coords[0]) 
+    + u*(coords[3]-coords[0]) 
+    + coords[0];
+
+}
+
 const std::string UnstructuredMesh::mesh_type = "unstructured";
 
 std::string UnstructuredMesh::get_mesh_type() const
@@ -327,6 +358,10 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const
   return ijk;
 }
 
+Position StructuredMesh::sample(uint64_t* seed) const {
+  fatal_error("Position sampling on structured meshes is not yet implemented");
+}
+
 int StructuredMesh::get_bin(Position r) const
 {
   // Determine indices
@@ -2049,6 +2084,36 @@ std::string MOABMesh::library() const
   return mesh_lib_type;
 }
 
+// Sample position within a tet for MOAB type tets
+Position MOABMesh::sample(uint64_t* seed) const {
+  // Get bin # assuming equal weights, IMP weigh by activity
+  int64_t tet_bin = round(n_bins()*prn(seed));
+
+  moab::EntityHandle tet_ent = get_ent_handle_from_bin(tet_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");
+  }
+  moab::CartVect p[4];
+  rval = mbi_->get_coords(conn1.data(), conn1.size(), p[0].array());
+  if (rval != moab::MB_SUCCESS) {
+    fatal_error("Failed to get tet coords");
+  }
+
+  std::array<Position, 4> tet_verts;
+  for (int i = 0; i < 4; i++) {
+    tet_verts[i] =  {p[i][0], p[i][1], p[i][2]};
+  }
+  // Samples position within tet using Barycentric stuff  
+  Position r = this->sample_tet(tet_verts, seed);
+
+  return r;
+}
+
+
 double MOABMesh::tet_volume(moab::EntityHandle tet) const
 {
   vector<moab::EntityHandle> conn;
@@ -2501,6 +2566,25 @@ void LibMesh::initialize()
   bbox_ = libMesh::MeshTools::create_bounding_box(*m_);
 }
 
+// Sample position within a tet for LibMesh type tets
+Position LibMesh::sample(uint64_t* seed) const {
+  // Get bin # assuming equal weights, IMP weigh by activity
+  int64_t tet_xi = round(n_bins()*prn(seed));
+
+  const auto& elem = get_element_from_bin(tet_xi);
+
+  // Get tet vertex coordinates from LibMesh
+  std::array<Position, 4> tet_verts;
+  for (int i = 0; i < elem.n_nodes(); i++) {
+    auto node_ref = elem.node_ref(i);
+    tet_verts[i] =  {node_ref(0), node_ref(1), node_ref(2)};
+  }
+  // Samples position within tet using Barycentric coordinates  
+  Position r = this->sample_tet(tet_verts, seed);
+
+  return r;
+}
+
 Position LibMesh::centroid(int bin) const
 {
   const auto& elem = this->get_element_from_bin(bin);