From 32eb58774df061c2198814ab99be6db7058288f6 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 10 May 2016 22:08:10 -0500
Subject: [PATCH] Add make_hexagon_region() function

---
 docs/source/pythonapi/index.rst | 10 ++++++++
 openmc/surface.py               | 45 ++++++++++++++++++++++++++++++++-
 2 files changed, 54 insertions(+), 1 deletion(-)

diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst
index 1631976e6..89d1b0508 100644
--- a/docs/source/pythonapi/index.rst
+++ b/docs/source/pythonapi/index.rst
@@ -122,6 +122,16 @@ Many of the above classes are derived from several abstract classes:
    openmc.Region
    openmc.Lattice
 
+One function is also available to create a hexagonal region defined by the
+intersection of six surface half-spaces.
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myfunction.rst
+
+   openmc.make_hexagon_region
+
 Constructing Tallies
 --------------------
 
diff --git a/openmc/surface.py b/openmc/surface.py
index 37e7c2ffd..84028c1af 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -2,11 +2,12 @@ from abc import ABCMeta
 from numbers import Real, Integral
 from xml.etree import ElementTree as ET
 import sys
+from math import sqrt
 
 import numpy as np
 
 from openmc.checkvalue import check_type, check_value, check_greater_than
-from openmc.region import Region
+from openmc.region import Region, Intersection
 
 if sys.version_info[0] >= 3:
     basestring = str
@@ -1503,3 +1504,45 @@ class Halfspace(Region):
     def __str__(self):
         return '-' + str(self.surface.id) if self.side == '-' \
             else str(self.surface.id)
+
+
+def make_hexagon_region(edge_length=1., orientation='y'):
+    """Create a hexagon region from six surface planes.
+
+    Parameters
+    ----------
+    edge_length : float
+        Length of a side of the hexagon in cm
+    orientation : {'x', 'y'}
+        An 'x' orientation means that two sides of the hexagon are parallel to
+        the x-axis and a 'y' orientation means that two sides of the hexagon are
+        parallel to the y-axis.
+
+    Returns
+    -------
+    openmc.Region
+        The inside of a hexagonal prism
+
+    """
+
+    l = edge_length
+
+    if orientation == 'x':
+        right = XPlane(x0=sqrt(3.)/2.*l)
+        left = XPlane(x0=-sqrt(3.)/2.*l)
+        c = sqrt(3.)/3.
+        ur = Plane(A=c, B=1., D=l)  # y = -x/sqrt(3) + a
+        ul = Plane(A=-c, B=1., D=l)  # y = x/sqrt(3) + a
+        lr = Plane(A=-c, B=1., D=-l) # y = x/sqrt(3) - a
+        ll = Plane(A=c, B=1., D=-l)  # y = -x/sqrt(3) - a
+        return Intersection(-right, +left, -ur, -ul, +lr, +ll)
+
+    elif orientation == 'y':
+        top = YPlane(y0=sqrt(3.)/2.*l)
+        bottom = YPlane(y0=-sqrt(3.)/2.*l)
+        c = sqrt(3.)
+        ur = Plane(A=c, B=1., D=c*l)  # y = -sqrt(3)*(x - a)
+        lr = Plane(A=-c, B=1., D=-c*l)  # y = sqrt(3)*(x + a)
+        ll = Plane(A=c, B=1., D=-c*l) # y = -sqrt(3)*(x + a)
+        ul = Plane(A=-c, B=1., D=c*l)  # y = sqrt(3)*(x + a)
+        return Intersection(-top, +bottom, -ur, +lr, +ll, -ul)