diff --git a/geometry.py b/geometry.py
index 0bfb70a..3dfcc5b 100755
--- a/geometry.py
+++ b/geometry.py
@@ -88,15 +88,33 @@ class Vertex(object):
     '''
     A vertex is a position along with other information such as color, normal vector and texture coordinates.
     '''
-    def __init__(self, x, y, z):
+    def __init__(self, x=0, y=0, z=0):
         self.x = x
         self.y = y
         self.z = z
         self.edges = []
-    
+        
     def __repr__(self):
         return "<%.2f, %.2f, %.2f>" % (self.x, self.y, self.z)
-
+        
+    def __add__(self, other):
+        # for now just assume type(other) = Vertex... bad, I know
+        return Vertex(self.x + other.x, self.y + other.y, self.z + other.z) 
+    
+    def __radd__(self, other):
+        return self.__add__(other)
+        
+    def __mul__(self, other):  
+        # for now just assume type(other) = int or float
+        return Vertex(self.x * other, self.y * other, self.z * other)
+    
+    def __rmul__(self, other):
+        return self.__mul__(other)
+    
+    def __div__(self, other):
+        # same assumption as __mult__
+        return Vertex(self.x / other, self.y / other, self.z / other)
+    
 class Edge(object):
     '''
     '''
@@ -104,34 +122,49 @@ class Edge(object):
         self.vertices = []
         self.faces = []
         self.edges = []
+        self.__edgeVertex = None
+        self.__subEdges = []
     
     def neighborFace(self, neighborFace):
         if neighborFace == self.faces[0]:
             return self.faces[1]
         else:
             return self.faces[0]
+    def __getMidPoint(self):
+        return sum(self.vertices, Vertex())/len(self.vertices)
+    midPoint = property(fget=__getMidPoint)
     
-    def edgeVertex(self):
+    def __getSubEdges(self):
+        if not self.__subEdges:
+            self.__subEdges = [Edge(), Edge()]
+            self.__subEdges[0].vertices = [self.vertices[0], self.edgeVertex]
+            self.__subEdges[1].vertices = [self.edgeVertex, self.vertices[1]]
+        return self.__subEdges
+    subEdges = property(fget=__getSubEdges)
+    
+    def __getEdgeVertex(self):
         '''
         Set each edge vertices to be the average of the two neighboring
             face vertices and its two original end vertices.
         '''
-        
-        # two neighboring face vertices:
-        neighborVertices = [face.centroid for face in self.faces]
-        
-        vertices = list(self.vertices)
-        vertices.extend(neighborVertices)
-        return self.__averageVertices(vertices)
+        if not self.__edgeVertex:
+            # two neighboring face vertices:
+            neighboringFaceVertices = [face.centroid for face in self.faces]
+            neighboringFaceVertices.extend(self.vertices)
+            xs = [vertex.x for vertex in neighboringFaceVertices]
+            ys = [vertex.y for vertex in neighboringFaceVertices]
+            zs = [vertex.z for vertex in neighboringFaceVertices]
+            x = sum(xs)/len(xs)
+            y = sum(ys)/len(ys)
+            z = sum(zs)/len(zs)
+            self.__edgeVertex = Vertex(x, y, z)
+            self.__edgeVertex.edges.extend(self.__subEdges)
+        return self.__edgeVertex
+    edgeVertex = property(fget=__getEdgeVertex)
     
     def __averageVertices(self, vertices):
-        xs = [vertex.x for vertex in vertices]
-        ys = [vertex.y for vertex in vertices]
-        zs = [vertex.z for vertex in vertices]
-        x = sum(xs)/len(xs)
-        y = sum(ys)/len(ys)
-        z = sum(zs)/len(zs)
-        return Vertex(x, y, z)
+        
+        return 
     
 class Face(object):
     '''
@@ -139,40 +172,105 @@ class Face(object):
     '''
     def __init__(self):
         self.edges = []
+        self.__centroid = None
+        self.__interiorEdges = []
+        self.__subFaces = []
         
     def __getCentroid(self):
