diff --git a/include/openmc/position.h b/include/openmc/position.h
index 2be72b27cb..3046f61636 100644
--- a/include/openmc/position.h
+++ b/include/openmc/position.h
@@ -63,6 +63,11 @@ struct Position {
     return std::sqrt(x*x + y*y + z*z);
   }
 
+  //! Reflect a direction across a normal vector
+  //! \param[in] other Vector to reflect across
+  //! \result Reflected vector
+  Position reflect(Position n) const;
+
   //! Rotate the position based on a rotation matrix
   Position rotate(const std::vector<double>& rotation) const;
 
@@ -89,6 +94,13 @@ inline Position operator/(Position a, Position b) { return a /= b; }
 inline Position operator/(Position a, double b)   { return a /= b; }
 inline Position operator/(double a, Position b)   { return b /= a; }
 
+inline Position Position::reflect(Position n) const {
+  const double projection = n.dot(*this);
+  const double magnitude = n.dot(n);
+  n *= (2.0 * projection / magnitude);
+  return *this - n;
+}
+
 inline bool operator==(Position a, Position b)
 {return a.x == b.x && a.y == b.y && a.z == b.z;}
 
diff --git a/src/cell.cpp b/src/cell.cpp
index 6c2bbf6bd7..6ca39749fd 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -779,7 +779,7 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
 // DAGMC Cell implementation
 //==============================================================================
 #ifdef DAGMC
-DAGCell::DAGCell() : Cell{} {};
+DAGCell::DAGCell() : Cell{} { simple_ = true; };
 
 std::pair<double, int32_t>
 DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const
diff --git a/src/surface.cpp b/src/surface.cpp
index f5a7f9f5f3..f2faa3efee 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -172,11 +172,9 @@ Surface::reflect(Position r, Direction u, Particle* p) const
   // Determine projection of direction onto normal and squared magnitude of
   // normal.
   Direction n = normal(r);
-  const double projection = n.dot(u);
-  const double magnitude = n.dot(n);
 
   // Reflect direction according to normal.
-  return u -= (2.0 * projection / magnitude) * n;
+  return u.reflect(n);
 }
 
 Direction
@@ -279,7 +277,13 @@ Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const
 {
   Expects(p);
   p->history_.reset_to_last_intersection();
-  p->last_dir_ = Surface::reflect(r, u, p);
+  moab::ErrorCode rval;
+  moab::EntityHandle surf = dagmc_ptr_->entity_by_index(2, dag_index_);
+  double pnt[3] = {r.x, r.y, r.z};
+  double dir[3];
+  rval = dagmc_ptr_->get_angle(surf, pnt, dir, &p->history_);
+  MB_CHK_ERR_CONT(rval);
+  p->last_dir_ = u.reflect(dir);
   return p->last_dir_;
 }