diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index afe40fe433..8329bc5101 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -114,6 +114,9 @@ public:
   //! \param group_id An HDF5 group id.
   virtual void to_hdf5(hid_t group_id) const = 0;
 
+  //! Get the BoundingBox for this cell.
+  virtual BoundingBox bounding_box() const = 0;
+
   //----------------------------------------------------------------------------
   // Accessors
 
@@ -192,6 +195,8 @@ public:
 
   void to_hdf5(hid_t group_id) const;
 
+  BoundingBox bounding_box() const;
+
 protected:
   bool contains_simple(Position r, Direction u, int32_t on_surface) const;
   bool contains_complex(Position r, Direction u, int32_t on_surface) const;
diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 78e02d994b..3c5ece4f9e 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -44,12 +44,33 @@ namespace model {
 
 struct BoundingBox
 {
-  double xmin;
-  double xmax;
-  double ymin;
-  double ymax;
-  double zmin;
-  double zmax;
+  double xmin = INFTY;
+  double xmax = -INFTY;
+  double ymin = INFTY;
+  double ymax = -INFTY;
+  double zmin = INFTY;
+  double zmax = -INFTY;
+
+  // in-place update
+  inline void update(const BoundingBox& other) {
+    xmin = std::min(xmin, other.xmin);
+    xmax = std::max(xmax, other.xmax);
+    ymin = std::min(ymin, other.ymin);
+    ymax = std::max(ymax, other.ymax);
+    zmin = std::min(zmin, other.zmin);
+    zmax = std::max(zmax, other.zmax);
+  };
+
+  // in-place intersection
+  inline void intersect(const BoundingBox& other) {
+    xmin = std::max(xmin, other.xmin);
+    xmax = std::min(xmax, other.xmax);
+    ymin = std::max(ymin, other.ymin);
+    ymax = std::min(ymax, other.ymax);
+    zmin = std::max(zmin, other.zmin);
+    zmax = std::min(zmax, other.zmax);
+  }
+
 };
 
 //==============================================================================
@@ -107,7 +128,7 @@ public:
   virtual void to_hdf5(hid_t group_id) const = 0;
 
   //! Get the BoundingBox for this surface.
-  virtual BoundingBox bounding_box() const = 0;
+  virtual BoundingBox bounding_box(bool pos_side) const = 0;
 };
 
 class CSGSurface : public Surface
@@ -136,7 +157,7 @@ public:
   Direction normal(Position r) const;
   Direction reflect(Position r, Direction u) const;
   //! Get the bounding box of this surface.
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
 
   moab::DagMC* dagmc_ptr_;
   int32_t dag_index_;
@@ -185,7 +206,7 @@ public:
   void to_hdf5_inner(hid_t group_id) const;
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
 
   double x0_;
 };
@@ -206,7 +227,7 @@ public:
   void to_hdf5_inner(hid_t group_id) const;
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
 
   double y0_;
 };
@@ -227,7 +248,7 @@ public:
   void to_hdf5_inner(hid_t group_id) const;
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
 
   double z0_;
 };
@@ -248,7 +269,7 @@ public:
   void to_hdf5_inner(hid_t group_id) const;
   bool periodic_translate(const PeriodicSurface* other, Position& r,
                           Direction& u) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
 
   double A_, B_, C_, D_;
 };
@@ -268,7 +289,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   double y0_, z0_, radius_;
 };
 
@@ -287,7 +308,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   double x0_, z0_, radius_;
 };
 
@@ -306,7 +327,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   double x0_, y0_, radius_;
 };
 
@@ -325,7 +346,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   double x0_, y0_, z0_, radius_;
 };
 
@@ -344,7 +365,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   double x0_, y0_, z0_, radius_sq_;
 };
 
@@ -363,7 +384,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   double x0_, y0_, z0_, radius_sq_;
 };
 
@@ -382,7 +403,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   double x0_, y0_, z0_, radius_sq_;
 };
 
@@ -400,7 +421,7 @@ public:
   double distance(Position r, Direction u, bool coincident) const;
   Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  BoundingBox bounding_box() const;
+  BoundingBox bounding_box(bool pos_side) const;
   // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
   double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_;
 };
diff --git a/src/cell.cpp b/src/cell.cpp
index 04273eded0..9afae2d60c 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -567,6 +567,13 @@ CSGCell::to_hdf5(hid_t cell_group) const
   close_group(group);
 }
 
+BoundingBox CSGCell::bounding_box() const {
+  BoundingBox bbox;
+  for (int32_t token : rpn_) {
+    bbox.update(model::surfaces[abs(token)-1]->bounding_box(token > 0));
+  }
+}
+
 //==============================================================================
 
 bool
diff --git a/src/surface.cpp b/src/surface.cpp
index 7b610031f4..e9d2f964ec 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -271,7 +271,7 @@ Direction DAGSurface::reflect(Position r, Direction u) const
   return simulation::last_dir;
 }
 
-BoundingBox DAGSurface::bounding_box() const
+BoundingBox DAGSurface::bounding_box(bool pos_side) const
 {
   moab::ErrorCode rval;
   moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
@@ -366,9 +366,13 @@ bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other,
 }
 
 BoundingBox
