Many changes:

* Rewrite surface handles in curves, so that Booleans beyond
      the first don't screw up.

    * If an intersection curve is identical to an existing curve
      (as happens when faces are coincident), take the piecewise
      linearization of the existing curve; this stops us from
      screwing up when different shells are pwl'd at different
      chord tols.

    * Hook up the plane faces again.

    * Remove coincident (parallel or anti-parallel) edges from the
      coincident-face edge lists when doing Booleans; those may
      happen if two faces are coincident with ours.

    * Miscellaneous bugfixes.

It doesn't seem to screw up very much now, although tangent edges
(and insufficient pwl resolution) may still cause problems.

[git-p4: depot-paths = "//depot/solvespace/": change = 1929]

draw.cpp

@@ -1023,14 +1023,8 @@ void GraphicsWindow::Paint(int w, int h) {
 
     // And the naked edges, if the user did Analyze -> Show Naked Edges.
     glLineWidth(7);
-    glEnable(GL_LINE_STIPPLE);
-    glLineStipple(1, 0x5555);
     glColor3d(1, 0, 0);
     glxDrawEdges(&(SS.nakedEdges));
-    glLineStipple(1, 0xaaaa);
-    glColor3d(0, 0, 0);
-    glxDrawEdges(&(SS.nakedEdges));
-    glDisable(GL_LINE_STIPPLE);
 
     // Then redraw whatever the mouse is hovering over, highlighted.
     glDisable(GL_DEPTH_TEST); 

dsc.h

@@ -15,6 +15,8 @@ public:
     // a + (vx)*i + (vy)*j + (vz)*k
     double w, vx, vy, vz;
 
+    static const Quaternion IDENTITY;
+
     static Quaternion From(double w, double vx, double vy, double vz);
     static Quaternion From(hParam w, hParam vx, hParam vy, hParam vz);
     static Quaternion From(Vector u, Vector v);
@@ -54,8 +56,7 @@ public:
                                                bool *parallel);
     static Vector AtIntersectionOfPlanes(Vector na, double da,
                                          Vector nb, double db,
-                                         Vector nc, double dc,
+                                         Vector nc, double dc, bool *parallel);
-                                         bool *parallel);
 
     double Element(int i);
     bool Equals(Vector v, double tol=LENGTH_EPS);

graphicswin.cpp

@@ -307,8 +307,8 @@ void GraphicsWindow::ZoomToFit(void) {
         if(dy != 0) scaley = 0.9*height/dy;
         scale = min(scalex, scaley);
 
-        scale = min(100, scale);
+        scale = min(300, scale);
-        scale = max(0.001, scale);
+        scale = max(0.003, scale);
     }
 
     // Then do another run, considering the perspective.

groupmesh.cpp

