Rewrite equiations generated for pt-on-line constraints.

Before this commit, pt-on-line constraints are buggy. To reproduce,
extrude a circle, then add a datum point and constrain it to the
axis of the circle, then move it. The cylinder will collapse.

To quote Jonathan:

> On investigation, I (a) confirm that the problem is
> the unconstrained extrusion depth going to zero, and (b) retract
> my earlier statement blaming extrude and other similar non-entity
> parameter treatment for this problem; you can easily reproduce it
> with a point in 3d constrained to lie on any line whose length
> is free.
>
> PT_ON_LINE is written using VectorsParallel, for no obvious reason.
> Rewriting that constraint to work on two projected distances (using
> any two basis vectors perpendicular to the line) should fix that
> problem, since replacing the "point on line in 3d" constraint with
> two "point on line in 2d" constraints works. That still has
> the hairy ball problem of choosing the basis vectors, which you
> can't do with a continuous function; you'd need Vector::Normal()
> or equivalent.
>
> You could write three equations and make the constraint itself
> introduce one new parameter for t. I don't know how well that
> would work numerically, but it would avoid the hairy ball problem,
> perhaps elegant at the cost of speed.

Indeed, this commit implements the latter solution: it introduces
an additional parameter, and migrates any old files that do not
include this parameter.

CHANGELOG.md

@@ -38,6 +38,10 @@ Other new features:
     workplane is displayed.
   * The "=" key is bound to "Zoom In", like "+" key.
 
+Bugs fixed:
+  * A point in 3d constrained to any line whose length is free no longer
+    causes the line length to collapse.
+
 2.3
 ---
 

src/constrainteq.cpp

@@ -182,7 +182,7 @@ void ConstraintBase::ModifyToSatisfy() {
         // that means no extra work.
         IdList<Equation,hEquation> l = {};
         // Generate the equations even if this is a reference dimension
-        GenerateReal(&l);
+        GenerateEquations(&l, /*forReference=*/true);
         ssassert(l.n == 1, "Expected constraint to generate a single equation");
 
         // These equations are written in the form f(...) - d = 0, where
@@ -201,14 +201,25 @@ void ConstraintBase::AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index)
     l->Add(&eq);
 }
 
-void ConstraintBase::Generate(IdList<Equation,hEquation> *l) const {
-    if(!reference) {
-        GenerateReal(l);
+void ConstraintBase::Generate(IdList<Param,hParam> *l) const {
+    switch(type) {
+        case Type::PT_ON_LINE: {
+            Param p = {};
+            p.h = h.param(0);
+            l->Add(&p);
+            break;
+        }
+
+        default:
+            break;
     }
 }
-void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) const {
-    Expr *exA = Expr::From(valA);
 
+void ConstraintBase::GenerateEquations(IdList<Equation,hEquation> *l,
+                                       bool forReference) const {
+    if(reference && !forReference) return;
+
+    Expr *exA = Expr::From(valA);
     switch(type) {
         case Type::PT_PT_DISTANCE:
             AddEq(l, Distance(workplane, ptA, ptB)->Minus(exA), 0);
@@ -382,37 +393,24 @@ void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) const {
             return;
         }
 
-        case Type::PT_ON_LINE:
-            if(workplane.v == EntityBase::FREE_IN_3D.v) {
-                EntityBase *ln = SK.GetEntity(entityA);
-                EntityBase *a = SK.GetEntity(ln->point[0]);
-                EntityBase *b = SK.GetEntity(ln->point[1]);
-                EntityBase *p = SK.GetEntity(ptA);
-
-                ExprVector ep = p->PointGetExprs();
-                ExprVector ea = a->PointGetExprs();
-                ExprVector eb = b->PointGetExprs();
-                ExprVector eab = ea.Minus(eb);
-
-                // Construct a vector from the point to either endpoint of
-                // the line segment, and choose the longer of these.
-                ExprVector eap = ea.Minus(ep);
-                ExprVector ebp = eb.Minus(ep);
-                ExprVector elp =
-                    (ebp.Magnitude()->Eval() > eap.Magnitude()->Eval()) ?
-                        ebp : eap;
-
-                if(p->group.v == group.v) {
-                    AddEq(l, VectorsParallel(0, eab, elp), 0);
-                    AddEq(l, VectorsParallel(1, eab, elp), 1);
-                } else {
-                    AddEq(l, VectorsParallel(0, elp, eab), 0);
-                    AddEq(l, VectorsParallel(1, elp, eab), 1);
-                }
-            } else {
-                AddEq(l, PointLineDistance(workplane, ptA, entityA), 0);
-            }
+        case Type::PT_ON_LINE: {
+            EntityBase *ln = SK.GetEntity(entityA);
+            EntityBase *a = SK.GetEntity(ln->point[0]);
+            EntityBase *b = SK.GetEntity(ln->point[1]);
+            EntityBase *p = SK.GetEntity(ptA);
+
+            ExprVector ep = p->PointGetExprsInWorkplane(workplane);
+            ExprVector ea = a->PointGetExprsInWorkplane(workplane);
+            ExprVector eb = b->PointGetExprsInWorkplane(workplane);
+
+            ExprVector ptOnLine = ea.Plus(eb.Minus(ea).ScaledBy(Expr::From(h.param(0))));
+            ExprVector eq = ptOnLine.Minus(ep);
+
+            AddEq(l, eq.x, 0);
+            AddEq(l, eq.y, 1);
+            if(workplane.v == EntityBase::FREE_IN_3D.v) AddEq(l, eq.z, 2);
             return;
+        }
 
         case Type::PT_ON_CIRCLE: {
             // This actually constrains the point to lie on the cylinder.

src/file.cpp

@@ -556,12 +556,46 @@ void SolveSpaceUI::UpgradeLegacyData() {
                 }
                 entity.Clear();
                 param.Clear();
+                break;
             }
 
             default:
                 break;
         }
     }
+
+    IdList<Param,hParam> oldParam = {};
+    SK.param.DeepCopyInto(&oldParam);
+    SS.GenerateAll(SS.Generate::REGEN);
+    for(Constraint &c : SK.constraint) {
+        switch(c.type) {
+            case Constraint::Type::PT_ON_LINE: {
+                IdList<Param,hParam> param = {};
+                c.Generate(&param);
+                bool allParamsExist = true;
+                for(Param &p : param) {
+                    if(oldParam.FindByIdNoOops(p.h) != NULL) continue;
+                    allParamsExist = false;
+                }
+                if(!allParamsExist) {
+                    EntityBase *ln = SK.GetEntity(c.entityA);
+                    EntityBase *a = SK.GetEntity(ln->point[0]);
+                    EntityBase *b = SK.GetEntity(ln->point[1]);
+                    EntityBase *p = SK.GetEntity(c.ptA);
+
+                    ExprVector ep = p->PointGetExprsInWorkplane(c.workplane);
+                    ExprVector ea = a->PointGetExprsInWorkplane(c.workplane);
+                    ExprVector eb = b->PointGetExprsInWorkplane(c.workplane);
+                    ExprVector eba = eb.Minus(ea);
+                    Param *param = SK.GetParam(c.h.param(0));
+                    param->val = eba.Dot(ep.Minus(ea))->Eval() / eba.Dot(eba)->Eval();
+                }
+                break;
+            }
+            default:
+                break;
+        }
+    }
 }
 
 bool SolveSpaceUI::LoadEntitiesFromFile(const std::string &filename, EntityList *le,

src/generate.cpp

@@ -239,6 +239,12 @@ void SolveSpaceUI::GenerateAll(Generate type, bool andFindFree, bool genForBBox)
 
             r->Generate(&(SK.entity), &(SK.param));
         }
