Rewrite equations 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 free parameter. The point being coincident with
the start of the line corresponds to the parameter being zero, and
point being coincident with the end corresponds to one).

In effect, instead of constraining two of three degrees of freedom
(for which the equations do not exist because of the hairy ball
theorem), it constrains three and adds one more.

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 {
+void ConstraintBase::Generate(IdList<Param,hParam> *l) const {
-    if(!reference) {
+    switch(type) {
-        GenerateReal(l);
+        case Type::PT_ON_LINE: {
+            Param p = {};
+            p.h = h.param(0);
+            l->Add(&p);
+            break;
         }
-}
-void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) const {
-    Expr *exA = Expr::From(valA);
 
+        default:
+            break;
+    }
+}
+
+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:
+        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 ep = p->PointGetExprsInWorkplane(workplane);
-                ExprVector ea = a->PointGetExprs();
+            ExprVector ea = a->PointGetExprsInWorkplane(workplane);
-                ExprVector eb = b->PointGetExprs();
+            ExprVector eb = b->PointGetExprsInWorkplane(workplane);
-                ExprVector eab = ea.Minus(eb);
 
-                // Construct a vector from the point to either endpoint of
+            ExprVector ptOnLine = ea.Plus(eb.Minus(ea).ScaledBy(Expr::From(h.param(0))));
-                // the line segment, and choose the longer of these.
+            ExprVector eq = ptOnLine.Minus(ep);
-                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, eq.x, 0);
-                    AddEq(l, VectorsParallel(0, eab, elp), 0);
+            AddEq(l, eq.y, 1);
-                    AddEq(l, VectorsParallel(1, eab, elp), 1);
+            if(workplane.v == EntityBase::FREE_IN_3D.v) AddEq(l, eq.z, 2);
-                } else {
-                    AddEq(l, VectorsParallel(0, elp, eab), 0);
-                    AddEq(l, VectorsParallel(1, elp, eab), 1);
-                }
-            } else {
-                AddEq(l, PointLineDistance(workplane, ptA, entityA), 0);
-            }
             return;
+        }
 
         case Type::PT_ON_CIRCLE: {
             // This actually constrains the point to lie on the cylinder.

src/file.cpp

@@ -556,12 +556,49 @@ void SolveSpaceUI::UpgradeLegacyData() {
                 }
                 entity.Clear();
                 param.Clear();
+                break;
             }
 
             default:
                 break;
         }
     }
+
+    IdList<Param,hParam> oldParam = {};
+    SK.param.DeepCopyInto(&oldParam);
+    SS.GenerateAll(SolveSpaceUI::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;
+                }
+                param.Clear();
+
+                if(!allParamsExist) {
+                    EntityBase *eln = SK.GetEntity(c.entityA);
+                    EntityBase *ea = SK.GetEntity(eln->point[0]);
+                    EntityBase *eb = SK.GetEntity(eln->point[1]);
+                    EntityBase *ep = SK.GetEntity(c.ptA);
+
+                    ExprVector exp = ep->PointGetExprsInWorkplane(c.workplane);
+                    ExprVector exa = ea->PointGetExprsInWorkplane(c.workplane);
+                    ExprVector exb = eb->PointGetExprsInWorkplane(c.workplane);
+                    ExprVector exba = exb.Minus(exa);
+                    Param *p = SK.GetParam(c.h.param(0));
+                    p->val = exba.Dot(exp.Minus(exa))->Eval() / exba.Dot(exba)->Eval();
+                }
+                break;
+            }
+            default:
+                break;
+        }
+    }
+    oldParam.Clear();
 }
 
 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 @@ public:
     uint32_t v;
 
     inline hEquation equation(int i) const;
+    inline hParam param(int i) const;
 };
 
 class ConstraintBase {
@@ -638,8 +639,10 @@ public:
 
     bool HasLabel() const;
 
-    void Generate(IdList<Equation,hEquation> *l) const;
+    void Generate(IdList<Param, hParam> *param) const;
-    void GenerateReal(IdList<Equation,hEquation> *l) 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++) {

test/constraint/pt_on_line/left_free_in_3d.slvs

@@ -147,6 +147,10 @@ AddParam
 Param.h.v.=00050012
 AddParam
 
+Param.h.v.=40000001
+Param.val=0.33333333333333331483
+AddParam
+
 Request.h.v=00000001
 Request.type=100
 Request.group.v=00000001

test/constraint/pt_on_line/normal.slvs

@@ -138,6 +138,10 @@ Param.h.v.=00050011
 Param.val=5.00000000000000000000
 AddParam
 
+Param.h.v.=40000001
+Param.val=0.33333333333333331483
+AddParam
+
 Request.h.v=00000001
 Request.type=100
 Request.group.v=00000001

test/constraint/pt_on_line/right_free_in_3d.slvs

@@ -147,6 +147,10 @@ AddParam
 Param.h.v.=00050012
 AddParam
 
+Param.h.v.=40000001
+Param.val=0.66666666666666662966
+AddParam
+
 Request.h.v=00000001
 Request.type=100
 Request.group.v=00000001