Implementing duplicate surface detection and removal

Co-Authored-By: Andrew Johnson <drewej@protonmail.com>
Co-Authored-By: Patrick Shriwise <pshriwise@gmail.com>
Co-Authored-By: Paul Romano <paul.k.romano@gmail.com>

openmc/geometry.py

@@ -77,15 +77,22 @@ class Geometry(object):
                 if universe.id in volume_calc.volumes:
                     universe.add_volume_information(volume_calc)
 
-    def export_to_xml(self, path='geometry.xml'):
+    def export_to_xml(self, path='geometry.xml', remove_surfs=False):
         """Export geometry to an XML file.
 
         Parameters
         ----------
         path : str
             Path to file to write. Defaults to 'geometry.xml'.
+        remove_surfs : bool
+            Whether or not to remove redundant surfaces from the geometry when
+            exporting
 
         """
+        # Find and remove redundant surfaces from the geometry
+        if remove_surfs:
+            self.remove_redundant_surfaces()
+
         # Create XML representation
         root_element = ET.Element("geometry")
         self.root_universe.create_xml_subelement(root_element, memo=set())
@@ -382,6 +389,26 @@ class Geometry(object):
             surfaces = cell.region.get_surfaces(surfaces)
         return surfaces
 
+    def get_redundant_surfaces(self):
+        """Return all of the topologically redundant surface ids
+
+        Returns
+        -------
+        dict
+            Dictionary whose keys are the ID of a redundant surface and whose
+            values are the topologically equivalent :class:`openmc.Surface`
+            that should replace it.
+
+        """
+        tally = defaultdict(list)
+        for surf in self.get_all_surfaces().values():
+            coeffs = tuple(surf._coefficients[k] for k in surf._coeff_keys)
+            key = (surf._type,) + coeffs
+            tally[key].append(surf)
+        return {replace.id: keep 
+                for keep, *redundant in tally.values()
+                for replace in redundant}
+
     def get_materials_by_name(self, name, case_sensitive=False, matching=False):
         """Return a list of materials with matching names.
 
@@ -579,6 +606,17 @@ class Geometry(object):
 
         return sorted(lattices, key=lambda x: x.id)
 
+    def remove_redundant_surfaces(self):
+        """Remove redundant surfaces from the geometry"""
+
+        # Get redundant surfaces
+        redundant_surfaces = self.get_redundant_surfaces()
+
+        # Iterate through all cells contained in the geometry 
+        for cell in self.get_all_cells().values():
+            # Recursively remove redundant surfaces from regions
+            cell.region.remove_redundant_surfaces(redundant_surfaces)
+
     def determine_paths(self, instances_only=False):
         """Determine paths through CSG tree for cells and materials.
 

openmc/region.py

@@ -45,8 +45,7 @@ class Region(metaclass=ABCMeta):
         return not self == other
 
     def get_surfaces(self, surfaces=None):
-        """
-        Recursively find all the surfaces referenced by a region and return them
+        """Recursively find all surfaces referenced by a region and return them
 
         Parameters
         ----------
@@ -65,6 +64,19 @@ class Region(metaclass=ABCMeta):
             surfaces = region.get_surfaces(surfaces)
         return surfaces
 
+    def remove_redundant_surfaces(self, redundant_surfaces):
+        """Recursively remove all redundant surfaces referenced by this region
+
+        Parameters
+        ----------
+        redundant_surfaces : dict
+            Dictionary mapping redundant surface IDs to class:`openmc.Surface`
+            instances that should replace them.
+
+        """
+        for region in self:
+            region.remove_redundant_surfaces(redundant_surfaces)
+
     @staticmethod
     def from_expression(expression, surfaces):
         """Generate a region given an infix expression.
@@ -597,8 +609,7 @@ class Complement(Region):
         return temp_region.bounding_box
 
     def get_surfaces(self, surfaces=None):
-        """
-        Recursively find and return all the surfaces referenced by the node
+        """Recursively find and return all the surfaces referenced by the node
 
         Parameters
         ----------
@@ -617,6 +628,19 @@ class Complement(Region):
             surfaces = region.get_surfaces(surfaces)
         return surfaces
 
+    def remove_redundant_surfaces(self, redundant_surfaces):
+        """Recursively remove all redundant surfaces referenced by this region
+
+        Parameters
+        ----------
+        redundant_surfaces : dict
+            Dictionary mapping redundant surface IDs to class:`openmc.Surface`
+            instances that should replace them.
+
+        """
+        for region in self.node:
+            region.remove_redundant_surfaces(redundant_surfaces)
+
     def clone(self, memo=None):
         """Create a copy of this region - each of the surfaces in the
         complement's node will be cloned and will have new unique IDs.