+        for(j = 0; j < SK.constraint.n; j++) {
+            Constraint *c = &SK.constraint.elem[j];
+            if(c->group.v != g->h.v) continue;
+
+            c->Generate(&(SK.param));
+        }
         g->Generate(&(SK.entity), &(SK.param));
 
         // The requests and constraints depend on stuff in this or the
@@ -509,6 +515,12 @@ void SolveSpaceUI::SolveGroup(hGroup hg, bool andFindFree) {
 
         r->Generate(&(sys.entity), &(sys.param));
     }
+    for(i = 0; i < SK.constraint.n; i++) {
+        Constraint *c = &SK.constraint.elem[i];
+        if(c->group.v != hg.v) continue;
+
+        c->Generate(&(sys.param));
+    }
     // And for the group itself
     Group *g = SK.GetGroup(hg);
     g->Generate(&(sys.entity), &(sys.param));

src/sketch.h

@@ -569,6 +569,7 @@ class hConstraint {
     uint32_t v;
 
     inline hEquation equation(int i) const;
+    inline hParam param(int i) const;
 };
 
 class ConstraintBase {
@@ -638,8 +639,10 @@ class ConstraintBase {
 
     bool HasLabel() const;
 
-    void Generate(IdList<Equation,hEquation> *l) const;
-    void GenerateReal(IdList<Equation,hEquation> *l) const;
+    void Generate(IdList<Param, hParam> *param) const;
+
+    void GenerateEquations(IdList<Equation,hEquation> *entity,
+                           bool forReference = false) const;
     // Some helpers when generating symbolic constraint equations
     void ModifyToSatisfy();
     void AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index) const;
@@ -874,6 +877,8 @@ inline hRequest hParam::request() const
 
 inline hEquation hConstraint::equation(int i) const
     { hEquation r; r.v = (v << 16) | (uint32_t)i; return r; }
+inline hParam hConstraint::param(int i) const
+    { hParam r; r.v = v | 0x40000000 | (uint32_t)i; return r; }
 
 inline bool hEquation::isFromConstraint() const
     { if(v & 0xc0000000) return false; else return true; }

src/system.cpp

@@ -346,7 +346,7 @@ void System::WriteEquationsExceptFor(hConstraint hc, Group *g) {
             continue;
         }
 
-        c->Generate(&eq);
+        c->GenerateEquations(&eq);
     }
     // And the equations from entities
     for(i = 0; i < SK.entity.n; i++) {