@@ -55,6 +55,12 @@ void Group::GenerateShellForStepAndRepeat(void) {
     Group *src = SS.GetGroup(opA);
     SShell *srcs = &(src->thisShell); // the shell to step and repeat
 
+    SShell workA, workB;
+    ZERO(&workA);
+    ZERO(&workB);
+    SShell *soFar = &workA, *scratch = &workB;
+    soFar->MakeFromCopyOf(src->PreviousGroupShell());
+
     int n = (int)valA, a0 = 0;
     if(subtype == ONE_SIDED && skipFirst) {
         a0++; n++;
@@ -64,25 +70,37 @@ void Group::GenerateShellForStepAndRepeat(void) {
         int ap = a*2 - (subtype == ONE_SIDED ? 0 : (n-1));
         int remap = (a == (n - 1)) ? REMAP_LAST : a;
 
+        SShell transd;
+        ZERO(&transd);
         if(type == TRANSLATE) {
             Vector trans = Vector::From(h.param(0), h.param(1), h.param(2));
             trans = trans.ScaledBy(ap);
-
+            Quaternion q = Quaternion::From(1, 0, 0, 0);
+            transd.MakeFromTransformationOf(srcs, trans, q);
         } else {
             Vector trans = Vector::From(h.param(0), h.param(1), h.param(2));
             double theta = ap * SS.GetParam(h.param(3))->val;
             double c = cos(theta), s = sin(theta);
             Vector axis = Vector::From(h.param(4), h.param(5), h.param(6));
             Quaternion q = Quaternion::From(c, s*axis.x, s*axis.y, s*axis.z);
-
+            // Rotation is centered at t; so A(x - t) + t = Ax + (t - At)
+            transd.MakeFromTransformationOf(srcs,
+                trans.Minus(q.Rotate(trans)), q);
         }
 
         if(src->meshCombine == COMBINE_AS_DIFFERENCE) {
-
+            scratch->MakeFromDifferenceOf(soFar, &transd);
         } else {
+            scratch->MakeFromUnionOf(soFar, &transd);
+        }
+        SWAP(SShell *, scratch, soFar);
 
+        scratch->Clear();
+        transd.Clear();
     }
-    }
+
+    runningShell.Clear();
+    runningShell = *soFar;
 }
 
 void Group::GenerateShellAndMesh(void) {
@@ -105,6 +123,55 @@ void Group::GenerateShellAndMesh(void) {
         }
         
         thisShell.MakeFromExtrusionOf(&(src->bezierLoopSet), tbot, ttop, color);
+        Vector onOrig = src->bezierLoopSet.point;
+        // And for any plane faces, annotate the model with the entity for
+        // that face, so that the user can select them with the mouse.
+        int i;
+        for(i = 0; i < thisShell.surface.n; i++) {
+            SSurface *ss = &(thisShell.surface.elem[i]);
+            hEntity face = Entity::NO_ENTITY;
+
+            Vector p = ss->PointAt(0, 0),
+                   n = ss->NormalAt(0, 0).WithMagnitude(1);
+            double d = n.Dot(p);
+
+            if(i == 0 || i == 1) {
+                // These are the top and bottom of the shell.
+                if(fabs((onOrig.Plus(ttop)).Dot(n) - d) < LENGTH_EPS) {
+                    face = Remap(Entity::NO_ENTITY, REMAP_TOP);
+                    ss->face = face.v;
+                }
+                if(fabs((onOrig.Plus(tbot)).Dot(n) - d) < LENGTH_EPS) {
+                    face = Remap(Entity::NO_ENTITY, REMAP_BOTTOM);
+                    ss->face = face.v;
+                }
+                continue;
+            }
+
+            // So these are the sides
+            if(ss->degm != 1 || ss->degn != 1) continue;
+
+            Entity *e;
+            for(e = SS.entity.First(); e; e = SS.entity.NextAfter(e)) {
+                if(e->group.v != opA.v) continue;
+                if(e->type != Entity::LINE_SEGMENT) continue;
+
+                Vector a = SS.GetEntity(e->point[0])->PointGetNum(),
+                       b = SS.GetEntity(e->point[1])->PointGetNum();
+                a = a.Plus(ttop);
+                b = b.Plus(ttop);
+                // Could get taken backwards, so check all cases.
+                if((a.Equals(ss->ctrl[0][0]) && b.Equals(ss->ctrl[1][0])) ||
+                   (b.Equals(ss->ctrl[0][0]) && a.Equals(ss->ctrl[1][0])) ||
+                   (a.Equals(ss->ctrl[0][1]) && b.Equals(ss->ctrl[1][1])) ||
+                   (b.Equals(ss->ctrl[0][1]) && a.Equals(ss->ctrl[1][1])))
+                {
+                    face = Remap(e->h, REMAP_LINE_TO_FACE);
+                    ss->face = face.v;
+                    break;
+                }
+            }
+        }
     } else if(type == LATHE) {
         Group *src = SS.GetGroup(opA);
 
@@ -138,7 +205,6 @@ void Group::GenerateShellAndMesh(void) {
         }
     }
 
-    runningMesh.Clear();
     runningShell.Clear();
 
     // If this group contributes no new mesh, then our running mesh is the
@@ -174,6 +240,7 @@ void Group::GenerateShellAndMesh(void) {
     }
 
 done:
+    runningMesh.Clear();
     runningShell.TriangulateInto(&runningMesh);
     emphEdges.Clear();
     if(h.v == SS.GW.activeGroup.v && SS.edgeColor != 0) {

polygon.cpp

@@ -222,6 +222,31 @@ intersects:
     return cnt;
 }
 