-SurfaceXPlane::bounding_box() const
+SurfaceXPlane::bounding_box(bool pos_side) const
 {
-  return {x0_, x0_, -INFTY, INFTY, -INFTY, INFTY};
+  if (pos_side) {
+    return {x0_, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+  } else {
+    return {-INFTY, x0_, -INFTY, INFTY, -INFTY, INFTY};
+  }
 }
 
 //==============================================================================
@@ -428,9 +432,13 @@ bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other,
 }
 
 BoundingBox
-SurfaceYPlane::bounding_box() const
+SurfaceYPlane::bounding_box(bool pos_side) const
 {
-  return {-INFTY, INFTY, y0_, y0_, -INFTY, INFTY};
+  if (pos_side) {
+    return {-INFTY, INFTY, y0_, INFTY, -INFTY, INFTY};
+  } else {
+    return {-INFTY, INFTY, -INFTY, y0_, -INFTY, INFTY};
+  }
 }
 
 //==============================================================================
@@ -474,9 +482,13 @@ bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other,
 }
 
 BoundingBox
-SurfaceZPlane::bounding_box() const
+SurfaceZPlane::bounding_box(bool pos_side) const
 {
-  return {-INFTY, INFTY, -INFTY, INFTY, z0_, z0_};
+  if (pos_side) {
+    return {-INFTY, INFTY, -INFTY, INFTY, z0_, INFTY};
+  } else {
+    return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, z0_};
+  }
 }
 
 //==============================================================================
@@ -539,9 +551,34 @@ bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r,
 }
 
 BoundingBox
-SurfacePlane::bounding_box() const
+SurfacePlane::bounding_box(bool pos_side) const
 {
-  return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+  BoundingBox bbox = {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+  if (A_ == 0.0 && B_ == 0.0) {
+    double val = D_ / C_;
+    if (pos_side) {
+      bbox.zmin = val;
+    } else {
+      bbox.zmax = val;
+    }
+  } else if (A_ == 0.0 && C_ == 0.0) {
+    double val = D_ / B_;
+    if (pos_side) {
+    bbox.ymin = val;
+    } else {
+    bbox.ymax = val;
+    }
+  } else if (B_ == 0.0 && C_ == 0.0) {
+    double val = D_ / A_;
+    if (pos_side) {
+    bbox.xmin = val;
+    } else {
+    bbox.xmax = val;
+    }
+  }
+
+  return bbox;
+
 }
 
 //==============================================================================
@@ -652,8 +689,12 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceXCylinder::bounding_box() const {
-  return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_};
+BoundingBox SurfaceXCylinder::bounding_box(bool pos_side) const {
+  if (pos_side) {
+    return {-INFTY, INFTY, y0_ - radius_, y0_ + radius_, z0_ - radius_, z0_ + radius_};
+  } else {
+    return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+  }
 }
 //==============================================================================
 // SurfaceYCylinder implementation
@@ -688,8 +729,12 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceYCylinder::bounding_box() const {
-  return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_};
+BoundingBox SurfaceYCylinder::bounding_box(bool pos_side) const {
+  if (pos_side) {
+    return {x0_ - radius_, x0_ + radius_, -INFTY, INFTY, z0_ - radius_, z0_ + radius_};
+  } else {
+    return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+  }
 }
 
 //==============================================================================
@@ -725,8 +770,12 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceZCylinder::bounding_box() const {
-  return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY};
+BoundingBox SurfaceZCylinder::bounding_box(bool pos_side) const {
+  if (pos_side) {
+    return {x0_ - radius_, x0_ + radius_, y0_ - radius_, y0_ + radius_, -INFTY, INFTY};
+  } else {
+    return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+  }
 }
 
 
@@ -798,10 +847,14 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceSphere::bounding_box() const {
-  return {x0_ - radius_, x0_ + radius_,
-          y0_ - radius_, y0_ + radius_,
-          z0_ - radius_, z0_ + radius_};
+BoundingBox SurfaceSphere::bounding_box(bool pos_side) const {
+  if (pos_side) {
+    return {x0_ - radius_, x0_ + radius_,
+            y0_ - radius_, y0_ + radius_,
+            z0_ - radius_, z0_ + radius_};
+  } else {
+    return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
+  }
 }
 
 //==============================================================================
@@ -922,7 +975,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceXCone::bounding_box() const {
+BoundingBox SurfaceXCone::bounding_box(bool pos_side) const {
   return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
 }
 
@@ -959,7 +1012,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceYCone::bounding_box() const {
+BoundingBox SurfaceYCone::bounding_box(bool pos_side) const {
   return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
 }
 
@@ -996,7 +1049,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceZCone::bounding_box() const {
+BoundingBox SurfaceZCone::bounding_box(bool pos_side) const {
   return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
 }
 
@@ -1093,7 +1146,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-BoundingBox SurfaceQuadric::bounding_box() const {
+BoundingBox SurfaceQuadric::bounding_box(bool pos_side) const {
   return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
 }
 
@@ -1194,7 +1247,8 @@ void read_surfaces(pugi::xml_node node)
       }
 
       // See if this surface makes part of the global bounding box.
-      BoundingBox bb = surf->bounding_box();
+      BoundingBox bb = surf->bounding_box(true);
+      bb.intersect(surf->bounding_box(false));
       if (bb.xmin > -INFTY && bb.xmin < xmin) {
         xmin = bb.xmin;
         i_xmin = i_surf;