@@ -110,7 +110,13 @@ class MeshIndependent : public SpatialDistribution {
 
 private:
   UnstructuredMesh* umesh_ptr_;
-  Position sampled_;
+  int32_t mesh_map_idx_;
+  std::string sample_scheme_;
+  double total_weight_;
+  std::vector<double> mesh_CDF_;
+  std::vector<double> mesh_weights_;
+  int64_t tot_bins_;
+  
 };
 
 //==============================================================================
 
@@ -79,7 +79,7 @@ class Mesh {
   //
   //! \param[in] seed Seed to use for random sampling
   //! \param[out] r Position within tet
-  virtual Position sample(uint64_t* seed) const=0;
+  virtual Position sample(uint64_t* seed, int32_t tet_bin) const=0;
 
   //! Determine which bins were crossed by a particle
   //
@@ -166,7 +166,7 @@ class StructuredMesh : public Mesh {
     }
   };
 
-  Position sample(uint64_t* seed) const override;
+  Position sample(uint64_t* seed, int32_t tet_bin) const override;
 
   int get_bin(Position r) const override;
 
@@ -565,7 +565,7 @@ class MOABMesh : public UnstructuredMesh {
 
   // Overridden Methods
 
-  Position sample(uint64_t* seed) const override;
+  Position sample(uint64_t* seed, int32_t tet_bin) const override;
 
   void bins_crossed(Position r0, Position r1, const Direction& u,
     vector<int>& bins, vector<double>& lengths) const override;
@@ -729,7 +729,7 @@ class LibMesh : public UnstructuredMesh {
   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;
+  Position sample(uint64_t* seed, int32_t tet_bin) const override;
 
   int get_bin(Position r) const override;
 
 
@@ -11,10 +11,12 @@
 from xml.etree import ElementTree as ET
 
 import openmc.checkvalue as cv
+from openmc.stats.multivariate import MeshIndependent
 
 from . import RegularMesh, Source, VolumeCalculation, WeightWindows
 from ._xml import clean_indentation, get_text, reorder_attributes
 from openmc.checkvalue import PathLike
+from .mesh import MeshBase
 
 
 class RunMode(Enum):
@@ -975,6 +977,8 @@ def _create_max_order_subelement(self, root):
     def _create_source_subelement(self, root):
         for source in self.source:
             root.append(source.to_xml_element())
+            if isinstance(source.space, MeshIndependent):
+                root.append(source.space.mesh.to_xml_element())
 
     def _create_volume_calcs_subelement(self, root):
         for calc in self.volume_calculations:
@@ -1323,6 +1327,12 @@ def _keff_trigger_from_xml_element(self, root):
     def _source_from_xml_element(self, root):
         for elem in root.findall('source'):
             self.source.append(Source.from_xml_element(elem))
+            if isinstance(Source.from_xml_element(elem), MeshIndependent):
+                mesh_id = int(get_text(elem, 'mesh'))
+                path = f"./mesh[@id='{mesh_id}']"
+                mesh_elem = root.find(path)
+                if mesh_elem is not None:
+                    self.source.space.mesh = MeshBase.from_xml_element(mesh_elem)
 
     def _volume_calcs_from_xml_element(self, root):
         volume_elems = root.findall("volume_calc")
 
@@ -16,6 +16,7 @@
 from openmc.stats.multivariate import UnitSphere, Spatial
 from openmc.stats.univariate import Univariate
 from ._xml import get_text
+from .mesh import MeshBase
 
 
 class Source:
 
@@ -3,12 +3,14 @@
 from math import pi, cos
 from numbers import Real
 from xml.etree import ElementTree as ET
+from xmlrpc.client import boolean
 
 import numpy as np
 
 import openmc.checkvalue as cv
 from .._xml import get_text
 from .univariate import Univariate, Uniform, PowerLaw
+from ..mesh import MeshBase
 
 
 class UnitSphere(ABC):
@@ -269,6 +271,8 @@ def from_xml_element(cls, elem):
             return CylindricalIndependent.from_xml_element(elem)
         elif distribution == 'spherical':
             return SphericalIndependent.from_xml_element(elem)
+        elif distribution == 'mesh':
+            return MeshIndependent.from_xml_element(elem)
         elif distribution == 'box' or distribution == 'fission':
             return Box.from_xml_element(elem)
         elif distribution == 'point':
@@ -616,6 +620,120 @@ def from_xml_element(cls, elem):
         origin = [float(x) for x in elem.get('origin').split()]
         return cls(r, phi, z, origin=origin)
 
+class MeshIndependent(Spatial):
+    """Spatial distribution for a mesh.
+
+    This distribution allows one to specify a mesh to sample over and the scheme to sample over the entire mesh. 
+
+    .. versionadded:: 0.13
+
+    Parameters
+    ----------
+    elem_weight_scheme : str, optional
+        The scheme for weighting and sampling elements from the mesh. Options are 'file' and 'volume' based weights.
+    mesh : openmc.MeshBase
+        The mesh instance used for sampling, mesh is written into settings.xml, mesh.id is written into the source distribution
+    weights_from_file : Iterable of Real, optional
+        A list of values which represent the weights of each element
+
+
+    Attributes
+    ----------
+    elem_weight_scheme : str, optional
+        Weighting scheme for sampling element from mesh, defaults to volume
+    mesh : openmc.MeshBase
+        The mesh instance used for sampling, mesh is written into settings.xml, mesh.id is written into the source distribution
+    weights_from_file : Iterable of Real, optional
+        A list of values which represent the weights of each element
+
+    """
+
+    def __init__(self, mesh, weights_from_file=None, volume_weighting=False):
+        self.mesh = mesh        
+        self.weights_from_file = weights_from_file
+        self.volume_weighting = volume_weighting
+
+    @property
+    def mesh(self):
+        return self._mesh
+
+    @mesh.setter
+    def mesh(self, mesh):
+        if mesh != None:
+            cv.check_type('Unstructured Mesh', mesh, MeshBase)
+        self._mesh = mesh
+
+    @property
+    def volume_weighting(self):
+        return self._volume_weighting
+
+    @volume_weighting.setter
+    def volume_weighting(self, volume_weighting):
+        cv.check_type('Scheme for sampling an element from the mesh', volume_weighting, bool)
+        self._volume_weighting = volume_weighting
+
+    @property
+    def weights_from_file(self):
+        return self._weights_from_file
+
+    @weights_from_file.setter
+    def weights_from_file(self, given_weights):
+        if given_weights is not None:
+            self._weights_from_file = (np.array(given_weights, dtype=float)).flatten()
+        else:
+            self._weights_from_file = None
+
+    @property
+    def num_weight_bins(self):
+        if self.weights_from_file == None:
+            raise ValueError('Weight bins are not set')
+        return self.weights_from_file.size
+                
+    def to_xml_element(self):
+        """Return XML representation of the spatial distribution
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing spatial distribution data
+
+        """
+        element = ET.Element('space')
+        element.set('type', 'mesh')
+        element.set("mesh_id", str(self.mesh.id))
+        element.set("volume_weighting", str(self.volume_weighting))
+
+        if self.weights_from_file is not None:
+            subelement = ET.SubElement(element, 'weights_from_file')
+            subelement.text = ' '.join(str(e) for e in self.weights_from_file)
+
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate spatial distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.MeshIndependent
+            Spatial distribution generated from XML element
+
+        """
+
+        mesh_id = int(elem.get('mesh_id'))
+        volume_weighting = elem.get("volume_weighting")
+        volume_weighting = get_text(elem, 'volume_weighting').lower() == 'true'
+        if elem.get('weights_from_file') is not None:
+            weights_from_file = [float(b) for b in get_text(elem, 'weights_from_file').split()]
+        else:
+            weights_from_file = None
+        return cls(volume_weighting, mesh_id, weights_from_file)
+
 
 class Box(Spatial):
     """Uniform distribution of coordinates in a rectangular cuboid.
 
@@ -3,6 +3,11 @@
 #include "openmc/error.h"
 #include "openmc/random_lcg.h"
 #include "openmc/xml_interface.h"
+#include "openmc/mesh.h"
+#include "openmc/search.h"
+
+#include <iostream>
+#include <fstream>
 
 namespace openmc {
 
@@ -184,24 +189,71 @@ Position SphericalIndependent::sample(uint64_t* seed) const
 // 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
+  // Read in unstructured mesh from mesh_id value
+  int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id"));
+  // Get pointer to spatial distribution
+  mesh_map_idx_ = model::mesh_map.at(mesh_id);
+  const auto& mesh_ptr = model::meshes[mesh_map_idx_];
 
-  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());
+  // Check whether mesh pointer points to an unstructured mesh
+  umesh_ptr_ = dynamic_cast<UnstructuredMesh*>(mesh_ptr.get());
   if (!umesh_ptr_) {fatal_error("Mesh passed to spatial distribution is not an unstructured mesh object"); }
+
+  // Create CDF based on weighting scheme
+  tot_bins_ = umesh_ptr_->n_bins();
+  std::vector<double> weights = {};
+  weights.resize(tot_bins_);
+  double temp_total_weight = 0.0;
+  mesh_CDF_.resize(tot_bins_+1);
+  mesh_CDF_[0] = {0.0};
+  total_weight_ = 0.0;
+  mesh_weights_.resize(tot_bins_);
+
+  // Create cdfs for sampling for an element over a mesh
+  // Volume scheme is weighted by the volume of each tet
+  // File scheme is weighted by an array given in the xml file
+
+  if (check_for_node(node, "weights_from_file")) {
+    mesh_weights_ = get_node_array<double>(node, "weights_from_file");
+    if (mesh_weights_.size() != tot_bins_){
+      write_message("WARNING: The size of the weights array from the xml file does not equal the number of elements in the mesh.");
+    } if (get_node_value_bool(node, "volume_weighting")){
+      for (int i = 0; i < tot_bins_; i++){
+        weights[i] = mesh_weights_[i]*umesh_ptr_->volume(i);
+      } 
+    } else if (!get_node_value_bool(node, "volume_weighting")){
+      for (int i = 0; i < tot_bins_; i++){
+        weights[i] = mesh_weights_[i];  
+      }
+    }
+  } else if (get_node_value_bool(node, "volume_weighting")){
+    for (int i = 0; i<tot_bins_; i++){
+      weights[i] = umesh_ptr_->volume(i);
+    }
+  } else {
+    for (int i = 0; i<tot_bins_; i++){
+      weights[i] = 1;
+    }
+  }
+  for (int i = 0; i<tot_bins_; i++){
+    temp_total_weight = temp_total_weight + weights[i];
+  }
+  total_weight_ = temp_total_weight;
+  for (int i = 0; i<tot_bins_; i++){
+    mesh_CDF_[i+1] = mesh_CDF_[i] + weights[i]/total_weight_;
+  }
 }
 
 Position MeshIndependent::sample(uint64_t* seed) const
 { 
-  Position sampled_ = umesh_ptr_->sample(seed);
-  return sampled_;
+  // Create random variable for sampling element from mesh
+  float eta = prn(seed);
+  // Sample over the CDF defined in initialization above
+  int32_t tet_bin = lower_bound_index(mesh_CDF_.begin(), mesh_CDF_.end(), eta);
+  return umesh_ptr_->sample(seed, tet_bin);
 }
 
 
 
@@ -405,7 +405,6 @@ void read_separate_xml_files()
 
   // Finalize cross sections having assigned temperatures
   finalize_cross_sections();
-
   read_tallies_xml();
 
   // Initialize distribcell_filters
 
@@ -358,7 +358,7 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const
   return ijk;
 }
 
-Position StructuredMesh::sample(uint64_t* seed) const {
+Position StructuredMesh::sample(uint64_t* seed, int32_t tet_bin) const {
   fatal_error("Position sampling on structured meshes is not yet implemented");
 }
 
@@ -2085,10 +2085,7 @@ 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
-  // This may underweigh first and last bins
-  int64_t tet_bin = trunc(n_bins()*prn(seed));
+Position MOABMesh::sample(uint64_t* seed, int32_t tet_bin) const {
 
   moab::EntityHandle tet_ent = get_ent_handle_from_bin(tet_bin);
 
@@ -2568,18 +2565,13 @@ 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
-  // 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);
-
+Position LibMesh::sample(uint64_t* seed, int32_t tet_bin) const {
+  const auto& elem = get_element_from_bin(tet_bin);
   // 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)};
+    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);