diff --git a/openmc/geometry.py b/openmc/geometry.py
index d682d7ee2..a9b222d11 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -270,35 +270,17 @@ class Geometry(object):
         """
         surfaces = OrderedDict()
         
-        for cell in self.get_all_cells().values():
-            self._get_surfaces_from_region(surfaces, cell.region)
+        root_cells = self._root_universe.cells
+        for cell in root_cells.values():
+            reg = cell.region
+            surfaces = reg.get_surfaces_from_region(surfaces)
+            
+        #surfaces = self._root_universe.
+        #for cell in self.get_all_cells().values():
+        #    self.get_surfaces_from_region(surfaces, cell.region)
         return surfaces
     
-    def _get_surfaces_from_region(self, surfaces, region):
-        """
-        Recursively find all the surfaces referenced by a region and return them
-        
-        Parameters
-        ----------
-        surfaces: collections.OrderedDict
-            Dictionary mapping surface IDs to :class:`openmc.Surface` instances
-            
-        region: openmc.surface.Region
-            The region of space defined by Surfaces
-        Returns
-        -------
-        collections.OrderedDict
-            Dictionary mapping surface IDs to :class:`openmc.Surface` instances
-
-        """
-        if isinstance(region, openmc.Halfspace):
-            s = region.surface
-            if s.id not in surfaces:
-                surfaces[s.id] = s
-        else:
-            for reg in region:
-                surfaces = self._get_surfaces_from_region(surfaces, reg)
-        return surfaces
+    
                 
     def get_materials_by_name(self, name, case_sensitive=False, matching=False):
         """Return a list of materials with matching names.
diff --git a/openmc/region.py b/openmc/region.py
index 4494a3f62..012faca75 100644
--- a/openmc/region.py
+++ b/openmc/region.py
@@ -1,5 +1,5 @@
 from abc import ABCMeta, abstractmethod
-from collections import Iterable
+from collections import Iterable, OrderedDict
 
 from six import add_metaclass
 import numpy as np
@@ -46,6 +46,25 @@ class Region(object):
     def __ne__(self, other):
         return not self == other
 
+    def get_surfaces_from_region(self, surfaces = OrderedDict()):
+        """
+        Recursively find all the surfaces referenced by a region and return them
+
+        Parameters
+        ----------
+        surfaces: collections.OrderedDict, optional
+            Dictionary mapping surface IDs to :class:`openmc.Surface` instances
+            
+        Returns
+        -------
+        surfaces: collections.OrderedDict
+            Dictionary mapping surface IDs to :class:`openmc.Surface` instances
+
+        """
+        for region in self:
+            surfaces = region.get_surfaces_from_region(surfaces)
+        return surfaces
+
     @staticmethod
     def from_expression(expression, surfaces):
         """Generate a region given an infix expression.
diff --git a/openmc/surface.py b/openmc/surface.py
index 4f20c856c..894f2f6a1 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -1,5 +1,5 @@
 from abc import ABCMeta
-from collections import Iterable
+from collections import Iterable, OrderedDict
 from numbers import Real, Integral
 from xml.etree import ElementTree as ET
 from math import sqrt
@@ -1880,6 +1880,26 @@ class Halfspace(Region):
     def __str__(self):
         return '-' + str(self.surface.id) if self.side == '-' \
             else str(self.surface.id)
+    
+    def get_surfaces_from_region(self, surfaces = OrderedDict()):
+        """
+        Returns the surface that this is a halfspace of.
+        Overwrites method Region.get_surfaces_from_region()
+    
+        Parameters
+        ----------
+        surfaces: collections.OrderedDict, optional
+            Dictionary mapping surface IDs to :class:`openmc.Surface` instances
+    
+        Returns
+        -------
+        surfaces: collections.OrderedDict
+            Dictionary mapping surface IDs to :class:`openmc.Surface` instances
+    
+        """
+        if self.surface.id not in surfaces:
+            surfaces[self.surface.id] = self.surface
+        return surfaces
 
 
 def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),