diff --git a/roundedcube.scad b/roundedcube.scad
new file mode 100644
index 0000000..da0baea
--- /dev/null
+++ b/roundedcube.scad
@@ -0,0 +1,148 @@
+use<functional.scad>
+
+/* Rounded cube using the cube/octahedron duality. */
+
+function rounded_cube(size, r=0, center=false, fn=12)=
+   (!is_array(size))?
+     rounded_cube(size=[size,size,size],r=r,center=center,fn=fn)
+   :(r>=min(size[0],size[1],size[2])/2)?
+      rounded_cube(size=size,r=min(size[0],size[1],size[2])/2-.0001,center=center,fn=fn)
+   :(r<=0)?
+      cube(size,center=center)
+   :
+      let(
+         steps=(fn>0)?ceil(fn/4):1,
+         sx=size[0]/2-r,
+         sy=size[1]/2-r,
+         sz=size[2]/2-r,
+         da=90/steps,
+         count=4*(steps+1)*(steps+2)-1,
+         rcube=[
+              flatten(concat(
+                  [for (i=[0:steps])
+                      let(
+                         z=r*(-cos(da*i))-sz,
+                         l=r*sin(da*i),
+                         dp=(i>0)?(90/i):0
+                      )
+                      concat(
+                         [for(j=[0:i]) [sx+l*cos(dp*j),sy+l*sin(dp*j),z]],
+                         [for(j=[0:i]) [-sx-l*sin(dp*j),sy+l*cos(dp*j),z]],
+                         [for(j=[0:i]) [-sx-l*cos(dp*j),-sy-l*sin(dp*j),z]],
+                         [for(j=[0:i]) [sx+l*sin(dp*j),-sy-l*cos(dp*j),z]]
+                      )
+                  ],
+                  [for (i=[steps:-1:0])
+                      let(
+                         z=sz+r*cos(da*i),
+                         l=r*sin(da*i),
+                         dp=(i>0)?90/i:0
+                      )
+                      concat(
+                         [for(j=[0:i]) [sx+l*cos(dp*j),sy+l*sin(dp*j),z]],
+                         [for(j=[0:i]) [-sx-l*sin(dp*j),sy+l*cos(dp*j),z]],
+                         [for(j=[0:i]) [-sx-l*cos(dp*j),-sy-l*sin(dp*j),z]],
+                         [for(j=[0:i]) [sx+l*sin(dp*j),-sy-l*cos(dp*j),z]]
+                      )
+                  ]
+            )),
+//Each side on the i'th row has i+1 points
+//The rows start with 0,4,8,12 ... 2*i*(i+1)
+//Each side on the ith row has i+1 points
+                 concat(
+                    [[0,1,2,3]],
+                    [for(i=[1:steps]) let(start=2*i*(i+1),laststart=2*i*(i-1),
+                                          ppl=4*(i+1),lppl=4*i,
+                                          pps=(i+1),lpps=i
+                                          )
+                       for(side=[0:3])
+                       [
+                          laststart+lpps*side,
+                          laststart+(lpps*side+lppl-1)%lppl,
+                          start+(pps*side+ppl-1)%ppl,
+                          start+pps*side
+                       ]
+                    ],
+                    [for(i=[1:steps]) let(start=2*i*(i+1),laststart=2*i*(i-1),
+                                          ppl=4*(i+1),lppl=4*i,
+                                          pps=(i+1),lpps=i
+                                          )
+                       for(side=[0:3]) for(j=[0:i-1])
+                       [
+                          laststart+lpps*side+j,
+                          start+pps*side+j,
+                          start+pps*side+1+j
+                       ]
+                    ],
+                    (steps<2)?[]: [for(i=[2:steps]) let(start=2*i*(i+1),laststart=2*i*(i-1),
+                                          ppl=4*(i+1),lppl=4*i,
+                                          pps=(i+1),lpps=i
+                                          )
+                       for(side=[0:3]) for(j=[0:i-2])
+                       [
+                          laststart+lpps*side+j,
+                          start+pps*side+j+1,
+                          laststart+lpps*side+j+1
+                       ]
+                    ],
+//These faces fuse the top & bottom
+//The steps'th row starts with 2*steps*(steps+1)
+//The (steps+1)'th row starts with 2*steps*(steps+1)+4*steps
+                    [for(side=[0:3]) let(start=2*steps*(steps+1)+4*(steps+1),
+                                          laststart=2*steps*(steps+1),
+                                          ppl=4*(steps+1),lppl=4*(steps+1),
+                                          pps=(steps+1),lpps=(steps+1)
+                                          ) for(j=[0:steps])
+                       
+                       [
+                          laststart+lpps*side+j,
+                          laststart+(lpps*side+lppl-1+j)%lppl,
+                          start+(pps*side+ppl-1+j)%ppl,
+                          start+pps*side+j
+                       ]
+                    ],
+//And the top is just the bottom wound in reverse.
+                    [[count-0,count-1,count-2,count-3]],
+                    [for(i=[1:steps]) let(start=2*i*(i+1),laststart=2*i*(i-1),
+                                          ppl=4*(i+1),lppl=4*i,
+                                          pps=(i+1),lpps=i
+                                          )
+                       for(side=[0:3])
+                       [
+                          count-(laststart+lpps*side),
+                          count-(laststart+(lpps*side+lppl-1)%lppl),
+                          count-(start+(pps*side+ppl-1)%ppl),
+                          count-(start+pps*side)
+                       ]
+                    ],
+                    [for(i=[1:steps]) let(start=2*i*(i+1),laststart=2*i*(i-1),
+                                          ppl=4*(i+1),lppl=4*i,
+                                          pps=(i+1),lpps=i
+                                          )
+                       for(side=[0:3]) for(j=[0:i-1])
+                       [
+                          count-(laststart+lpps*side+j),
+                          count-(start+pps*side+j),
+                          count-(start+pps*side+1+j)
+                       ]
+                    ],
+                    (steps<2)?[]: [for(i=[2:steps]) let(start=2*i*(i+1),laststart=2*i*(i-1),
+                                          ppl=4*(i+1),lppl=4*i,
+                                          pps=(i+1),lpps=i
+                                          )
+                       for(side=[0:3]) for(j=[0:i-2])
+                       [
+                          count-(laststart+lpps*side+j),
+                          count-(start+pps*side+j+1),
+                          count-(laststart+lpps*side+j+1)
+                       ]
+                    ]
+                 )
+          ]
+   )
+   center?rcube:translate(size/2,rcube)
+;
+
+rcube=(rounded_cube([3,2,1],1/3,fn=24,center=false));
+poly3d(rcube);
+
diff --git a/subdivision.scad b/subdivision.scad
index 27f9f43..e417f46 100644
--- a/subdivision.scad
+++ b/subdivision.scad
@@ -1,6 +1,7 @@
 use <functional.scad>
 
 // faces better be triangles... or else!