+//-----------------------------------------------------------------------------
+// Remove unnecessary edges: if two are anti-parallel then remove both, and if
+// two are parallel then remove one.
+//-----------------------------------------------------------------------------
+void SEdgeList::CullExtraneousEdges(void) {
+    l.ClearTags();
+    int i, j;
+    for(i = 0; i < l.n; i++) {
+        SEdge *se = &(l.elem[i]);
+        for(j = i+1; j < l.n; j++) {
+            SEdge *set = &(l.elem[j]);
+            if((set->a).Equals(se->a) && (set->b).Equals(se->b)) {
+                // Two parallel edges exist; so keep only the first one.
+                set->tag = 1;
+            }
+            if((set->a).Equals(se->b) && (set->b).Equals(se->a)) {
+                // Two anti-parallel edges exist; so keep neither.
+                se->tag = 1;
+                set->tag = 1;
+            }
+        }
+    }
+    l.RemoveTagged();
+}
+
 void SContour::AddPoint(Vector p) {
     SPoint sp;
     sp.tag = 0;

polygon.h

@@ -26,6 +26,7 @@ public:
     bool AssembleContour(Vector first, Vector last, SContour *dest,
                             SEdge *errorAt, bool keepDir);
     int AnyEdgeCrossings(Vector a, Vector b, Vector *pi=NULL);
+    void CullExtraneousEdges(void);
 };
 
 class SPoint {

srf/boolean.cpp

@@ -302,8 +302,12 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
     SEdgeList sameNormal, oppositeNormal;
     ZERO(&sameNormal);
     ZERO(&oppositeNormal);
-    agnst->MakeCoincidentEdgesInto(&ret, true, &sameNormal);
+    agnst->MakeCoincidentEdgesInto(&ret, true, &sameNormal, into);
-    agnst->MakeCoincidentEdgesInto(&ret, false, &oppositeNormal);
+    agnst->MakeCoincidentEdgesInto(&ret, false, &oppositeNormal, into);
+    // and cull parallel or anti-parallel pairs, which may occur if multiple
+    // surfaces are coincident with ours
+    sameNormal.CullExtraneousEdges();
+    oppositeNormal.CullExtraneousEdges();
     // and build the trees for quick in-polygon testing
     SBspUv *sameBsp = SBspUv::From(&sameNormal);
     SBspUv *oppositeBsp = SBspUv::From(&oppositeNormal);
@@ -318,9 +322,9 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
         for(sc = into->curve.First(); sc; sc = into->curve.NextAfter(sc)) {
             if(sc->source != SCurve::FROM_INTERSECTION) continue;
             if(opA) {
-                if(sc->surfB.v != h.v || sc->surfA.v != ss->h.v) continue;
-            } else {
                 if(sc->surfA.v != h.v || sc->surfB.v != ss->h.v) continue;
+            } else {
+                if(sc->surfB.v != h.v || sc->surfA.v != ss->h.v) continue;
             }
            
             int i;
@@ -403,11 +407,6 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
 
         if(KeepEdge(type, opA, tag)) {
             final.AddEdge(se->a, se->b, se->auxA, se->auxB);
-        } else {
-            if(I == 1) {
-                dbp("orig vs. shell: %d (l=%g)",
-                    c_shell, ((se->b).Minus(se->a)).Magnitude());
-            }
         }
     }
 
@@ -436,33 +435,23 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
         }
     }
 
-    // Cull extraneous edges; duplicates or anti-parallel pairs
+    // Cull extraneous edges; duplicates or anti-parallel pairs. In particular,
-    final.l.ClearTags();
+    // we can get duplicate edges if our surface intersects the other shell
-    int i, j;
+    // at an edge, so that both surfaces intersect coincident (and both
-    for(i = 0; i < final.l.n; i++) {
+    // generate an intersection edge).
-        se = &(final.l.elem[i]);
+    final.CullExtraneousEdges();
-        for(j = i+1; j < final.l.n; j++) {
-            SEdge *set = &(final.l.elem[j]);
-            if((set->a).Equals(se->a) && (set->b).Equals(se->b)) {
-                // Two parallel edges exist; so keep only the first one. This
-                // can happen if our surface intersects the shell at an edge,
-                // so that we get two copies of the intersection edge.
-                set->tag = 1;
-            }
-            if((set->a).Equals(se->b) && (set->b).Equals(se->a)) {
-                // Two anti-parallel edges exist; so keep neither.
-                se->tag = 1;
-                set->tag = 1;
-            }
-        }
-    }
-    final.l.RemoveTagged();
-
-//    if(I == 0) DEBUGEDGELIST(&final, &ret);
 
     // Use our reassembled edges to trim the new surface.
     ret.TrimFromEdgeList(&final);
 
