diff --git a/Scratch.ipynb b/Scratch.ipynb
index 007a7a0..bb1db6a 100644
--- a/Scratch.ipynb
+++ b/Scratch.ipynb
@@ -14,7 +14,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -23,11 +23,34 @@
     "    [1, 0, 0],\n",
     "    [1, 1, 0],\n",
     "    [0, 1, 0],\n",
-    "], dtype=numpy.float64)\n",
-    "b2 = numpy.array(b1 + [0,0,1])\n",
+    "], dtype=numpy.float64) - 0.5\n",
+    "b2 = numpy.array(b1 + [0,0,1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ang = 0.15\n",
+    "xform = meshutil.Transform().translate(0,0,0.5)\n",
     "\n",
-    "testmesh = meshutil.join_boundary_simple(b1, b2)\n",
-    "testmesh.to_stl_mesh().save(\"simple.stl\")"
+    "mesh = meshutil.join_boundary_simple(b1, b2)\n",
+    "for i,bound in enumerate(meshutil.subdivide_boundary(b2)):\n",
+    "    m = xform.apply_to(bound) #xform.rotate([0,0,1], i*numpy.pi/2).apply_to(bound)\n",
+    "    mesh = mesh.concat(meshutil.join_boundary_simple(bound, m))\n",
+    "\n",
+    "# One problem: without an actual transform in subdivide_boundary,\n",
+    "# the transforms aren't inherited, and there is no transformed space.\n",
+    "# I'm just dealing directly with transformed geometry.\n",
+    "#\n",
+    "# My rotation around z in the loop above won't work the way I intended.\n",
+    "#\n",
+    "# In other words: I don't need to subdivide *geometry*.\n",
+    "# I need to subdivide *space* and then put geometry in it.\n",
+    "\n",
+    "mesh.to_stl_mesh().save(\"simple.stl\")"
    ]
   },
   {
diff --git a/meshutil.py b/meshutil.py
index 24779fe..597fa08 100644
--- a/meshutil.py
+++ b/meshutil.py
@@ -55,6 +55,9 @@ class Transform(object):
         return self._compose(mtx_translate(*a, **kw))
     def rotate(self, *a, **kw):
         return self._compose(mtx_rotate(*a, **kw))
+    def apply_to(self, vs):
+        vh = numpy.hstack([vs, numpy.ones((vs.shape[0], 1), dtype=vs.dtype)])
+        return vh.dot(self.mtx.T)[:,0:3]
 
 def mtx_scale(sx, sy=None, sz=None):
     if sy is None:
@@ -136,6 +139,21 @@ def cube_distort(angle, open_xz=False):
         # winding order?
     return FaceVertexMesh(verts, faces)
 
+def subdivide_boundary(bound):
+    # assume bound1 has shape (4,3).
+    # Midpoints of every segment:
+    mids = (bound + numpy.roll(bound, 1, axis=0)) / 2
+    mids_adj = numpy.roll(mids, -1, axis=0)
+    # Centroid:
+    centroid = numpy.mean(bound, axis=0)
+    # Now, every single new boundary has: one vertex of 'bound', an
+    # adjacent midpoint, a centroid, and the other adjacent midpoint.
+    bounds = [
+        numpy.array([bound[i,:], mids[i,:], centroid, mids_adj[i,:]])
+        for i in range(4)
+    ]
+    return bounds
+
 def join_boundary_simple(bound1, bound2):
     # bound1 & bound2 are both arrays of shape (N,3), representing
     # the points of a boundary.  This joins the two boundaries by