openmc/surface.py

@@ -2209,6 +2209,21 @@ class Halfspace(Region):
         surfaces[self.surface.id] = self.surface
         return surfaces
 
+    def remove_redundant_surfaces(self, redundant_surfaces):
+        """Recursively remove all redundant surfaces referenced by this region
+
+        Parameters
+        ----------
+        redundant_surfaces : dict
+            Dictionary mapping redundant surface IDs to surface IDs for the 
+            :class:`openmc.Surface` instances that should replace them.
+
+        """
+
+        surf = redundant_surfaces.get(self.surface.id)
+        if surf is not None:
+            self.surface = surf
+
     def clone(self, memo=None):
         """Create a copy of this halfspace, with a cloned surface with a
         unique ID.

tests/unit_tests/test_geometry.py

@@ -301,3 +301,50 @@ def test_rotation_matrix():
     assert geom.find((0.0, 0.5, 0.0))[-1] == c3
     assert geom.find((0.0, -0.5, 0.0))[-1] == c1
     assert geom.find((0.0, -1.5, 0.0))[-1] == c2
+
+def test_remove_redundant_surfaces():
+    """Test ability to remove redundant surfaces"""
+
+    m1 = openmc.Material()
+    m1.add_nuclide('U235', 1.0, 'wo')
+    m1.add_nuclide('O16', 2.0, 'wo')
+    m1.set_density('g/cm3', 10.0)
+
+    m2 = openmc.Material()
+    m2.add_element('Zr', 1.0)
+    m2.set_density('g/cm3', 2.0)
+
+    m3 = openmc.Material()
+    m3.add_nuclide('H1', 2.0)
+    m3.add_nuclide('O16', 1.0)
+    m3.set_density('g/cm3', 1.0)
+
+    def get_cyl_cell(r1, r2, z1, z2, fill):
+        """Create a finite height cylindrical cell with a fill"""
+
+        cyl2 = openmc.ZCylinder(r=r2)
+        zplane1 = openmc.ZPlane(z1)
+        zplane2 = openmc.ZPlane(z2)
+        if np.isclose(r1, 0):
+            region = -cyl2 & +zplane1 & -zplane2
+        else:
+            cyl1 = openmc.ZCylinder(r=r1)
+            region = +cyl1 & -cyl2 & +zplane1 & -zplane2
+        return openmc.Cell(region=region, fill=fill)
+
+    r1, r2, r3 = 1., 2., 3.
+    z1, z2 = -2., 2.
+    fuel = get_cyl_cell(0, r1, z1, z2, m1)
+    clad = get_cyl_cell(r1, r2, z1, z2, m2)
+    water = get_cyl_cell(r2, r3, z1, z2, m3)
+    root = openmc.Universe(cells=[fuel, clad, water])
+    geom = openmc.Geometry(root)
+    
+    # There should be 6 redundant surfaces in this geometry
+    n_redundant_surfs = len(geom.get_redundant_surfaces().keys())
+    assert n_redundant_surfs == 6
+    # Remove redundant surfaces
+    geom.remove_redundant_surfaces()
+    # There should be 0 remaining redundant surfaces
+    n_redundant_surfs = len(geom.get_redundant_surfaces().keys())
+    assert n_redundant_surfs == 0