-        # gather all vertex coords
-        faceVertices = self.__getFaceVertices()
-        xs = [vertex.x for vertex in faceVertices]
-        ys = [vertex.y for vertex in faceVertices]
-        zs = [vertex.z for vertex in faceVertices]
-        
-        # average each vertex component
-        x = sum(xs) / len(xs)
-        y = sum(ys) / len(ys)
-        z = sum(zs) / len(zs)
-        
-        return Vertex(x, y, z)
+        if not self.__centroid:
+            # gather all face vertex coords
+            faceVertices = list(set([vertex for edge in self.edges for vertex in edge.vertices]))
+            xs = [vertex.x for vertex in faceVertices]
+            ys = [vertex.y for vertex in faceVertices]
+            zs = [vertex.z for vertex in faceVertices]
+            
+            # average each vertex component
+            x = sum(xs) / len(xs)
+            y = sum(ys) / len(ys)
+            z = sum(zs) / len(zs)
+            
+            self.__centroid = Vertex(x, y, z)
+        return self.__centroid
     centroid = property(fget=__getCentroid)
     
-    def __getFaceVertices(self):
-        return list([vertex for edge in self.edges for vertex in edge.vertices])
-    vertices = property(fget=__getFaceVertices)
+    def __getSubFaces(self):
+        self.__setupSubDivisions()
+        return self.__subFaces
+    subFaces = property(fget=__getSubFaces)
     
+    def __getInteriorEdges(self):
+        self.__setupSubDivisions()
+        return self.__interiorEdges
+    interiorEdges = property(fget=__getInteriorEdges    )
+    
+    def __setupSubDivisions(self):
+        '''
+           v0      ev0    v1
+            *------e0-----*
+            |      |      |
+            |      |      |
+       ev3 e|11----f5----e|1 ev1
+            |      |      |
+            |      |      |
+            *------e2-----*
+           v3      ev2    v2
+        '''
+        if not self.__subFaces:
+            # create empty subFaces that will be filled with edge references below
+            #   these need to at least exist so the interior edges have something to reference
+            self.__subFaces = [Face() for edge in self.edges]
+            
+        if not self.__interiorEdges:
+            # set up empty edge objects to be filled below
+            self.__interiorEdges = [Edge() for edge in self.edges]
+            
+            # each interior edge connects the exterior edge vertex (mid-point) to the faceVertex (centroid)
+            for index in range(len(self.edges)):
+                prevIndex = (index - 1) % len(self.edges)
+                nextIndex = (index + 1) % len(self.edges)
+                
+                # end vertices are face centroid and currEdge edgeVertex
+                self.__interiorEdges[index].vertices = [self.edges[index].edgeVertex,
+                                                        self.centroid]
+                
+                # wing edges are the current edge's subEdges (ordered same as vertex order) and the prev and next interior edges
+                self.__interiorEdges[index].edges = [self.edges[index].subEdges[0],
+                                                     self.edges[index].subEdges[1],
+                                                     self.__interiorEdges[prevIndex],
+                                                     self.__interiorEdges[nextIndex]]
+                
+                # edge faces are the new subFaces (current and next faces), the current will be define below
+                #   and the next will be defined on the next iteration (or already defined on the last iteration)
+                self.__interiorEdges[index].faces = [self.__subFaces[index],
+                                                     self.__subFaces[nextIndex]]
+                
+                # now reference the current edge back into the faces,
+                #   and the edge.subEdges, and the edge.edgeVertex
+                
+                # current subFace (same index as current interior edge)
+                #   set its edges to reference the same edges used to setup the interior edge
+                #   order will be pretty important on these steps...
+                self.__subFaces[index].edges = [self.edges[index].subEdges[0],
+                                                self.__interiorEdges[index],
+                                                self.__interiorEdges[prevIndex],
+                                                self.edges[prevIndex].subEdges[1]]
+                
+                # just set one of the vertex edges, the other belongs to another face and will get added when that face is run
+                self.edges[index].edgeVertex.edges.append(self.__interiorEdges[index])
+                
+                self.edges[index].subEdges[0].faces.append(self.__subFaces[index])
+                self.edges[index].subEdges[0].faces.append(self.__subFaces[index])
+               
 class Polygon(object):
     '''
     Face splitting should happend on the polygon level(?).  It doesn't make sense to split just one face since
         it needs to average vertices with all adjoinging faces
     '''
     
