diff --git a/Makefile b/Makefile
index 1d906d9..3df3ed4 100644
--- a/Makefile
+++ b/Makefile
@@ -40,6 +40,7 @@ SSOBJS   = $(OBJDIR)\solvespace.obj \
            $(OBJDIR)\export.obj \
 
 SRFOBJS =  $(OBJDIR)\ratpoly.obj \
+           $(OBJDIR)\triangulate.obj \
 
 
 RES = $(OBJDIR)\resource.res
@@ -55,7 +56,7 @@ all: $(OBJDIR)/solvespace.exe
 clean:
 	rm -f obj/*
 
-$(OBJDIR)/solvespace.exe: $(SSOBJS) $(SRFOBJS) $(W32OBJS) $(FREEZE) $(RES)
+$(OBJDIR)/solvespace.exe: $(SRFOBJS) $(SSOBJS) $(W32OBJS) $(FREEZE) $(RES)
     @$(CC) $(DEFINES) $(CFLAGS) -Fe$(OBJDIR)/solvespace.exe $(SSOBJS) $(SRFOBJS) $(W32OBJS) $(FREEZE) $(RES) $(LIBS)
     editbin /nologo /STACK:8388608 $(OBJDIR)/solvespace.exe
     @echo solvespace.exe
diff --git a/dsc.h b/dsc.h
index 54b9229..f3f9692 100644
--- a/dsc.h
+++ b/dsc.h
@@ -105,6 +105,15 @@ public:
         elem[n++] = *t;
     }
 
+    T *First(void) {
+        return (n == 0) ? NULL : &(elem[0]);
+    }
+    T *NextAfter(T *prev) {
+        if(!prev) return NULL;
+        if(prev - elem == (n - 1)) return NULL;
+        return prev + 1;
+    }
+
     void ClearTags(void) {
         int i;
         for(i = 0; i < n; i++) {
@@ -218,6 +227,15 @@ public:
         return NULL;
     }
 
+    T *First(void) {
+        return (n == 0) ? NULL : &(elem[0]);
+    }
+    T *NextAfter(T *prev) {
+        if(!prev) return NULL;
+        if(prev - elem == (n - 1)) return NULL;
+        return prev + 1;
+    }
+
     void ClearTags(void) {
         int i;
         for(i = 0; i < n; i++) {
diff --git a/glhelper.cpp b/glhelper.cpp
index 9c381b0..675b114 100644
--- a/glhelper.cpp
+++ b/glhelper.cpp
@@ -166,8 +166,6 @@ void glxFillMesh(int specColor, SMesh *m, DWORD h, DWORD s1, DWORD s2)
     glBegin(GL_TRIANGLES);
     for(int i = 0; i < m->l.n; i++) {
         STriangle *tr = &(m->l.elem[i]);
-        Vector n = tr->Normal();
-        glNormal3d(n.x, n.y, n.z);
 
         int color;
         if(specColor < 0) {
@@ -183,9 +181,24 @@ void glxFillMesh(int specColor, SMesh *m, DWORD h, DWORD s1, DWORD s2)
             glBegin(GL_TRIANGLES);
         }
 
-        glxVertex3v(tr->a);
-        glxVertex3v(tr->b);
-        glxVertex3v(tr->c);
+        if(tr->an.EqualsExactly(Vector::From(0, 0, 0))) {
+            // Compute the normal from the vertices
+            Vector n = tr->Normal();
+            glNormal3d(n.x, n.y, n.z);
+            glxVertex3v(tr->a);
+            glxVertex3v(tr->b);
+            glxVertex3v(tr->c);
+        } else {
+            // Use the exact normals that are specified
+            glNormal3d((tr->an).x, (tr->an).y, (tr->an).z);
+            glxVertex3v(tr->a);
+
+            glNormal3d((tr->bn).x, (tr->bn).y, (tr->bn).z);
+            glxVertex3v(tr->b);
+
+            glNormal3d((tr->cn).x, (tr->cn).y, (tr->cn).z);
+            glxVertex3v(tr->c);
+        }
 
         if((s1 != 0 && tr->meta.face == s1) || 
            (s2 != 0 && tr->meta.face == s2))
diff --git a/graphicswin.cpp b/graphicswin.cpp
index 32d546e..f3a1391 100644
--- a/graphicswin.cpp
+++ b/graphicswin.cpp
@@ -105,6 +105,7 @@ const GraphicsWindow::MenuEntry GraphicsWindow::menu[] = {
 
 { 0, "&Analyze",                            0,                          NULL  },
 { 1, "Measure &Volume\tCtrl+Shift+V",       MNU_VOLUME,         'V'|S|C,mAna  },
+{ 1, "Show &Naked Edges\tCtrl+Shift+N",     MNU_NAKED_EDGES,    'N'|S|C,mAna  },
 { 1, NULL,                                  0,                          NULL  },
 { 1, "Show Degrees of &Freedom\tCtrl+Shift+F", MNU_SHOW_DOF,    'F'|S|C,mAna  },
 { 1, NULL,                                  0,                          NULL  },
diff --git a/groupmesh.cpp b/groupmesh.cpp
index 1ec4361..12698b3 100644
--- a/groupmesh.cpp
+++ b/groupmesh.cpp
@@ -140,6 +140,7 @@ void Group::GenerateMeshForStepAndRepeat(void) {
 
 void Group::GenerateMesh(void) {
     thisMesh.Clear();
+    thisShell.Clear();
     STriMeta meta = { 0, color };
 
     if(type == TRANSLATE || type == ROTATE) {
@@ -160,6 +161,10 @@ void Group::GenerateMesh(void) {
         } else {
             tbot = translate.ScaledBy(-1); ttop = translate.ScaledBy(1);
         }
+        
+        thisShell = SShell::FromExtrusionOf(&(src->bezierLoopSet), tbot, ttop);
+        thisShell.TriangulateInto(&thisMesh);
+/*
         bool flipBottom = translate.Dot(src->poly.normal) > 0;
 
         // Get a triangulation of the source poly; this is not a closed mesh.
@@ -221,7 +226,7 @@ void Group::GenerateMesh(void) {
                 thisMesh.AddTriangle(meta, bbot, btop, atop);
             }
         }
-        edges.Clear();
+        edges.Clear(); */
     } else if(type == LATHE) {
         SEdgeList edges;
         ZERO(&edges);
@@ -291,6 +296,7 @@ void Group::GenerateMesh(void) {
     }
 
     runningMesh.Clear();
+    runningShell.Clear();
 
     // If this group contributes no new mesh, then our running mesh is the
     // same as last time, no combining required. Likewise if we have a mesh
@@ -308,6 +314,8 @@ void Group::GenerateMesh(void) {
     SMesh *a = PreviousGroupMesh();
     if(meshCombine == COMBINE_AS_UNION) {
         runningMesh.MakeFromUnion(a, &thisMesh);
+        runningMesh = thisMesh;
+        ZERO(&thisMesh);
     } else if(meshCombine == COMBINE_AS_DIFFERENCE) {
         runningMesh.MakeFromDifference(a, &thisMesh);
     } else {
diff --git a/polygon.cpp b/polygon.cpp
index b7cbd6e..ab60e6b 100644
--- a/polygon.cpp
+++ b/polygon.cpp
@@ -39,10 +39,16 @@ bool STriangle::ContainsPointProjd(Vector n, Vector p) {
 
 void STriangle::FlipNormal(void) {
     SWAP(Vector, a, b);
+    SWAP(Vector, an, bn);
 }
 
 STriangle STriangle::From(STriMeta meta, Vector a, Vector b, Vector c) {
-    STriangle tr = { 0, meta, a, b, c };
+    STriangle tr;
+    ZERO(&tr);
+    tr.meta = meta;
+    tr.a = a;
+    tr.b = b;
+    tr.c = c;
     return tr;
 }
 
diff --git a/polygon.h b/polygon.h
index 71b201a..742ab86 100644
--- a/polygon.h
+++ b/polygon.h
@@ -77,6 +77,7 @@ public:
     int         tag;
     STriMeta    meta;
     Vector      a, b, c;
+    Vector      an, bn, cn;
 
     static STriangle From(STriMeta meta, Vector a, Vector b, Vector c);
     Vector Normal(void);
diff --git a/sketch.h b/sketch.h
index 226e9dc..5af35bc 100644
--- a/sketch.h
+++ b/sketch.h
@@ -154,6 +154,8 @@ public:
         bool            yes;
     }               meshError;
     SEdgeList       emphEdges;
+    SShell          thisShell;
+    SShell          runningShell;
 
     static const int COMBINE_AS_UNION           = 0;
     static const int COMBINE_AS_DIFFERENCE      = 1;
diff --git a/solvespace.cpp b/solvespace.cpp
index 8a8c11e..9c3f5c8 100644
--- a/solvespace.cpp
+++ b/solvespace.cpp
@@ -399,6 +399,10 @@ void SolveSpace::MenuAnalyze(int id) {
             }
             break;
 
+        case GraphicsWindow::MNU_NAKED_EDGES: {
+            break;
+        }
+
         case GraphicsWindow::MNU_VOLUME: {
             SMesh *m = &(SS.GetGroup(SS.GW.activeGroup)->runningMesh);
            
diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp
index 27078d9..1ac80b9 100644
--- a/srf/ratpoly.cpp
+++ b/srf/ratpoly.cpp
@@ -2,7 +2,13 @@
 
 double Bernstein(int k, int deg, double t)
 {
+    if(k > deg || k < 0) return 0;
+
     switch(deg) {
+        case 0:
+            return 1;
+            break;
+
         case 1:
             if(k == 0) {
                 return (1 - t);
@@ -36,6 +42,11 @@ double Bernstein(int k, int deg, double t)
     oops();
 }
 
+double BernsteinDerivative(int k, int deg, double t)
+{
+    return deg*(Bernstein(k-1, deg-1, t) - Bernstein(k, deg-1, t));
+}
+
 SBezier SBezier::From(Vector p0, Vector p1) {
     SBezier ret;
     ZERO(&ret);
@@ -92,11 +103,15 @@ Vector SBezier::PointAt(double t) {
 }
 
 void SBezier::MakePwlInto(List<Vector> *l) {
-    l->Add(&(ctrl[0]));
-    MakePwlWorker(l, 0.0, 1.0);
+    MakePwlInto(l, Vector::From(0, 0, 0));
 }
+void SBezier::MakePwlInto(List<Vector> *l, Vector offset) {
+    Vector p = (ctrl[0]).Plus(offset);
+    l->Add(&p);
 
-void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb) {
+    MakePwlWorker(l, 0.0, 1.0, offset);
+}
+void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb, Vector off) {
     Vector pa = PointAt(ta);
     Vector pb = PointAt(tb);
 
@@ -115,11 +130,12 @@ void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb) {
     double step = 1.0/SS.maxSegments;
     if((tb - ta) < step || d < tol) {
         // A previous call has already added the beginning of our interval.
+        pb = pb.Plus(off);
         l->Add(&pb);
     } else {
         double tm = (ta + tb) / 2;
-        MakePwlWorker(l, ta, tm);
-        MakePwlWorker(l, tm, tb);
+        MakePwlWorker(l, ta, tm, off);
+        MakePwlWorker(l, tm, tb, off);
     }
 }
 
@@ -131,10 +147,12 @@ void SBezier::Reverse(void) {
     }
 }
 
+
 void SBezierList::Clear(void) {
     l.Clear();
 }
 
+
 SBezierLoop SBezierLoop::FromCurves(SBezierList *sbl,
                                     bool *allClosed, SEdge *errorAt)
 {
@@ -194,6 +212,12 @@ SBezierLoop SBezierLoop::FromCurves(SBezierList *sbl,
 
 void SBezierLoop::Reverse(void) {
     l.Reverse();
+    SBezier *sb;
+    for(sb = l.First(); sb; sb = l.NextAfter(sb)) {
+        // If we didn't reverse each curve, then the next curve in list would
+        // share your start, not your finish.
+        sb->Reverse();
+    }
 }
 
 void SBezierLoop::MakePwlInto(SContour *sc) {
@@ -216,6 +240,7 @@ void SBezierLoop::MakePwlInto(SContour *sc) {
     sc->l.elem[sc->l.n - 1] = sc->l.elem[0];
 }
 
+
 SBezierLoopSet SBezierLoopSet::From(SBezierList *sbl, SPolygon *poly,
                                     bool *allClosed, SEdge *errorAt)
 {
@@ -240,6 +265,11 @@ SBezierLoopSet SBezierLoopSet::From(SBezierList *sbl, SPolygon *poly,
 
     poly->normal = poly->ComputeNormal();
     ret.normal = poly->normal;
+    if(poly->l.n > 0) {
+        ret.point = poly->AnyPoint();
+    } else {
+        ret.point = Vector::From(0, 0, 0);
+    }
     poly->FixContourDirections();
 
     for(i = 0; i < poly->l.n; i++) {
@@ -262,6 +292,21 @@ void SBezierLoopSet::Clear(void) {
     l.Clear();
 }
 
+void SCurve::Clear(void) {
+    pts.Clear();
+}
+
+STrimBy STrimBy::EntireCurve(SShell *shell, hSCurve hsc) {
+    STrimBy stb;
+    ZERO(&stb);
+    stb.curve = hsc;
+    SCurve *sc = shell->curve.FindById(hsc);
+    stb.start = sc->pts.elem[0];
+    stb.finish = sc->pts.elem[sc->pts.n - 1];
+
+    return stb;
+}
+
 SSurface SSurface::FromExtrusionOf(SBezier *sb, Vector t0, Vector t1) {
     SSurface ret;
     ZERO(&ret);
@@ -281,25 +326,293 @@ SSurface SSurface::FromExtrusionOf(SBezier *sb, Vector t0, Vector t1) {
     return ret;
 }
 
-SShell SShell::FromExtrusionOf(SBezierList *sbl, Vector t0, Vector t1) {
-    SShell ret;
+SSurface SSurface::FromPlane(Vector pt, Vector n) {
+    SSurface ret;
     ZERO(&ret);
 
-    // Group the input curves into loops, not necessarily in the right order.
+    ret.degm = 1;
+    ret.degn = 1;
 
+    Vector u = n.Normal(0), v = n.Normal(1);
 
-    // Find the plane that contains our input section.
+    ret.weight[0][0] = ret.weight[0][1] = 1; 
+    ret.weight[1][0] = ret.weight[1][1] = 1;
 
-    // Generate a polygon from the curves, and use this to test how many
-    // times each loop is enclosed. Then set the direction (cw/ccw) to
-    // be correct for outlines/holes, so that we generate correct normals.
-
-    // Now generate all the surfaces, top/bottom planes plus extrusions.
-   
-    // And now all the curves, trimming the top and bottom and their extrusion
-
-    // And the lines, trimming adjacent extrusion surfaces.
+    ret.ctrl[0][0] = pt;
+    ret.ctrl[0][1] = pt.Plus(u);
+    ret.ctrl[1][0] = pt.Plus(v);
+    ret.ctrl[1][1] = pt.Plus(v).Plus(u);
+    
     return ret;
 }
 
+Vector SSurface::PointAt(double u, double v) {
+    Vector num = Vector::From(0, 0, 0);
+    double den = 0;
+
+    int i, j;
+    for(i = 0; i <= degm; i++) {
+        for(j = 0; j <= degn; j++) {
+            double Bi = Bernstein(i, degm, u),
+                   Bj = Bernstein(j, degn, v);
+
+            num = num.Plus(ctrl[i][j].ScaledBy(Bi*Bj*weight[i][j]));
+            den += weight[i][j]*Bi*Bj;
+        }
+    }
+    num = num.ScaledBy(1.0/den);
+    return num;
+}
+
+void SSurface::TangentsAt(double u, double v, Vector *tu, Vector *tv) {
+    Vector num   = Vector::From(0, 0, 0),
+           num_u = Vector::From(0, 0, 0),
+           num_v = Vector::From(0, 0, 0);
+    double den   = 0,
+           den_u = 0,
+           den_v = 0;
+
+    int i, j;
+    for(i = 0; i <= degm; i++) {
+        for(j = 0; j <= degn; j++) {
+            double Bi  = Bernstein(i, degm, u),
+                   Bj  = Bernstein(j, degn, v),
+                   Bip = BernsteinDerivative(i, degm, u),
+                   Bjp = BernsteinDerivative(j, degn, v);
+
+            num = num.Plus(ctrl[i][j].ScaledBy(Bi*Bj*weight[i][j]));
+            den += weight[i][j]*Bi*Bj;
+
+            num_u = num_u.Plus(ctrl[i][j].ScaledBy(Bip*Bj*weight[i][j]));
+            den_u += weight[i][j]*Bip*Bj;
+
+            num_v = num_v.Plus(ctrl[i][j].ScaledBy(Bi*Bjp*weight[i][j]));
+            den_v += weight[i][j]*Bi*Bjp;
+        }
+    }
+    // Quotient rule; f(t) = n(t)/d(t), so f' = (n'*d - n*d')/(d^2)
+    *tu = ((num_u.ScaledBy(den)).Minus(num.ScaledBy(den_u)));
+    *tu = tu->ScaledBy(1.0/(den*den));
+
+    *tv = ((num_v.ScaledBy(den)).Minus(num.ScaledBy(den_v)));
+    *tv = tv->ScaledBy(1.0/(den*den));
+}
+
+Vector SSurface::NormalAt(double u, double v) {
+    Vector tu, tv;
+    TangentsAt(u, v, &tu, &tv);
+    return tu.Cross(tv);
+}
+
+void SSurface::ClosestPointTo(Vector p, double *u, double *v) {
+    int i, j;
+    double minDist = 1e10;
+    double res = 7.0;
+    for(i = 0; i < (int)res; i++) {
+        for(j = 0; j <= (int)res; j++) {
+            double tryu = (i/res), tryv = (j/res);
+            
+            Vector tryp = PointAt(tryu, tryv);
+            double d = (tryp.Minus(p)).Magnitude();
+            if(d < minDist) {
+                *u = tryu;
+                *v = tryv;
+                minDist = d;
+            }
+        }
+    }
+
+    // Initial guess is in u, v
+    Vector p0;
+    for(i = 0; i < 50; i++) {
+        p0 = PointAt(*u, *v);
+        if(p0.Equals(p)) {
+            return;
+        }
+
+        Vector tu, tv;
+        TangentsAt(*u, *v, &tu, &tv);
+        
+        // Project the point into a plane through p0, with basis tu, tv; a
+        // second-order thing would converge faster but needs second
+        // derivatives.
+        Vector dp = p.Minus(p0);
+        double du = dp.Dot(tu), dv = dp.Dot(tv);
+        *u += du / (tu.MagSquared());
+        *v += dv / (tu.MagSquared());
+    }
+    dbp("didn't converge");
+    dbp("have %.3f %.3f %.3f", CO(p0));
+    dbp("want %.3f %.3f %.3f", CO(p));
+    if(isnan(*u) || isnan(*v)) {
+        *u = *v = 0;
+    }
+}
+
+void SSurface::TriangulateInto(SShell *shell, SMesh *sm) {
+    SEdgeList el;
+    ZERO(&el);
+
+    STrimBy *stb;
+    for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
+        SCurve *sc = shell->curve.FindById(stb->curve);
+
+        Vector prevuv, ptuv;
+        bool inCurve = false;
+        Vector *pt;
+        double u = 0, v = 0;
+        for(pt = sc->pts.First(); pt; pt = sc->pts.NextAfter(pt)) {
+            ClosestPointTo(*pt, &u, &v);
+            ptuv = Vector::From(u, v, 0);
+
+            if(inCurve) {
+                el.AddEdge(prevuv, ptuv);
+            }
+            prevuv = ptuv;
+
+            if(pt->EqualsExactly(stb->start)) inCurve = true;
+            if(pt->EqualsExactly(stb->finish)) inCurve = false;
+        }
+    }
+    
+    SPolygon poly;
+    ZERO(&poly);
+    if(!el.AssemblePolygon(&poly, NULL)) {
+        dbp("failed to assemble polygon to trim nurbs surface in uv space");
+    }
+
+    int i, start = sm->l.n;
+    STriMeta meta = { 0, 0x888888 };
+    poly.normal = Vector::From(0, 0, 1);
+    poly.TriangulateInto(sm, meta);
+    for(i = start; i < sm->l.n; i++) {
+        STriangle *st = &(sm->l.elem[i]);
+        st->an = NormalAt(st->a.x, st->a.y);
+        st->bn = NormalAt(st->b.x, st->b.y);
+        st->cn = NormalAt(st->c.x, st->c.y);
+        st->a = PointAt(st->a.x, st->a.y);
+        st->b = PointAt(st->b.x, st->b.y);
+        st->c = PointAt(st->c.x, st->c.y);
+        if((st->Normal()).Dot(st->an) < 0) {
+            // Have to get the vertices in the right order
+            st->FlipNormal();
+        }
+    }
+    
+    el.Clear();
+    poly.Clear();
+}
+
+void SSurface::Clear(void) {
+    trim.Clear();
+}
+
+SShell SShell::FromExtrusionOf(SBezierLoopSet *sbls, Vector t0, Vector t1) {
+    SShell ret;
+    ZERO(&ret);
+
+    // Make the extrusion direction consistent with respect to the normal
+    // of the sketch we're extruding.
+    if((t0.Minus(t1)).Dot(sbls->normal) < 0) {
+        SWAP(Vector, t0, t1);
+    }
+
+    // First, generate the top and bottom surfaces of the extrusion; just
+    // planes.
+    SSurface s0, s1;
+    s0 = SSurface::FromPlane(sbls->point.Plus(t0), sbls->normal.ScaledBy(-1));
+    s1 = SSurface::FromPlane(sbls->point.Plus(t1), sbls->normal.ScaledBy( 1));
+    hSSurface hs0 = ret.surface.AddAndAssignId(&s0),
+              hs1 = ret.surface.AddAndAssignId(&s1);
+    
+    // Now go through the input curves. For each one, generate its surface
+    // of extrusion, its two translated trim curves, and one trim line. We
+    // go through by loops so that we can assign the lines correctly.
+    SBezierLoop *sbl;
+    for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
+        SBezier *sb;
+
+        typedef struct {
+            STrimBy     trim;
+            hSSurface   hs;
+        } TrimLine;
+        List<TrimLine> trimLines;
+        ZERO(&trimLines);
+
+        for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
+            // Generate the surface of extrusion of this curve, and add
+            // it to the list
+            SSurface ss = SSurface::FromExtrusionOf(sb, t0, t1);
+            hSSurface hsext = ret.surface.AddAndAssignId(&ss);
+
+            // Translate the curve by t0 and t1 to produce two trim curves
+            SCurve sc;
+            ZERO(&sc);
+            sb->MakePwlInto(&(sc.pts), t0);
+            hSCurve hc0 = ret.curve.AddAndAssignId(&sc);
+            STrimBy stb0 = STrimBy::EntireCurve(&ret, hc0);
+
+            ZERO(&sc);
+            sb->MakePwlInto(&(sc.pts), t1);
+            hSCurve hc1 = ret.curve.AddAndAssignId(&sc);
+            STrimBy stb1 = STrimBy::EntireCurve(&ret, hc1);
+
+            // The translated curves trim the flat top and bottom surfaces.
+            (ret.surface.FindById(hs0))->trim.Add(&stb0);
+            (ret.surface.FindById(hs1))->trim.Add(&stb1);
+
+            // The translated curves also trim the surface of extrusion.
+            (ret.surface.FindById(hsext))->trim.Add(&stb0);
+            (ret.surface.FindById(hsext))->trim.Add(&stb1);
+
+            // And form the trim line
+            Vector pt = sb->Finish();
+            Vector p0 = pt.Plus(t0), p1 = pt.Plus(t1);
+            ZERO(&sc);
+            sc.pts.Add(&p0);
+            sc.pts.Add(&p1);
+            hSCurve hl = ret.curve.AddAndAssignId(&sc);
+            // save this for later
+            TrimLine tl;
+            tl.trim = STrimBy::EntireCurve(&ret, hl);
+            tl.hs = hsext;
+            trimLines.Add(&tl);
+        }
+
+        int i;
+        for(i = 0; i < trimLines.n; i++) {
+            TrimLine *tl = &(trimLines.elem[i]);
+            SSurface *ss = ret.surface.FindById(tl->hs);
+
+            TrimLine *tlp = &(trimLines.elem[WRAP(i-1, trimLines.n)]);
+
+            ss->trim.Add(&(tl->trim));
+            ss->trim.Add(&(tlp->trim));
+        }
+        trimLines.Clear();
+    }
+
+    return ret;
+}
+
+void SShell::TriangulateInto(SMesh *sm) {
+    SSurface *s;
+    for(s = surface.First(); s; s = surface.NextAfter(s)) {
+        s->TriangulateInto(this, sm);
+    }
+}
+
+void SShell::Clear(void) {
+    SSurface *s;
+    for(s = surface.First(); s; s = surface.NextAfter(s)) {
+        s->Clear();
+    }
+    surface.Clear();
+
+    SCurve *c;
+    for(c = curve.First(); c; c = curve.NextAfter(c)) {
+        c->Clear();
+    }
+    curve.Clear();
+}
 
diff --git a/srf/surface.h b/srf/surface.h
index c510db2..9564f0b 100644
--- a/srf/surface.h
+++ b/srf/surface.h
@@ -6,6 +6,8 @@
 double Bernstein(int k, int deg, double t);
 double BernsteinDerivative(int k, int deg, double t);
 
+class SShell;
+
 class hSSurface {
 public:
     DWORD v;
@@ -29,7 +31,8 @@ public:
     Vector Start(void);
     Vector Finish(void);
     void MakePwlInto(List<Vector> *l);
-    void MakePwlWorker(List<Vector> *l, double ta, double tb);
+    void MakePwlInto(List<Vector> *l, Vector offset);
+    void MakePwlWorker(List<Vector> *l, double ta, double tb, Vector offset);
 
     void Reverse(void);
 
@@ -61,6 +64,7 @@ class SBezierLoopSet {
 public:
     List<SBezierLoop> l;
     Vector normal;
+    Vector point;
 
     static SBezierLoopSet From(SBezierList *spcl, SPolygon *poly,
                           bool *allClosed, SEdge *errorAt);
@@ -73,10 +77,12 @@ class SCurve {
 public:
     hSCurve         h;
 
-    SBezier         exact; // or deg = 0 if we don't know the exact form
+    bool            isExact;
+    SBezier         exact;
+
     List<Vector>    pts;
-    hSSurface       srfA;
-    hSSurface       srfB;
+
+    void Clear(void);
 };
 
 // A segment of a curve by which a surface is trimmed: indicates which curve,
@@ -89,8 +95,11 @@ public:
     Vector      start;
     Vector      finish;
     Vector      out;
+
+    static STrimBy STrimBy::EntireCurve(SShell *shell, hSCurve hsc);
 };
 
+// A rational polynomial surface in Bezier form.
 class SSurface {
 public:
     hSSurface       h;
@@ -102,14 +111,16 @@ public:
     List<STrimBy>   trim;
 
     static SSurface FromExtrusionOf(SBezier *spc, Vector t0, Vector t1);
+    static SSurface FromPlane(Vector pt, Vector n);
 
     void ClosestPointTo(Vector p, double *u, double *v);
     Vector PointAt(double u, double v);
-    Vector TangentWrtUAt(double u, double v);
-    Vector TangentWrtVAt(double u, double v);
+    void TangentsAt(double u, double v, Vector *tu, Vector *tv);
     Vector NormalAt(double u, double v);
 
-    void TriangulateInto(SMesh *sm);
+    void TriangulateInto(SShell *shell, SMesh *sm);
+
+    void Clear(void);
 };
 
 class SShell {
@@ -117,9 +128,12 @@ public:
     IdList<SCurve,hSCurve>      curve;
     IdList<SSurface,hSSurface>  surface;
 
-    static SShell FromExtrusionOf(SBezierList *spcl, Vector t0, Vector t1);
+    static SShell FromExtrusionOf(SBezierLoopSet *sbls, Vector t0, Vector t1);
 
     static SShell FromUnionOf(SShell *a, SShell *b);
+
+    void TriangulateInto(SMesh *sm);
+    void Clear(void);
 };
 
 #endif
diff --git a/srf/triangulate.cpp b/srf/triangulate.cpp
new file mode 100644
index 0000000..e55a327
--- /dev/null
+++ b/srf/triangulate.cpp
@@ -0,0 +1,3 @@
+#include "../solvespace.h"
+
+
diff --git a/ui.h b/ui.h
index 3b81f99..9f78f5b 100644
--- a/ui.h
+++ b/ui.h
@@ -240,6 +240,7 @@ public:
         MNU_COMMENT,
         // Analyze
         MNU_VOLUME,
+        MNU_NAKED_EDGES,
         MNU_SHOW_DOF,
         MNU_TRACE_PT,
         MNU_STOP_TRACING,
diff --git a/undoredo.cpp b/undoredo.cpp
index 4648c0b..dd48ee7 100644
--- a/undoredo.cpp
+++ b/undoredo.cpp
@@ -53,6 +53,8 @@ void SolveSpace::PushFromCurrentOnto(UndoStack *uk) {
         ZERO(&(dest.polyError));
         ZERO(&(dest.thisMesh));
         ZERO(&(dest.runningMesh));
+        ZERO(&(dest.thisShell));
+        ZERO(&(dest.runningShell));
         ZERO(&(dest.meshError));
         ZERO(&(dest.emphEdges));
 
@@ -95,6 +97,8 @@ void SolveSpace::PopOntoCurrentFrom(UndoStack *uk) {
         g->bezierLoopSet.Clear();
         g->thisMesh.Clear();
         g->runningMesh.Clear();
+        g->thisShell.Clear();
+        g->runningShell.Clear();
         g->meshError.interferesAt.Clear();
         g->emphEdges.Clear();
         g->remap.Clear();