+    SPolygon poly;
+    ZERO(&poly);
+    final.l.ClearTags();
+    if(!final.AssemblePolygon(&poly, NULL, true)) {
+        DEBUGEDGELIST(&inter, &ret);
+    }
+    poly.Clear();
+
     sameNormal.Clear();
     oppositeNormal.Clear();
     final.Clear();
@@ -478,20 +467,20 @@ void SShell::CopySurfacesTrimAgainst(SShell *against, SShell *into,
     for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
         SSurface ssn;
         ssn = ss->MakeCopyTrimAgainst(against, this, into, type, opA);
-        into->surface.AddAndAssignId(&ssn);
+        ss->newH = into->surface.AddAndAssignId(&ssn);
         I++;
     }
 }
 
 void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) {
     SSurface *sa;
-    for(sa = agnst->surface.First(); sa; sa = agnst->surface.NextAfter(sa)) {
+    for(sa = surface.First(); sa; sa = surface.NextAfter(sa)) {
         SSurface *sb;
-        for(sb = surface.First(); sb; sb = surface.NextAfter(sb)) {
+        for(sb = agnst->surface.First(); sb; sb = agnst->surface.NextAfter(sb)){
             // Intersect every surface from our shell against every surface
             // from agnst; this will add zero or more curves to the curve
             // list for into.
-            sa->IntersectAgainst(sb, agnst, this, into);
+            sa->IntersectAgainst(sb, this, agnst, into);
         }
         FLAG++;
     }
@@ -517,17 +506,34 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
     // the surfaces in B (which is all of the intersection curves).
     a->MakeIntersectionCurvesAgainst(b, this);
     
-    I = 100;
     if(b->surface.n == 0 || a->surface.n == 0) {
         // Then trim and copy the surfaces
         a->CopySurfacesTrimAgainst(b, this, type, true);
         b->CopySurfacesTrimAgainst(a, this, type, false);
     } else {
-        I = 0;
         a->CopySurfacesTrimAgainst(b, this, type, true);
         b->CopySurfacesTrimAgainst(a, this, type, false);
     }
 
+    // Now that we've copied the surfaces, we know their new hSurfaces, so
+    // rewrite the curves to refer to the surfaces by their handles in the
+    // result.
+    SCurve *sc;
+    for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) {
+        if(sc->source == SCurve::FROM_A) {
+            sc->surfA = a->surface.FindById(sc->surfA)->newH;
+            sc->surfB = a->surface.FindById(sc->surfB)->newH;
+        } else if(sc->source == SCurve::FROM_B) {
+            sc->surfA = b->surface.FindById(sc->surfA)->newH;
+            sc->surfB = b->surface.FindById(sc->surfB)->newH;
+        } else if(sc->source == SCurve::FROM_INTERSECTION) {
+            sc->surfA = a->surface.FindById(sc->surfA)->newH;
+            sc->surfB = b->surface.FindById(sc->surfB)->newH;
+        } else {
+            oops();
+        }
+    }
+
     // And clean up the piecewise linear things we made as a calculation aid
     a->CleanupAfterBoolean();
     b->CleanupAfterBoolean();
@@ -682,6 +688,13 @@ int SBspUv::ClassifyPoint(Point2d p, Point2d eb) {
 }
 
 int SBspUv::ClassifyEdge(Point2d ea, Point2d eb) {
-    return ClassifyPoint((ea.Plus(eb)).ScaledBy(0.5), eb);
+    int ret = ClassifyPoint((ea.Plus(eb)).ScaledBy(0.5), eb);
+    if(ret == EDGE_OTHER) {
+        // Perhaps the edge is tangent at its midpoint (and we screwed up
+        // somewhere earlier and failed to split it); try a different
+        // point on the edge.
+        ret = ClassifyPoint(ea.Plus((eb.Minus(ea)).ScaledBy(0.294)), eb);
+    }
+    return ret;
 }
 

srf/ratpoly.cpp

@@ -144,15 +144,10 @@ Vector SBezier::PointAt(double t) {
 }
 
 void SBezier::MakePwlInto(List<Vector> *l) {
-    MakePwlInto(l, Vector::From(0, 0, 0));
+    l->Add(&(ctrl[0]));
+    MakePwlWorker(l, 0.0, 1.0);
 }
-void SBezier::MakePwlInto(List<Vector> *l, Vector offset) {
+void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb) {
-    Vector p = (ctrl[0]).Plus(offset);
-    l->Add(&p);
-
-    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);
 
@@ -169,12 +164,11 @@ void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb, Vector off) {
     double step = 1.0/SS.maxSegments;
     if((tb - ta) < step || d < SS.ChordTolMm()) {
         // 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, off);
+        MakePwlWorker(l, ta, tm);
-        MakePwlWorker(l, tm, tb, off);
+        MakePwlWorker(l, tm, tb);
     }
 }
 
@@ -206,6 +200,18 @@ SBezier SBezier::TransformedBy(Vector t, Quaternion q) {
     return ret;
 }
 
+bool SBezier::Equals(SBezier *b) {
+    // We just test of identical degree and control points, even though two
+    // curves could still be coincident (even sharing endpoints).
+    if(deg != b->deg) return false;
+    int i;
+    for(i = 0; i <= deg; i++) {
+        if(!(ctrl[i]).Equals(b->ctrl[i])) return false;
+        if(fabs(weight[i] - b->weight[i]) > LENGTH_EPS) return false;
+    }
+    return true;
+}
+
 void SBezierList::Clear(void) {
     l.Clear();
 }
@@ -375,6 +381,8 @@ SCurve SCurve::FromTransformationOf(SCurve *a, Vector t, Quaternion q) {
     ret.h = a->h;
     ret.isExact = a->isExact;
     ret.exact = (a->exact).TransformedBy(t, q);
+    ret.surfA = a->surfA;
+    ret.surfB = a->surfB;
     
     Vector *p;
     for(p = a->pts.First(); p; p = a->pts.NextAfter(p)) {
@@ -454,7 +462,7 @@ bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) {
 }
 
 bool SSurface::IsCylinder(Vector *center, Vector *axis, double *r,
-                             double *dtheta)
+                             Vector *start, Vector *finish)
 {
     SBezier sb;
     if(!IsExtrusion(&sb, axis)) return false;
@@ -480,10 +488,13 @@ bool SSurface::IsCylinder(Vector *center, Vector *axis, double *r,
             pb2 = (sb.ctrl[2]).Project2d(u, v).Minus(c2);
     
     double thetaa = atan2(pa2.y, pa2.x), // in fact always zero due to csys
-           thetab = atan2(pb2.y, pb2.x);
+           thetab = atan2(pb2.y, pb2.x),
-    *dtheta = WRAP_NOT_0(thetab - thetaa, 2*PI);
+           dtheta = WRAP_NOT_0(thetab - thetaa, 2*PI);
 
-    if(fabs(sb.weight[1] - cos(*dtheta/2)) > LENGTH_EPS) return false;
+    if(fabs(sb.weight[1] - cos(dtheta/2)) > LENGTH_EPS) return false;
+
+    *start  = (sb.ctrl[0]).Minus(*center);
+    *finish = (sb.ctrl[2]).Minus(*center);
 
     return true;
 }
@@ -743,11 +754,20 @@ void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) {
     }
 }
 