-    def __init__(self, vertices, edges, faces):
-        self.vertices = vertices
-        self.edges = edges
-        self.faces = faces
-    
-    def __getFaceNeighbors(self, testFace):
-        for face in self.faces:
-            testFace.isNeighbor(face)
+    def __init__(self):
+        self.vertices = []
+        self.edges = []
+        self.faces = []
     
     def catmullClarkRefine(self):
         '''
@@ -184,126 +282,170 @@ class Polygon(object):
         For each original vertex P, take the average F of all n face vertices for faces touching P, and take the average R of all n edge midvertices for edges touching P, where each edge midvertex is the average of its two endvertex vertices. Move each original vertex to the vertex
         '''
         
-        # this calculates and returns the averaged vertex point
-        print self.faces[0].edges[0].edgeVertex()
         
-        # this returns a vertex at the face centroid
-        print self.faces[0].centroid
         
-        # this will get to the first neighborFace...
-        #   each face knows its faceVertex, and all of its edgeVertices
-        neighborFace = self.faces[0].edges[0].neighborFace(self.faces[0])
         
-        print neighborFace.centroid
+            
+            
+        # each face knows how to subdivide and create a set of subfaces, including interior edges and setup their references correctly... <- not completely finished...
+        p = Polygon()
+        edges = []
+        vertices = []
+        faces = []
+        
+        for face in self.faces:
+            for subFace in face.subFaces:
+                faces.append(subFace)
+                for edge in subFace.edges:
+                    edges.append(edge)
+                    for vertex in edge.vertices:
+                        vertices.append(vertex)
+        
+        for vertex in self.vertices:
+            faceVertices = []
+            edgeMidPoints = []
+            for edge in vertex.edges:
+                edgeMidPoints.append(edge.midPoint)
+                for face in edge.faces:
+                    faceVertices.append(face.centroid)
+            
+            f = sum(list(set(faceVertices)), Vertex())/len(list(set(faceVertices)))
+            
+            r = sum(list(set(edgeMidPoints)), Vertex())/len(list(set(edgeMidPoints)))
+            p = vertex
+            n = len(vertex.edges)
+            v = (f + (2.0 * r) + (n - 3.0) * p) / n
+            # print v
+            vertex.x = v.x
+            vertex.y = v.y
+            vertex.z = v.z
         
         # so now what.........
+        # (F + 2R + (n-3) P) / n
+        # 
+        # F = average of all face vertices touching P
+        # R = average of all edge vertices touching P
+        # P original point
+        # n = number of edges connecting to P
+
+        p.faces = faces
+        p.vertices = vertices
+        p.edges = edges
         
+        
+        # plotting these in excel seems to show the correct values (at first glace...)
+        
+        # so now what.........        
         # (F + 2R + (n-3) P) / n
         # 
         # F = average of all face vertices touching P
         # R = average of all edge vertices touching P
         # P original point
         # n = face vertices or edge vertices (should be the same number)
-      
+        return p
         
-v = []
-v.append(Vertex(0.0, 1.0, 0.0))
-v.append(Vertex(1.0, 1.0, 0.0))
-v.append(Vertex(1.0, 0.0, 0.0))
-v.append(Vertex(0.0, 0.0, 0.0))
-v.append(Vertex(0.0, 1.0, 1.0))
-v.append(Vertex(1.0, 1.0, 1.0))
-v.append(Vertex(1.0, 0.0, 1.0))
-v.append(Vertex(0.0, 0.0, 1.0))
+def createPolygon():
+    v = []
+    v.append(Vertex(0.0, 1.0, 0.0))
+    v.append(Vertex(1.0, 1.0, 0.0))
+    v.append(Vertex(1.0, 0.0, 0.0))
+    v.append(Vertex(0.0, 0.0, 0.0))
+    
+    v.append(Vertex(0.0, 1.0, 1.0))
+    v.append(Vertex(1.0, 1.0, 1.0))
+    v.append(Vertex(1.0, 0.0, 1.0))
+    v.append(Vertex(0.0, 0.0, 1.0))
 
-e = []
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
-e.append(Edge())
+    e = []
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
+    e.append(Edge())
 
