diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index da2fca2c8..c458caa7d 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -24,6 +24,7 @@
 #endif
 
 #ifdef LIBMESH
+#include "libmesh/bounding_box.h"
 #include "libmesh/libmesh.h"
 #include "libmesh/elem.h"
 #include "libmesh/equation_systems.h"
@@ -32,6 +33,7 @@
 #include "libmesh/dof_map.h"
 #include "libmesh/mesh.h"
 #include "libmesh/point.h"
+#include "libmesh/sphere.h"
 #endif
 
 namespace openmc {
@@ -578,7 +580,8 @@ private:
   std::vector<std::unique_ptr<libMesh::PointLocatorBase>> point_locators_; //!< pointers to locators for each thread
   std::unique_ptr<libMesh::EquationSystems> equation_systems_; //!< pointer to the equation systems of the mesh (for result output)
   std::map<std::string, unsigned int> variable_map_; //!< mappint of variable names (scores) to their numbers on the mesh
-  BoundingBox bbox_; //!< bounding box of the mesh
+  libMesh::BoundingBox bbox_; //!< bounding box of the mesh
+  libMesh::Sphere bsphere_; //<! bounding sphere of the mesh
   std::string eq_system_name_; //!< name of the equation system holding OpenMC results
   libMesh::Elem* first_element_; //!< pointer to the first element in the mesh (maybe should be a key?)
   std::set<libMesh::Elem*> boundary_elements_; //<! all boundary elements in the mesh
diff --git a/src/mesh.cpp b/src/mesh.cpp
index a3b36e3de..a3985de1d 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -28,6 +28,10 @@
 #include "openmc/tallies/filter.h"
 #include "openmc/xml_interface.h"
 
+#ifdef LIBMESH
+#include "libmesh/mesh_tools.h"
+#endif
+
 namespace openmc {
 
 //==============================================================================
@@ -2083,19 +2087,13 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMeshBase(node) {
   first_element_ = *m_->elements_begin();
 
   // bounding box for the mesh
-  bbox_ = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
+  bbox_ = libMesh::MeshTools::create_bounding_box(*m_);
+  bsphere_ = libMesh::MeshTools::bounding_sphere(*m_);
 
   // determine boundary elements and create bounding box
   for (int i = 0; i < m_->n_elem(); i++) {
     auto e = m_->elem_ptr(i);
 
-    // update bounding box for each node
-    for (int j = 0; j < e->n_nodes(); j++) {
-      auto n = e->node_ref(j);
-      Position r(n(0), n(1), n(2));
-      bbox_.update(r);
-    }
-
     for (int k = 0; k < e->n_neighbors(); k++) {
       if (!e->neighbor_ptr(k)) {
         boundary_elements_.insert(e);
@@ -2210,11 +2208,12 @@ LibMesh::bins_crossed(const Particle* p,
 int
 LibMesh::get_bin(Position r) const
 {
-  if (!bbox_.contains(r)) { return -1; }
-
   // look-up a tet using the point locator
   libMesh::Point p(r.x, r.y, r.z);
 
+  // quick rejection check
+  //  if (bsphere_.above_surface(p) || !bbox_.contains_point(p)) { return -1; }
+  if (!bbox_.contains_point(p)) { return -1; }
   int thread = omp_get_thread_num();
   auto e = (*point_locators_[thread])(p);