-void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv) {
+void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv,
+                             SShell *useCurvesFrom) {
     STrimBy *stb;
     for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
         SCurve *sc = shell->curve.FindById(stb->curve);
 
+        // We have the option to use the curves from another shell; this
+        // is relevant when generating the coincident edges while doing the
+        // Booleans, since the curves from the output shell will be split
+        // against any intersecting surfaces (and the originals aren't).
+        if(useCurvesFrom) {
+            sc = useCurvesFrom->curve.FindById(sc->newH);
+        }
+
         Vector prev, prevuv, ptuv;
         bool inCurve = false, empty = true;
         double u = 0, v = 0;
@@ -783,9 +803,8 @@ void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv) {
             if(pt->Equals(stb->start)) inCurve = true;
             if(pt->Equals(stb->finish)) inCurve = false;
         }
-        if(inCurve || empty) {
+        if(inCurve) dbp("trim was unterminated");
-            dbp("trim was empty or unterminated");
+        if(empty)   dbp("trim was empty");
-        }
     }
 }
 
@@ -910,13 +929,17 @@ void SShell::MakeFromExtrusionOf(SBezierLoopSet *sbls, Vector t0, Vector t1,
             // Translate the curve by t0 and t1 to produce two trim curves
             SCurve sc;
             ZERO(&sc);
-            sb->MakePwlInto(&(sc.pts), t0);
+            sc.isExact = true;
+            sc.exact = sb->TransformedBy(t0, Quaternion::IDENTITY);
+            (sc.exact).MakePwlInto(&(sc.pts));
             sc.surfA = hs0;
             sc.surfB = hsext;
             hSCurve hc0 = curve.AddAndAssignId(&sc);
 
             ZERO(&sc);
-            sb->MakePwlInto(&(sc.pts), t1);
+            sc.isExact = true;
+            sc.exact = sb->TransformedBy(t1, Quaternion::IDENTITY);
+            (sc.exact).MakePwlInto(&(sc.pts));
             sc.surfA = hs1;
             sc.surfB = hsext;
             hSCurve hc1 = curve.AddAndAssignId(&sc);
@@ -936,10 +959,10 @@ void SShell::MakeFromExtrusionOf(SBezierLoopSet *sbls, Vector t0, Vector t1,
 
             // 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.isExact = true;
-            sc.pts.Add(&p1);
+            sc.exact = SBezier::From(pt.Plus(t0), pt.Plus(t1));
+            (sc.exact).MakePwlInto(&(sc.pts));
             hSCurve hl = curve.AddAndAssignId(&sc);
             // save this for later
             TrimLine tl;
@@ -969,10 +992,7 @@ void SShell::MakeFromExtrusionOf(SBezierLoopSet *sbls, Vector t0, Vector t1,
 }
 
 void SShell::MakeFromCopyOf(SShell *a) {
-    Vector t = Vector::From(0, 0, 0);
+    MakeFromTransformationOf(a, Vector::From(0, 0, 0), Quaternion::IDENTITY);
-    Quaternion q = Quaternion::From(1, 0, 0, 0);
-
-    MakeFromTransformationOf(a, t, q);
 }
 
 void SShell::MakeFromTransformationOf(SShell *a, Vector t, Quaternion q) {

srf/surface.h

@@ -61,9 +61,9 @@ public:
     Vector PointAt(double t);
     Vector Start(void);
     Vector Finish(void);
+    bool Equals(SBezier *b);
     void MakePwlInto(List<Vector> *l);
-    void MakePwlInto(List<Vector> *l, Vector offset);
+    void MakePwlWorker(List<Vector> *l, double ta, double tb);
-    void MakePwlWorker(List<Vector> *l, double ta, double tb, Vector offset);
 
     void GetBoundingProjd(Vector u, Vector orig, double *umin, double *umax);
     void Reverse(void);
@@ -170,6 +170,10 @@ class SSurface {
 public:
     hSSurface       h;
 
+    // Same as newH for the curves; record what a surface gets renamed to
+    // when I copy things over.
+    hSSurface       newH;
+
     int             color;
     DWORD           face;
 
@@ -222,10 +226,12 @@ public:
     bool CoincidentWithPlane(Vector n, double d);
     bool CoincidentWith(SSurface *ss, bool sameNormal);
     bool IsExtrusion(SBezier *of, Vector *along);
-    bool IsCylinder(Vector *center, Vector *axis, double *r, double *dtheta);
+    bool IsCylinder(Vector *center, Vector *axis, double *r,
+            Vector *start, Vector *finish);
 
     void TriangulateInto(SShell *shell, SMesh *sm);
-    void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv);
+    void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv,
+            SShell *useCurvesFrom=NULL);
     void MakeClassifyingBsp(SShell *shell);
     double ChordToleranceForEdge(Vector a, Vector b);
 
@@ -254,7 +260,7 @@ public:
     void AllPointsIntersecting(Vector a, Vector b, List<SInter> *il,
                                 bool seg, bool trimmed);
     void MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal,
-                                 SEdgeList *el);
+                                 SEdgeList *el, SShell *useCurvesFrom);
     void CleanupAfterBoolean(void);
 
     static const int INSIDE            = 100;