-f = []
-f.append(Face())
-f.append(Face())
-f.append(Face())
-f.append(Face())
-f.append(Face())
-f.append(Face())
+    f = []
+    f.append(Face())
+    f.append(Face())
+    f.append(Face())
+    f.append(Face())
+    f.append(Face())
+    f.append(Face())
 
-# vertex list
-v[0].edges = [e[0], e[3], e[4]]
-v[1].edges = [e[0], e[5], e[1]]
-v[2].edges = [e[1], e[6], e[2]]
-v[3].edges = [e[2], e[7], e[3]]
-v[4].edges = [e[4], e[11], e[8]]
-v[5].edges = [e[5], e[9], e[8]]
-v[6].edges = [e[2], e[9], e[10]]
-v[7].edges = [e[7], e[10], e[11]]
+    # vertex list
+    v[0].edges = [e[0], e[3], e[4]]
+    v[1].edges = [e[0], e[5], e[1]]
+    v[2].edges = [e[1], e[6], e[2]]
+    v[3].edges = [e[2], e[7], e[3]]
+    v[4].edges = [e[4], e[11], e[8]]
+    v[5].edges = [e[5], e[9], e[8]]
+    v[6].edges = [e[6], e[9], e[10]]
+    v[7].edges = [e[7], e[10], e[11]]
 
-# face list
-f[0].edges = [e[0], e[1], e[2], e[3]]
-f[1].edges = [e[1], e[5], e[9], e[6]]
-f[2].edges = [e[2], e[6], e[10], e[7]]
-f[3].edges = [e[4], e[3], e[7], e[11]]
-f[4].edges = [e[8], e[5], e[0], e[4]]
-f[5].edges = [e[10], e[9], e[8], e[11]]
+    # face list
+    f[0].edges = [e[0], e[1], e[2], e[3]]
+    f[1].edges = [e[1], e[5], e[9], e[6]]
+    f[2].edges = [e[2], e[6], e[10], e[7]]
+    f[3].edges = [e[4], e[3], e[7], e[11]]
+    f[4].edges = [e[8], e[5], e[0], e[4]]
+    f[5].edges = [e[10], e[9], e[8], e[11]]
 
-#winged edges ordered by face, then by vertex reference
-e[0].vertices, e[0].faces, e[0].edges = [v[0], v[1]], [f[0], f[4]], [e[3], e[1], e[4], e[5]]
-e[1].vertices, e[1].faces, e[1].edges = [v[1], v[2]], [f[0], f[1]], [e[0], e[2], e[5], e[6]]
-e[2].vertices, e[2].faces, e[2].edges = [v[2], v[3]], [f[0], f[2]], [e[3], e[7], e[1], e[6]]
-e[3].vertices, e[3].faces, e[3].edges = [v[3], v[0]], [f[3], f[0]], [e[4], e[7], e[0], e[2]]
-e[4].vertices, e[4].faces, e[4].edges = [v[0], v[4]], [f[3], f[4]], [e[11], e[3], e[0], e[8]]
-e[5].vertices, e[5].faces, e[5].edges = [v[5], v[1]], [f[4], f[1]], [e[8], e[0], e[9], e[1]]
-e[6].vertices, e[6].faces, e[6].edges = [v[2], v[6]], [f[1], f[2]], [e[1], e[9], e[2], e[10]]
-e[7].vertices, e[7].faces, e[7].edges = [v[7], v[3]], [f[3], f[2]], [e[11], e[3], e[10], e[2]]
-e[8].vertices, e[8].faces, e[8].edges = [v[4], v[5]], [f[4], f[5]], [e[4], e[5], e[11], e[9]]
-e[9].vertices, e[9].faces, e[9].edges = [v[5], v[6]], [f[1], f[5]], [e[5], e[6], e[8], e[10]]
-e[10].vertices, e[10].faces, e[10].edges = [v[7], v[6]], [f[2], f[5]], [e[7], e[6], e[11], e[9]]
-e[11].vertices, e[11].faces, e[11].edges = [v[4], v[7]], [f[3], f[5]], [e[4], e[7], e[8], e[10]]
+    #winged edges ordered by face, then by vertex reference
+    e[0].vertices, e[0].faces, e[0].edges = [v[0], v[1]],    [f[0], f[4]], [e[3], e[1], e[4], e[5]]
+    e[1].vertices, e[1].faces, e[1].edges = [v[1], v[2]],    [f[0], f[1]], [e[0], e[2], e[5], e[6]]
+    e[2].vertices, e[2].faces, e[2].edges = [v[2], v[3]],    [f[0], f[2]], [e[1], e[3], e[6], e[7]]
+    e[3].vertices, e[3].faces, e[3].edges = [v[3], v[0]],    [f[3], f[0]], [e[4], e[7], e[0], e[2]]
+    e[4].vertices, e[4].faces, e[4].edges = [v[0], v[4]],    [f[3], f[4]], [e[11], e[3], e[0], e[8]]
+    e[5].vertices, e[5].faces, e[5].edges = [v[5], v[1]],    [f[4], f[1]], [e[8], e[0], e[9], e[1]]
+    e[6].vertices, e[6].faces, e[6].edges = [v[2], v[6]],    [f[1], f[2]], [e[1], e[9], e[2], e[10]]
+    e[7].vertices, e[7].faces, e[7].edges = [v[7], v[3]],    [f[3], f[2]], [e[11], e[3], e[10], e[2]]
+    e[8].vertices, e[8].faces, e[8].edges = [v[4], v[5]],    [f[4], f[5]], [e[4], e[5], e[11], e[9]]
+    e[9].vertices, e[9].faces, e[9].edges = [v[5], v[6]],    [f[1], f[5]], [e[5], e[6], e[8], e[10]]
+    e[10].vertices, e[10].faces, e[10].edges = [v[7], v[6]], [f[2], f[5]], [e[7], e[6], e[11], e[9]]
+    e[11].vertices, e[11].faces, e[11].edges = [v[4], v[7]], [f[3], f[5]], [e[4], e[7], e[8], e[10]]
 