+
 function subdivide_faces(n=1, poly) = 
   let(
     points = poly[0],
@@ -28,3 +29,122 @@ function subdivide_faces(n=1, poly) =
     ])
   )
   n > 1 ? subdivide_faces(n-1, [allpoints, faces]) : n > 0 ? [allpoints,faces] : poly;
+
+//Utility function for splitting circular vectors.
+function circle_cut (list,start,end)=
+   (end<start)?
+      concat([for (i=[start:len(list)-1]) list[i]],[for (i=[0:end]) list[i]])
+   :
+      [for (i=[start:end]) list[i]]
+;
+
+//Ineffictient (Order N^2) naive face triangulation function that works
+//by ear clipping:
+//Finds a convex vertex (order N)
+//Checks if the 'ear' with that vertex contains any points.
+//If not, the ear can be clipped, otherwise there's a diagonal in the ear,
+
+function triangulate_poly3d(poly)=
+   [poly[0],triangulate_faces(poly[0],poly[1])]
+;
+
+function triangulate_faces(points,faces)=
+   flatten([for(i=[0:len(faces)-1]) triangulate_face(points,faces[i])])
+;
+
+function triangulate_face(points,face)=
+   let(count=len(face))
+   (3==count)?
+      [face]
+   :
+      let(wd=face_winding(points,face),cv=convex_vertex(wd,points,face),
+          pv=(count+cv-1)%count,nv=(cv+1)%count,
+          p0=points[face[pv]],p1=points[face[cv]],p2=points[face[nv]],
+          tests=[
+             [cross(wd,p0-p2),cross(wd,p0-p2)*p0],
+             [cross(wd,p1-p0),cross(wd,p1-p0)*p1],
+             [cross(wd,p2-p1),cross(wd,p2-p1)*p2]
+          ],
+          eartest=point_in_ear(points,face,tests),
+          clipableear=(eartest[0]<0),
+          diagonalpoint=eartest[1]
+       )
+       (clipableear)? //There is no point inside the ear.
+          flatten([
+             [circle_cut(face,pv,nv)],
+             triangulate_face(points,circle_cut(face,nv,pv))
+          ])
+       : //If there is a point inside the ear, make a diagonal and clip along that.
+          flatten([
+             triangulate_face(points,circle_cut(face,cv,diagonalpoint)),
+             triangulate_face(points,circle_cut(face,diagonalpoint,cv))
+          ])
+;
+
+function vsum(v_list,i=0)=
+  (i<len(v_list)-1)?
+     v_list[i]+vsum(v_list,i+1)
+  :
+     [0,0,0]
+;
+
+function face_winding(points,face)=
+   let(count=len(face))
+   vsum([
+      for(i=[0:count-3]) cross(
+         points[face[(i+1)]]-points[face[0]],
+         points[face[(i+2)]]-points[face[(i+1)]]
+      )
+   ])
+;
+
+function convex_vertex(wd,points,face,i=0)=
+   let(count=len(face),
+       p0=points[face[i]],p1=points[face[(i+1)%count]],p2=points[face[(i+2)%count]]
+   )
+   (len(face)>i)?
+      (cross(p1-p0,p2-p1)*wd>0)? 
+         (i+1)%count
+      :
+         convex_vertex(wd,points,face,i+1)
+   ://This should never happen since there is at least 1 convex vertex.
+      undef
+;
+
+//The order of tests matters - This assumes Test 0 indicates which point is furthest in 
+//the ear - and that point will form a diagonal with the tip of the ear.
+
+function point_in_ear(points,face,tests,i=0)=
+    (i<len(face)-1)?
+       let(
+           prev=point_in_ear(points,face,tests,i+1),
+           test=check_point_in_ear(points[face[i]],tests)
+       )
+       (test>prev[0])?
+          [test,i]
+       :
+          prev
+    :
+       [check_point_in_ear(points[face[i]],tests),i]
+;
+
+function check_point_in_ear (point,tests)=
+   let(
+      result=[
+         (point*tests[0][0])-tests[0][1],
+         (point*tests[1][0])-tests[1][1],
+         (point*tests[2][0])-tests[2][1]
+      ]
+   )
+   (result[0]>0 && result[1]>0 && result[2]>0)?
+      result[0]
+   :
+      -1
+;
+
+
+
+
+
+
+