srf/surfinter.cpp

@@ -7,13 +7,47 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
 {
     SCurve sc;
     ZERO(&sc);
+    // Important to keep the order of (surfA, surfB) consistent; when we later
+    // rewrite the identifiers, we rewrite surfA from A and surfB from B.
     sc.surfA = h;
     sc.surfB = srfB->h;
-    sb->MakePwlInto(&(sc.pts));
+    sc.exact = *sb;
+    sc.isExact = true;
 
-    // Now split the line where it intersects our existing surfaces
+    // Now we have to piecewise linearize the curve. If there's already an
-    SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, srfB);
+    // identical curve in the shell, then follow that pwl exactly, otherwise
+    // calculate from scratch.
+    SCurve split, *existing = NULL, *se;
+    SBezier sbrev = *sb;
+    sbrev.Reverse();
+    bool backwards = false;
+    for(se = into->curve.First(); se; se = into->curve.NextAfter(se)) {
+        if(se->isExact) {
+            if(sb->Equals(&(se->exact))) {
+                existing = se;
+                break;
+            }
+            if(sbrev.Equals(&(se->exact))) {
+                existing = se;
+                backwards = true;
+                break;
+            }
+        }
+    }
+    if(existing) {
+        Vector *v;
+        for(v = existing->pts.First(); v; v = existing->pts.NextAfter(v)) {
+            sc.pts.Add(v);
+        }
+        if(backwards) sc.pts.Reverse();
+        split = sc;
+        ZERO(&sc);
+    } else {
+        sb->MakePwlInto(&(sc.pts));
+        // and split the line where it intersects our existing surfaces
+        split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, srfB);
         sc.Clear();