-# just to prove to myself that the objects are the same, this is what years of pass by value have done to me...
+    # just to prove to myself that the objects are the same, this is what years of pass by value have done to me...
 
-# print id(v[0].x)
-# print id(e[0].vertices[0].x)
-# print id(f[0].edges[0].vertices[0].x)
+    # print id(v[0].x)
+    # print id(e[0].vertices[0].x)
+    # print id(f[0].edges[0].vertices[0].x)
 
-# v[0].x = 9
+    # v[0].x = 9
 
-# print id(v[0].x)
-# print id(e[0].vertices[0].x)
-# print id(f[0].edges[0].vertices[0].x)
+    # print id(v[0].x)
+    # print id(e[0].vertices[0].x)
+    # print id(f[0].edges[0].vertices[0].x)
 
-# print v[0].x
-# print e[0].vertices[0].x
-# print f[0].edges[0].vertices[0].x
+    # print v[0].x
+    # print e[0].vertices[0].x
+    # print f[0].edges[0].vertices[0].x
+    p = Polygon()
+    p.vertices = v
+    p.edges = e
+    p.faces = f
+    return p
 
-polygon = Polygon(v, e, f)
-polygon.catmullClarkRefine()
+polygon = createPolygon()
+newPolygon = polygon.catmullClarkRefine()
 
+from numpy import *
+import pylab
+import mpl_toolkits.mplot3d.axes3d as p3
 
-
-# from numpy import *
-# import pylab
-# import mpl_toolkits.mplot3d.axes3d as p3
-
-# fig = pylab.figure()
-# ax = p3.Axes3D(fig)
-# ax.plot_wireframe(xs, ys, zs)
-# ax.set_xlabel('X')
-# ax.set_ylabel('Y')
-# ax.set_zlabel('Z')
-# pylab.show()
-
-# v4 = Vertex(4, 0, 0, 1)
-# v5 = Vertex(5, 1, 0, 1)
-# v6 = Vertex(6, 0, 1, 1)
-# v7 = Vertex(7, 1, 1, 1)
+fig = pylab.figure()
+ax = p3.Axes3D(fig)
+for edge in newPolygon.edges:
+    xs = [vertex.x for vertex in edge.vertices]
+    ys = [vertex.y for vertex in edge.vertices]
+    zs = [vertex.z for vertex in edge.vertices]
+    ax.plot_wireframe(xs, ys, zs)
+ax.set_xlabel('X')
+ax.set_ylabel('Y')
+ax.set_zlabel('Z')
+pylab.show()
\ No newline at end of file