From 80077ed3fd04758a1fb10095648e977b515a5305 Mon Sep 17 00:00:00 2001
From: NybergWISC <mnyberg@wisc.edu>
Date: Tue, 21 Jun 2022 12:06:45 -0500
Subject: [PATCH] New changes to spatial distribution and source to include
 mesh sources. Made some edits suggested by Shriwise, now compiles
 successfully

---
 include/openmc/distribution_spatial.h | 19 +++++++++++++++++++
 src/distribution_spatial.cpp          | 25 +++++++++++++++++++++++++
 src/mesh.cpp                          |  6 ++++--
 src/source.cpp                        |  2 ++
 4 files changed, 50 insertions(+), 2 deletions(-)

diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h
index 5e426b878..8a77a4bd3 100644
--- a/include/openmc/distribution_spatial.h
+++ b/include/openmc/distribution_spatial.h
@@ -5,6 +5,7 @@
 
 #include "openmc/distribution.h"
 #include "openmc/position.h"
+#include "openmc/mesh.h"
 
 namespace openmc {
 
@@ -94,6 +95,24 @@ private:
   Position origin_;    //!< Cartesian coordinates of the sphere center
 };
 
+//==============================================================================
+//! Distribution of points within a mesh
+//==============================================================================
+
+class MeshIndependent : public SpatialDistribution {
+public:
+  explicit MeshIndependent(pugi::xml_node node);
+
+  //! Sample a position from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled position
+  Position sample(uint64_t* seed) const;
+
+private:
+  UnstructuredMesh* umesh_ptr_;
+  Position sampled_;
+};
+
 //==============================================================================
 //! Uniform distribution of points over a box
 //==============================================================================
diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp
index a06f84d80..f78c2c49b 100644
--- a/src/distribution_spatial.cpp
+++ b/src/distribution_spatial.cpp
@@ -180,6 +180,31 @@ Position SphericalIndependent::sample(uint64_t* seed) const
   return {x, y, z};
 }
 
+//==============================================================================
+// MeshIndependent implementation
+//==============================================================================
+
+// Loosely adapted Patrick's code shown here
+MeshIndependent::MeshIndependent(pugi::xml_node node)
+{
+  // No in-tet distributions implemented, could include distributions for the barycentric coords
+  
+  int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
+  const auto& mesh_ptr = model::meshes[mesh_id];
+
+  // line from Patrick code, unsure of use
+  // int bin = mesh_ptr->get_bin();
+  const UnstructuredMesh* umesh_ptr_ = dynamic_cast<UnstructuredMesh*>(mesh_ptr.get());
+  if (!umesh_ptr_) {fatal_error("Mesh passed to spatial distribution is not an unstructured mesh object"); }
+}
+
+Position MeshIndependent::sample(uint64_t* seed) const
+{ 
+  Position sampled_ = umesh_ptr_->sample(seed);
+  return sampled_;
+}
+
+
 //==============================================================================
 // SpatialBox implementation
 //==============================================================================
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 191f0af5d..56d793079 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -2087,7 +2087,8 @@ std::string MOABMesh::library() const
 // 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));
+  // This may underweigh first and last bins
+  int64_t tet_bin = trunc(n_bins()*prn(seed));
 
   moab::EntityHandle tet_ent = get_ent_handle_from_bin(tet_bin);
 
@@ -2569,7 +2570,8 @@ void LibMesh::initialize()
 // 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));
+  // This may underweight first and last bin by 1/2 TODO
+  int64_t tet_xi = trunc(n_bins()*prn(seed));
 
   const auto& elem = get_element_from_bin(tet_xi);
 
diff --git a/src/source.cpp b/src/source.cpp
index 3a80e1156..58aac711f 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -94,6 +94,8 @@ IndependentSource::IndependentSource(pugi::xml_node node)
         space_ = UPtrSpace {new CylindricalIndependent(node_space)};
       } else if (type == "spherical") {
         space_ = UPtrSpace {new SphericalIndependent(node_space)};
+      } else if (type =="mesh"){
+        space_ = UPtrSpace {new MeshIndependent(node_space)};
       } else if (type == "box") {
         space_ = UPtrSpace {new SpatialBox(node_space)};
       } else if (type == "fission") {