+    }
    
     if(0 && sb->deg == 1) {
         dbp(" ");
@@ -26,6 +60,8 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
             prev = v;
         }
     }
+    // Nothing should be generating zero-len edges.
+    if((sb->Start()).Equals(sb->Finish())) oops();
 
     split.source = SCurve::FROM_INTERSECTION;
     into->curve.AddAndAssignId(&split);
@@ -98,7 +134,7 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
             }
         }
 
-        if(tmax > tmin) {
+        if(tmax > tmin + LENGTH_EPS) {
             SBezier bezier = SBezier::From(p.Plus(dl.ScaledBy(tmin)),
                                            p.Plus(dl.ScaledBy(tmax)));
             AddExactIntersectionCurve(&bezier, b, agnstA, agnstB, into);
@@ -415,8 +451,8 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
 
     // All the intersections between the line and the surface; either special
     // cases that we can quickly solve in closed form, or general numerical.
-    Vector center, axis;
+    Vector center, axis, start, finish;
-    double radius, dtheta;
+    double radius;
     if(degm == 1 && degn == 1) {
         // Against a plane, easy.
         Vector n = NormalAt(0, 0).WithMagnitude(1);
@@ -432,7 +468,7 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
             ClosestPointTo(p, &(inter.p.x), &(inter.p.y));
             inters.Add(&inter);
         }
-    } else if(IsCylinder(&center, &axis, &radius, &dtheta) && 0) {
+    } else if(IsCylinder(&center, &axis, &radius, &start, &finish) && 0) {
         // XXX, cylinder is easy in closed form
     } else {
         // General numerical solution by subdivision, fallback
@@ -564,8 +600,9 @@ int SShell::ClassifyPoint(Vector p, Vector pout) {
         // then our ray always lies on edge, and that's okay. Otherwise
         // try again in a different random direction.
         if((edge_inters == 2) || !onEdge) break;
-        if(cnt++ > 20) {
+        if(cnt++ > 5) {
             dbp("can't find a ray that doesn't hit on edge!");
+            dbp("on edge = %d, edge_inters = %d", onEdge, edge_inters);
             break;
         }
     }
@@ -626,12 +663,12 @@ bool SSurface::CoincidentWithPlane(Vector n, double d) {
 // the prototype surface.
 //-----------------------------------------------------------------------------
 void SShell::MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal,
-                                     SEdgeList *el)
+                                     SEdgeList *el, SShell *useCurvesFrom)
 {
     SSurface *ss;
     for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
         if(proto->CoincidentWith(ss, sameNormal)) {
-            ss->MakeEdgesInto(this, el, false);
+            ss->MakeEdgesInto(this, el, false, useCurvesFrom);
         }
     }
 

util.cpp

@@ -100,6 +100,8 @@ void MakeMatrix(double *mat, double a11, double a12, double a13, double a14,
     mat[15] = a44;
 }
 
+const Quaternion Quaternion::IDENTITY = { 1, 0, 0, 0 };
+
 Quaternion Quaternion::From(double w, double vx, double vy, double vz) {
     Quaternion q;
     q.w  = w;

wishlist.txt

@@ -3,17 +3,16 @@ marching algorithm for surface intersection
 surfaces of revolution (lathed)
 cylinder-line special cases
 exact boundaries when near pwl trim
-step and repeat rotate/translate
 tangent intersections
 short pwl edge avoidance
-faces
-take consistent pwl with coincident faces
 exact curve export (at least for dxf)
 hidden line removal from mesh
 line styles (color, thickness)
+assembly
 
 -----
 loop detection
 incremental regen of entities?
+IGES and STEP export