Solvespace

1
//-----------------------------------------------------------------------------
2
// Given a constraint, generate one or more equations in our symbolic algebra
3
// system to represent that constraint; also various geometric helper
4
// functions for that.
5
//
6
// Copyright 2008-2013 Jonathan Westhues.
7
//-----------------------------------------------------------------------------
8
#include "solvespace.h"
9

10
const hConstraint ConstraintBase::NO_CONSTRAINT = { 0 };
11

12
bool ConstraintBase::HasLabel() const {
13
    switch(type) {
14
        case Type::PT_LINE_DISTANCE:
15
        case Type::PT_PLANE_DISTANCE:
16
        case Type::PT_FACE_DISTANCE:
17
        case Type::PT_PT_DISTANCE:
18
        case Type::PROJ_PT_DISTANCE:
19
        case Type::DIAMETER:
20
        case Type::LENGTH_RATIO:
21
        case Type::ARC_ARC_LEN_RATIO:  
22
        case Type::ARC_LINE_LEN_RATIO: 
23
        case Type::LENGTH_DIFFERENCE:
24
        case Type::ARC_ARC_DIFFERENCE: 
25
        case Type::ARC_LINE_DIFFERENCE:
26
        case Type::ANGLE:
27
        case Type::COMMENT:
28
            return true;
29

30
        default:
31
            return false;
32
    }
33
}
34

35
bool ConstraintBase::IsProjectible() const {
36
    switch(type) {
37
        case Type::POINTS_COINCIDENT:
38
        case Type::PT_PT_DISTANCE:
39
        case Type::PT_LINE_DISTANCE:
40
        case Type::PT_ON_LINE:
41
        case Type::EQUAL_LENGTH_LINES:
42
        case Type::EQ_LEN_PT_LINE_D:
43
        case Type::EQ_PT_LN_DISTANCES:
44
        case Type::EQUAL_ANGLE:
45
        case Type::LENGTH_RATIO:
46
        case Type::ARC_ARC_LEN_RATIO:  
47
        case Type::ARC_LINE_LEN_RATIO: 
48
        case Type::LENGTH_DIFFERENCE:
49
        case Type::ARC_ARC_DIFFERENCE: 
50
        case Type::ARC_LINE_DIFFERENCE:
51
        case Type::SYMMETRIC:
52
        case Type::SYMMETRIC_HORIZ:
53
        case Type::SYMMETRIC_VERT:
54
        case Type::SYMMETRIC_LINE:
55
        case Type::AT_MIDPOINT:
56
        case Type::HORIZONTAL:
57
        case Type::VERTICAL:
58
        case Type::ANGLE:
59
        case Type::PARALLEL:
60
        case Type::PERPENDICULAR:
61
        case Type::WHERE_DRAGGED:
62
        case Type::COMMENT:
63
            return true;
64

65
        case Type::PT_PLANE_DISTANCE:
66
        case Type::PT_FACE_DISTANCE:
67
        case Type::PROJ_PT_DISTANCE:
68
        case Type::PT_IN_PLANE:
69
        case Type::PT_ON_FACE:
70
        case Type::EQUAL_LINE_ARC_LEN:
71
        case Type::DIAMETER:
72
        case Type::PT_ON_CIRCLE:
73
        case Type::SAME_ORIENTATION:
74
        case Type::CUBIC_LINE_TANGENT:
75
        case Type::CURVE_CURVE_TANGENT:
76
        case Type::ARC_LINE_TANGENT:
77
        case Type::EQUAL_RADIUS:
78
            return false;
79
    }
80
    ssassert(false, "Impossible");
81
}
82

83
ExprVector ConstraintBase::VectorsParallel3d(ExprVector a, ExprVector b, hParam p) {
84
    return a.Minus(b.ScaledBy(Expr::From(p)));
85
}
86

87
Expr *ConstraintBase::PointLineDistance(hEntity wrkpl, hEntity hpt, hEntity hln)
88
{
89
    EntityBase *ln = SK.GetEntity(hln);
90
    EntityBase *a = SK.GetEntity(ln->point[0]);
91
    EntityBase *b = SK.GetEntity(ln->point[1]);
92

93
    EntityBase *p = SK.GetEntity(hpt);
94

95
    if(wrkpl == EntityBase::FREE_IN_3D) {
96
        ExprVector ep = p->PointGetExprs();
97

98
        ExprVector ea = a->PointGetExprs();
99
        ExprVector eb = b->PointGetExprs();
100
        ExprVector eab = ea.Minus(eb);
101
        Expr *m = eab.Magnitude();
102

103
        return ((eab.Cross(ea.Minus(ep))).Magnitude())->Div(m);
104
    } else {
105
        Expr *ua, *va, *ub, *vb;
106
        a->PointGetExprsInWorkplane(wrkpl, &ua, &va);
107
        b->PointGetExprsInWorkplane(wrkpl, &ub, &vb);
108

109
        Expr *du = ua->Minus(ub);
110
        Expr *dv = va->Minus(vb);
111

112
        Expr *u, *v;
113
        p->PointGetExprsInWorkplane(wrkpl, &u, &v);
114

115
        Expr *m = ((du->Square())->Plus(dv->Square()))->Sqrt();
116

117
        Expr *proj = (dv->Times(ua->Minus(u)))->Minus(
118
                     (du->Times(va->Minus(v))));
119

120
        return proj->Div(m);
121
    }
122
}
123

124
Expr *ConstraintBase::PointPlaneDistance(ExprVector p, hEntity hpl) {
125
    ExprVector n;
126
    Expr *d;
127
    SK.GetEntity(hpl)->WorkplaneGetPlaneExprs(&n, &d);
128
    return (p.Dot(n))->Minus(d);
129
}
130

131
Expr *ConstraintBase::Distance(hEntity wrkpl, hEntity hpa, hEntity hpb) {
132
    EntityBase *pa = SK.GetEntity(hpa);
133
    EntityBase *pb = SK.GetEntity(hpb);
134
    ssassert(pa->IsPoint() && pb->IsPoint(),
135
             "Expected two points to measure projected distance between");
136

137
    if(wrkpl == EntityBase::FREE_IN_3D) {
138
        // This is true distance
139
        ExprVector ea, eb, eab;
140
        ea = pa->PointGetExprs();
141
        eb = pb->PointGetExprs();
142
        eab = ea.Minus(eb);
143

144
        return eab.Magnitude();
145
    } else {
146
        // This is projected distance, in the given workplane.
147
        Expr *au, *av, *bu, *bv;
148

149
        pa->PointGetExprsInWorkplane(wrkpl, &au, &av);
150
        pb->PointGetExprsInWorkplane(wrkpl, &bu, &bv);
151

152
        Expr *du = au->Minus(bu);
153
        Expr *dv = av->Minus(bv);
154

155
        return ((du->Square())->Plus(dv->Square()))->Sqrt();
156
    }
157
}
158

159
//-----------------------------------------------------------------------------
160
// Return the cosine of the angle between two vectors. If a workplane is
161
// specified, then it's the cosine of their projections into that workplane.
162
//-----------------------------------------------------------------------------
163
Expr *ConstraintBase::DirectionCosine(hEntity wrkpl,
164
                                      ExprVector ae, ExprVector be)
165
{
166
    if(wrkpl == EntityBase::FREE_IN_3D) {
167
        Expr *mags = (ae.Magnitude())->Times(be.Magnitude());
168
        return (ae.Dot(be))->Div(mags);
169
    } else {
170
        EntityBase *w = SK.GetEntity(wrkpl);
171
        ExprVector u = w->Normal()->NormalExprsU();
172
        ExprVector v = w->Normal()->NormalExprsV();
173
        Expr *ua = u.Dot(ae);
174
        Expr *va = v.Dot(ae);
175
        Expr *ub = u.Dot(be);
176
        Expr *vb = v.Dot(be);
177
        Expr *maga = (ua->Square()->Plus(va->Square()))->Sqrt();
178
        Expr *magb = (ub->Square()->Plus(vb->Square()))->Sqrt();
179
        Expr *dot = (ua->Times(ub))->Plus(va->Times(vb));
180
        return dot->Div(maga->Times(magb));
181
    }
182
}
183

184
ExprVector ConstraintBase::PointInThreeSpace(hEntity workplane,
185
                                             Expr *u, Expr *v)
186
{
187
    EntityBase *w = SK.GetEntity(workplane);
188

189
    ExprVector ub = w->Normal()->NormalExprsU();
190
    ExprVector vb = w->Normal()->NormalExprsV();
191
    ExprVector ob = w->WorkplaneGetOffsetExprs();
192

193
    return (ub.ScaledBy(u)).Plus(vb.ScaledBy(v)).Plus(ob);
194
}
195

196
void ConstraintBase::ModifyToSatisfy() {
197
    if(type == Type::ANGLE) {
198
        Vector a = SK.GetEntity(entityA)->VectorGetNum();
199
        Vector b = SK.GetEntity(entityB)->VectorGetNum();
200
        if(other) a = a.ScaledBy(-1);
201
        if(workplane != EntityBase::FREE_IN_3D) {
202
            a = a.ProjectVectorInto(workplane);
203
            b = b.ProjectVectorInto(workplane);
204
        }
205
        double c = (a.Dot(b))/(a.Magnitude() * b.Magnitude());
206
        valA = acos(c)*180/PI;
207
    } else if(type == Type::PT_ON_LINE) {
208
        EntityBase *eln = SK.GetEntity(entityA);
209
        EntityBase *ea = SK.GetEntity(eln->point[0]);
210
        EntityBase *eb = SK.GetEntity(eln->point[1]);
211
        EntityBase *ep = SK.GetEntity(ptA);
212
        ExprVector exp = ep->PointGetExprsInWorkplane(workplane);
213
        ExprVector exa = ea->PointGetExprsInWorkplane(workplane);
214
        ExprVector exb = eb->PointGetExprsInWorkplane(workplane);
215
        ExprVector exba = exb.Minus(exa);
216
        SK.GetParam(valP)->val = exba.Dot(exp.Minus(exa))->Eval() / exba.Dot(exba)->Eval();
217
    } else {
218
        // We'll fix these ones up by looking at their symbolic equation;
219
        // that means no extra work.
220
        IdList<Equation,hEquation> l = {};
221
        // Generate the equations even if this is a reference dimension
222
        GenerateEquations(&l, /*forReference=*/true);
223
        ssassert(l.n == 1, "Expected constraint to generate a single equation");
224

225
        // These equations are written in the form f(...) - d = 0, where
226
        // d is the value of the valA.
227
        valA += (l[0].e)->Eval();
228

229
        l.Clear();
230
    }
231
}
232

233
void ConstraintBase::AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index) const
234
{
235
    Equation eq;
236
    eq.e = expr;
237
    eq.h = h.equation(index);
238
    l->Add(&eq);
239
}
240

241
void ConstraintBase::AddEq(IdList<Equation,hEquation> *l, const ExprVector &v,
242
                           int baseIndex) const {
243
    AddEq(l, v.x, baseIndex);
244
    AddEq(l, v.y, baseIndex + 1);
245
    if(workplane == EntityBase::FREE_IN_3D) {
246
        AddEq(l, v.z, baseIndex + 2);
247
    }
248
}
249

250
void ConstraintBase::Generate(IdList<Param,hParam> *l) {
251
    switch(type) {
252
        case Type::PARALLEL:
253
        case Type::CUBIC_LINE_TANGENT:
254
            // Add new parameter only when we operate in 3d space
255
            if(workplane != EntityBase::FREE_IN_3D) break;
256
            // fallthrough
257
        case Type::SAME_ORIENTATION:
258
        case Type::PT_ON_LINE: {
259
            Param p = {};
260
            valP = h.param(0);
261
            p.h = valP;
262
            l->Add(&p);
263
            break;
264
        }
265

266
        default:
267
            break;
268
    }
269
}
270

271
void ConstraintBase::GenerateEquations(IdList<Equation,hEquation> *l,
272
                                       bool forReference) const {
273
    if(reference && !forReference) return;
274

275
    Expr *exA = Expr::From(valA);
276
    switch(type) {
277
        case Type::PT_PT_DISTANCE:
278
            AddEq(l, Distance(workplane, ptA, ptB)->Minus(exA), 0);
279
            return;
280

281
        case Type::PROJ_PT_DISTANCE: {
282
            ExprVector pA = SK.GetEntity(ptA)->PointGetExprs(),
283
                       pB = SK.GetEntity(ptB)->PointGetExprs(),
284
                       dp = pB.Minus(pA);
285

286
            ExprVector pp = SK.GetEntity(entityA)->VectorGetExprs();
287
            pp = pp.WithMagnitude(Expr::From(1.0));
288

289
            AddEq(l, (dp.Dot(pp))->Minus(exA), 0);
290
            return;
291
        }
292

293
        case Type::PT_LINE_DISTANCE:
294
            AddEq(l,
295
                PointLineDistance(workplane, ptA, entityA)->Minus(exA), 0);
296
            return;
297

298
        case Type::PT_PLANE_DISTANCE: {
299
            ExprVector pt = SK.GetEntity(ptA)->PointGetExprs();
300
            AddEq(l, (PointPlaneDistance(pt, entityA))->Minus(exA), 0);
301
            return;
302
        }
303

304
        case Type::PT_FACE_DISTANCE: {
305
            ExprVector pt = SK.GetEntity(ptA)->PointGetExprs();
306
            EntityBase *f = SK.GetEntity(entityA);
307
            ExprVector p0 = f->FaceGetPointExprs();
308
            ExprVector n = f->FaceGetNormalExprs();
309
            AddEq(l, (pt.Minus(p0)).Dot(n)->Minus(exA), 0);
310
            return;
311
        }
312

313
        case Type::EQUAL_LENGTH_LINES: {
314
            EntityBase *a = SK.GetEntity(entityA);
315
            EntityBase *b = SK.GetEntity(entityB);
316
            AddEq(l, Distance(workplane, a->point[0], a->point[1])->Minus(
317
                     Distance(workplane, b->point[0], b->point[1])), 0);
318
            return;
319
        }
320

321
        // These work on distance squared, since the pt-line distances are
322
        // signed, and we want the absolute value.
323
        case Type::EQ_LEN_PT_LINE_D: {
324
            EntityBase *forLen = SK.GetEntity(entityA);
325
            Expr *d1 = Distance(workplane, forLen->point[0], forLen->point[1]);
326
            Expr *d2 = PointLineDistance(workplane, ptA, entityB);
327
            AddEq(l, (d1->Square())->Minus(d2->Square()), 0);
328
            return;
329
        }
330
        case Type::EQ_PT_LN_DISTANCES: {
331
            Expr *d1 = PointLineDistance(workplane, ptA, entityA);
332
            Expr *d2 = PointLineDistance(workplane, ptB, entityB);
333
            AddEq(l, (d1->Square())->Minus(d2->Square()), 0);
334
            return;
335
        }
336

337
        case Type::LENGTH_RATIO: {
338
            EntityBase *a = SK.GetEntity(entityA);
339
            EntityBase *b = SK.GetEntity(entityB);
340
            Expr *la = Distance(workplane, a->point[0], a->point[1]);
341
            Expr *lb = Distance(workplane, b->point[0], b->point[1]);
342
            AddEq(l, (la->Div(lb))->Minus(exA), 0);
343
            return;
344
        }
345
        
346
        case Type::ARC_ARC_LEN_RATIO: {
347
            EntityBase *arc1  = SK.GetEntity(entityA),
348
                       *arc2  = SK.GetEntity(entityB);
349

350
            // And get the arc1 radius, and the cosine of its angle
351
            EntityBase *ao1 = SK.GetEntity(arc1->point[0]),
352
                       *as1 = SK.GetEntity(arc1->point[1]),
353
                       *af1 = SK.GetEntity(arc1->point[2]);
354

355
            ExprVector aos1 = (as1->PointGetExprs()).Minus(ao1->PointGetExprs()),
356
                       aof1 = (af1->PointGetExprs()).Minus(ao1->PointGetExprs());
357
            Expr *r1 = aof1.Magnitude();
358

359
            ExprVector n1 = arc1->Normal()->NormalExprsN();
360
            ExprVector u1 = aos1.WithMagnitude(Expr::From(1.0));
361
            ExprVector v1 = n1.Cross(u1);
362
            // so in our new csys, we start at (1, 0, 0)
363
            Expr *costheta1 = aof1.Dot(u1)->Div(r1);
364
            Expr *sintheta1 = aof1.Dot(v1)->Div(r1);
365

366
            double thetas1, thetaf1, dtheta1;
367
            arc1->ArcGetAngles(&thetas1, &thetaf1, &dtheta1);
368
            Expr *theta1;
369
            if(dtheta1 < 3*PI/4) {
370
                theta1 = costheta1->ACos();
371
            } else if(dtheta1 < 5*PI/4) {
372
                // As the angle crosses pi, cos theta1 is not invertible;
373
                // so use the sine to stop blowing up
374
                theta1 = Expr::From(PI)->Minus(sintheta1->ASin());
375
            } else {
376
                theta1 = (Expr::From(2*PI))->Minus(costheta1->ACos());
377
            }
378
            
379
            // And get the arc2 radius, and the cosine of its angle
380
            EntityBase *ao2 = SK.GetEntity(arc2->point[0]),
381
                       *as2 = SK.GetEntity(arc2->point[1]),
382
                       *af2 = SK.GetEntity(arc2->point[2]);
383

384
            ExprVector aos2 = (as2->PointGetExprs()).Minus(ao2->PointGetExprs()),
385
                       aof2 = (af2->PointGetExprs()).Minus(ao2->PointGetExprs());
386
            Expr *r2 = aof2.Magnitude();
387

388
            ExprVector n2 = arc2->Normal()->NormalExprsN();
389
            ExprVector u2 = aos2.WithMagnitude(Expr::From(1.0));
390
            ExprVector v2 = n2.Cross(u2);
391
            // so in our new csys, we start at (1, 0, 0)
392
            Expr *costheta2 = aof2.Dot(u2)->Div(r2);
393
            Expr *sintheta2 = aof2.Dot(v2)->Div(r2);
394

395
            double thetas2, thetaf2, dtheta2;
396
            arc2->ArcGetAngles(&thetas2, &thetaf2, &dtheta2);
397
            Expr *theta2;
398
            if(dtheta2 < 3*PI/4) {
399
                theta2 = costheta2->ACos();
400
            } else if(dtheta2 < 5*PI/4) {
401
                // As the angle crosses pi, cos theta2 is not invertible;
402
                // so use the sine to stop blowing up
403
                theta2 = Expr::From(PI)->Minus(sintheta2->ASin());
404
            } else {
405
                theta2 = (Expr::From(2*PI))->Minus(costheta2->ACos());
406
            }
407
            // And write the equation; (r1*theta1) / ( r2*theta2) = some ratio
408
            AddEq(l, (r1->Times(theta1))->Div(r2->Times(theta2))->Minus(exA), 0);
409
            return;
410
        }
411
        
412
        case Type::ARC_LINE_LEN_RATIO: {
413
            EntityBase *line   = SK.GetEntity(entityA),
414
                       *arc1   = SK.GetEntity(entityB);
415
             
416
            Expr *ll = Distance(workplane, line->point[0], line->point[1]);
417
               
418
            // And get the arc1 radius, and the cosine of its angle
419
            EntityBase *ao1 = SK.GetEntity(arc1->point[0]),
420
                       *as1 = SK.GetEntity(arc1->point[1]),
421
                       *af1 = SK.GetEntity(arc1->point[2]);
422

423
            ExprVector aos1 = (as1->PointGetExprs()).Minus(ao1->PointGetExprs()),
424
                       aof1 = (af1->PointGetExprs()).Minus(ao1->PointGetExprs());
425
            Expr *r1 = aof1.Magnitude();
426
            ExprVector n1 = arc1->Normal()->NormalExprsN();
427
            ExprVector u1 = aos1.WithMagnitude(Expr::From(1.0));
428
            ExprVector v1 = n1.Cross(u1);
429
            // so in our new csys, we start at (1, 0, 0)
430
            Expr *costheta1 = aof1.Dot(u1)->Div(r1);
431
            Expr *sintheta1 = aof1.Dot(v1)->Div(r1);
432

433
            double thetas1, thetaf1, dtheta1;
434
            arc1->ArcGetAngles(&thetas1, &thetaf1, &dtheta1);
435
            Expr *theta1;
436
            if(dtheta1 < 3*PI/4) {
437
                theta1 = costheta1->ACos();
438
            } else if(dtheta1 < 5*PI/4) {
439
                // As the angle crosses pi, cos theta1 is not invertible;
440
                // so use the sine to stop blowing up
441
                theta1 = Expr::From(PI)->Minus(sintheta1->ASin());
442
            } else {
443
                theta1 = (Expr::From(2*PI))->Minus(costheta1->ACos());
444
            }
445
            // And write the equation; (r1*theta1) / ( length) = some ratio
446
            AddEq(l, (r1->Times(theta1))->Div(ll)->Minus(exA), 0);
447
            return;
448
        }
449

450
        case Type::LENGTH_DIFFERENCE: {
451
            EntityBase *a = SK.GetEntity(entityA);
452
            EntityBase *b = SK.GetEntity(entityB);
453
            Expr *la = Distance(workplane, a->point[0], a->point[1]);
454
            Expr *lb = Distance(workplane, b->point[0], b->point[1]);
455
            AddEq(l, (la->Minus(lb))->Minus(exA), 0);
456
            return;
457
        }
458
        
459
        case Type::ARC_ARC_DIFFERENCE: {
460
            EntityBase *arc1  = SK.GetEntity(entityA),
461
                       *arc2  = SK.GetEntity(entityB);
462

463
            // And get the arc1 radius, and the cosine of its angle
464
            EntityBase *ao1 = SK.GetEntity(arc1->point[0]),
465
                       *as1 = SK.GetEntity(arc1->point[1]),
466
                       *af1 = SK.GetEntity(arc1->point[2]);
467

468
            ExprVector aos1 = (as1->PointGetExprs()).Minus(ao1->PointGetExprs()),
469
                       aof1 = (af1->PointGetExprs()).Minus(ao1->PointGetExprs());
470
            Expr *r1 = aof1.Magnitude();
471

472
            ExprVector n1 = arc1->Normal()->NormalExprsN();
473
            ExprVector u1 = aos1.WithMagnitude(Expr::From(1.0));
474
            ExprVector v1 = n1.Cross(u1);
475
            // so in our new csys, we start at (1, 0, 0)
476
            Expr *costheta1 = aof1.Dot(u1)->Div(r1);
477
            Expr *sintheta1 = aof1.Dot(v1)->Div(r1);
478

479
            double thetas1, thetaf1, dtheta1;
480
            arc1->ArcGetAngles(&thetas1, &thetaf1, &dtheta1);
481
            Expr *theta1;
482
            if(dtheta1 < 3*PI/4) {
483
                theta1 = costheta1->ACos();
484
            } else if(dtheta1 < 5*PI/4) {
485
                // As the angle crosses pi, cos theta1 is not invertible;
486
                // so use the sine to stop blowing up
487
                theta1 = Expr::From(PI)->Minus(sintheta1->ASin());
488
            } else {
489
                theta1 = (Expr::From(2*PI))->Minus(costheta1->ACos());
490
            }
491
            
492
            // And get the arc2 radius, and the cosine of its angle
493
            EntityBase *ao2 = SK.GetEntity(arc2->point[0]),
494
                       *as2 = SK.GetEntity(arc2->point[1]),
495
                       *af2 = SK.GetEntity(arc2->point[2]);
496

497
            ExprVector aos2 = (as2->PointGetExprs()).Minus(ao2->PointGetExprs()),
498
                       aof2 = (af2->PointGetExprs()).Minus(ao2->PointGetExprs());
499
            Expr *r2 = aof2.Magnitude();
500

501
            ExprVector n2 = arc2->Normal()->NormalExprsN();
502
            ExprVector u2 = aos2.WithMagnitude(Expr::From(1.0));
503
            ExprVector v2 = n2.Cross(u2);
504
            // so in our new csys, we start at (1, 0, 0)
505
            Expr *costheta2 = aof2.Dot(u2)->Div(r2);
506
            Expr *sintheta2 = aof2.Dot(v2)->Div(r2);
507

508
            double thetas2, thetaf2, dtheta2;
509
            arc2->ArcGetAngles(&thetas2, &thetaf2, &dtheta2);
510
            Expr *theta2;
511
            if(dtheta2 < 3*PI/4) {
512
                theta2 = costheta2->ACos();
513
            } else if(dtheta2 < 5*PI/4) {
514
                // As the angle crosses pi, cos theta2 is not invertible;
515
                // so use the sine to stop blowing up
516
                theta2 = Expr::From(PI)->Minus(sintheta2->ASin());
517
            } else {
518
                theta2 = (Expr::From(2*PI))->Minus(costheta2->ACos());
519
            }
520
            // And write the equation; (r1*theta1) - ( r2*theta2) = some difference
521
            AddEq(l, (r1->Times(theta1))->Minus(r2->Times(theta2))->Minus(exA), 0);
522
            return;
523
        }
524
        
525
        case Type::ARC_LINE_DIFFERENCE: {
526
            EntityBase *line   = SK.GetEntity(entityA),
527
                       *arc1   = SK.GetEntity(entityB);
528
             
529
            Expr *ll = Distance(workplane, line->point[0], line->point[1]);
530
               
531
            // And get the arc1 radius, and the cosine of its angle
532
            EntityBase *ao1 = SK.GetEntity(arc1->point[0]),
533
                       *as1 = SK.GetEntity(arc1->point[1]),
534
                       *af1 = SK.GetEntity(arc1->point[2]);
535

536
            ExprVector aos1 = (as1->PointGetExprs()).Minus(ao1->PointGetExprs()),
537
                       aof1 = (af1->PointGetExprs()).Minus(ao1->PointGetExprs());
538
            Expr *r1 = aof1.Magnitude();
539
            ExprVector n1 = arc1->Normal()->NormalExprsN();
540
            ExprVector u1 = aos1.WithMagnitude(Expr::From(1.0));
541
            ExprVector v1 = n1.Cross(u1);
542
            // so in our new csys, we start at (1, 0, 0)
543
            Expr *costheta1 = aof1.Dot(u1)->Div(r1);
544
            Expr *sintheta1 = aof1.Dot(v1)->Div(r1);
545

546
            double thetas1, thetaf1, dtheta1;
547
            arc1->ArcGetAngles(&thetas1, &thetaf1, &dtheta1);
548
            Expr *theta1;
549
            if(dtheta1 < 3*PI/4) {
550
                theta1 = costheta1->ACos();
551
            } else if(dtheta1 < 5*PI/4) {
552
                // As the angle crosses pi, cos theta1 is not invertible;
553
                // so use the sine to stop blowing up
554
                theta1 = Expr::From(PI)->Minus(sintheta1->ASin());
555
            } else {
556
                theta1 = (Expr::From(2*PI))->Minus(costheta1->ACos());
557
            }
558
            // And write the equation; (r1*theta1) - ( length) = some difference
559
            AddEq(l, (r1->Times(theta1))->Minus(ll)->Minus(exA), 0);
560
            return;
561
        }
562
        
563
        case Type::DIAMETER: {
564
            EntityBase *circle = SK.GetEntity(entityA);
565
            Expr *r = circle->CircleGetRadiusExpr();
566
            AddEq(l, (r->Times(Expr::From(2)))->Minus(exA), 0);
567
            return;
568
        }
569

570
        case Type::EQUAL_RADIUS: {
571
            EntityBase *c1 = SK.GetEntity(entityA);
572
            EntityBase *c2 = SK.GetEntity(entityB);
573
            AddEq(l, (c1->CircleGetRadiusExpr())->Minus(
574
                      c2->CircleGetRadiusExpr()), 0);
575
            return;
576
        }
577

578
        case Type::EQUAL_LINE_ARC_LEN: {
579
            EntityBase *line = SK.GetEntity(entityA),
580
                       *arc  = SK.GetEntity(entityB);
581

582
            // Get the line length
583
            ExprVector l0 = SK.GetEntity(line->point[0])->PointGetExprs(),
584
                       l1 = SK.GetEntity(line->point[1])->PointGetExprs();
585
            Expr *ll = (l1.Minus(l0)).Magnitude();
586

587
            // And get the arc radius, and the cosine of its angle
588
            EntityBase *ao = SK.GetEntity(arc->point[0]),
589
                       *as = SK.GetEntity(arc->point[1]),
590
                       *af = SK.GetEntity(arc->point[2]);
591

592
            ExprVector aos = (as->PointGetExprs()).Minus(ao->PointGetExprs()),
593
                       aof = (af->PointGetExprs()).Minus(ao->PointGetExprs());
594
            Expr *r = aof.Magnitude();
595

596
            ExprVector n = arc->Normal()->NormalExprsN();
597
            ExprVector u = aos.WithMagnitude(Expr::From(1.0));
598
            ExprVector v = n.Cross(u);
599
            // so in our new csys, we start at (1, 0, 0)
600
            Expr *costheta = aof.Dot(u)->Div(r);
601
            Expr *sintheta = aof.Dot(v)->Div(r);
602

603
            double thetas, thetaf, dtheta;
604
            arc->ArcGetAngles(&thetas, &thetaf, &dtheta);
605
            Expr *theta;
606
            if(dtheta < 3*PI/4) {
607
                theta = costheta->ACos();
608
            } else if(dtheta < 5*PI/4) {
609
                // As the angle crosses pi, cos theta is not invertible;
610
                // so use the sine to stop blowing up
611
                theta = Expr::From(PI)->Minus(sintheta->ASin());
612
            } else {
613
                theta = (Expr::From(2*PI))->Minus(costheta->ACos());
614
            }
615

616
            // And write the equation; r*theta = L
617
            AddEq(l, (r->Times(theta))->Minus(ll), 0);
618
            return;
619
        }
620

621
        case Type::POINTS_COINCIDENT: {
622
            EntityBase *a = SK.GetEntity(ptA);
623
            EntityBase *b = SK.GetEntity(ptB);
624
            if(workplane == EntityBase::FREE_IN_3D) {
625
                ExprVector pa = a->PointGetExprs();
626
                ExprVector pb = b->PointGetExprs();
627
                AddEq(l, pa.x->Minus(pb.x), 0);
628
                AddEq(l, pa.y->Minus(pb.y), 1);
629
                AddEq(l, pa.z->Minus(pb.z), 2);
630
            } else {
631
                Expr *au, *av;
632
                Expr *bu, *bv;
633
                a->PointGetExprsInWorkplane(workplane, &au, &av);
634
                b->PointGetExprsInWorkplane(workplane, &bu, &bv);
635
                AddEq(l, au->Minus(bu), 0);
636
                AddEq(l, av->Minus(bv), 1);
637
            }
638
            return;
639
        }
640

641
        case Type::PT_IN_PLANE:
642
            // This one works the same, whether projected or not.
643
            AddEq(l, PointPlaneDistance(
644
                        SK.GetEntity(ptA)->PointGetExprs(), entityA), 0);
645
            return;
646

647
        case Type::PT_ON_FACE: {
648
            // a plane, n dot (p - p0) = 0
649
            ExprVector p = SK.GetEntity(ptA)->PointGetExprs();
650
            EntityBase *f = SK.GetEntity(entityA);
651
            ExprVector p0 = f->FaceGetPointExprs();
652
            ExprVector n = f->FaceGetNormalExprs();
653
            AddEq(l, (p.Minus(p0)).Dot(n), 0);
654
            return;
655
        }
656

657
        case Type::PT_ON_LINE: {
658
            EntityBase *ln = SK.GetEntity(entityA);
659
            EntityBase *a = SK.GetEntity(ln->point[0]);
660
            EntityBase *b = SK.GetEntity(ln->point[1]);
661
            EntityBase *p = SK.GetEntity(ptA);
662

663
            ExprVector ep = p->PointGetExprsInWorkplane(workplane);
664
            ExprVector ea = a->PointGetExprsInWorkplane(workplane);
665
            ExprVector eb = b->PointGetExprsInWorkplane(workplane);
666

667
            ExprVector ptOnLine = ea.Plus(eb.Minus(ea).ScaledBy(Expr::From(valP)));
668
            ExprVector eq = ptOnLine.Minus(ep);
669

670
            AddEq(l, eq);
671
            return;
672
        }
673

674
        case Type::PT_ON_CIRCLE: {
675
            // This actually constrains the point to lie on the cylinder.
676
            EntityBase *circle = SK.GetEntity(entityA);
677
            ExprVector center = SK.GetEntity(circle->point[0])->PointGetExprs();
678
            ExprVector pt     = SK.GetEntity(ptA)->PointGetExprs();
679
            EntityBase *normal = SK.GetEntity(circle->normal);
680
            ExprVector u = normal->NormalExprsU(),
681
                       v = normal->NormalExprsV();
682

683
            Expr *du = (center.Minus(pt)).Dot(u),
684
                 *dv = (center.Minus(pt)).Dot(v);
685

686
            Expr *r = circle->CircleGetRadiusExpr();
687

688
            AddEq(l, du->Square()->Plus(dv->Square())->Sqrt()->Minus(r), 0);
689
            return;
690
        }
691

692
        case Type::AT_MIDPOINT:
693
            if(workplane == EntityBase::FREE_IN_3D) {
694
                EntityBase *ln = SK.GetEntity(entityA);
695
                ExprVector a = SK.GetEntity(ln->point[0])->PointGetExprs();
696
                ExprVector b = SK.GetEntity(ln->point[1])->PointGetExprs();
697
                ExprVector m = (a.Plus(b)).ScaledBy(Expr::From(0.5));
698

699
                if(ptA.v) {
700
                    ExprVector p = SK.GetEntity(ptA)->PointGetExprs();
701
                    AddEq(l, (m.x)->Minus(p.x), 0);
702
                    AddEq(l, (m.y)->Minus(p.y), 1);
703
                    AddEq(l, (m.z)->Minus(p.z), 2);
704
                } else {
705
                    AddEq(l, PointPlaneDistance(m, entityB), 0);
706
                }
707
            } else {
708
                EntityBase *ln = SK.GetEntity(entityA);
709
                EntityBase *a = SK.GetEntity(ln->point[0]);
710
                EntityBase *b = SK.GetEntity(ln->point[1]);
711

712
                Expr *au, *av, *bu, *bv;
713
                a->PointGetExprsInWorkplane(workplane, &au, &av);
714
                b->PointGetExprsInWorkplane(workplane, &bu, &bv);
715
                Expr *mu = Expr::From(0.5)->Times(au->Plus(bu));
716
                Expr *mv = Expr::From(0.5)->Times(av->Plus(bv));
717

718
                if(ptA.v) {
719
                    EntityBase *p = SK.GetEntity(ptA);
720
                    Expr *pu, *pv;
721
                    p->PointGetExprsInWorkplane(workplane, &pu, &pv);
722
                    AddEq(l, pu->Minus(mu), 0);
723
                    AddEq(l, pv->Minus(mv), 1);
724
                } else {
725
                    ExprVector m = PointInThreeSpace(workplane, mu, mv);
726
                    AddEq(l, PointPlaneDistance(m, entityB), 0);
727
                }
728
            }
729
            return;
730
        
731
        case Type::SYMMETRIC:
732
            if(workplane == EntityBase::FREE_IN_3D) {
733
                EntityBase *plane = SK.GetEntity(entityA);
734
                EntityBase *ea = SK.GetEntity(ptA);
735
                EntityBase *eb = SK.GetEntity(ptB);
736
                ExprVector a = ea->PointGetExprs();
737
                ExprVector b = eb->PointGetExprs();
738

739
                // The midpoint of the line connecting the symmetric points
740
                // lies on the plane of the symmetry.
741
                ExprVector m = (a.Plus(b)).ScaledBy(Expr::From(0.5));
742
                AddEq(l, PointPlaneDistance(m, plane->h), 0);
743

744
                // And projected into the plane of symmetry, the points are
745
                // coincident.
746
                Expr *au, *av, *bu, *bv;
747
                ea->PointGetExprsInWorkplane(plane->h, &au, &av);
748
                eb->PointGetExprsInWorkplane(plane->h, &bu, &bv);
749
                AddEq(l, au->Minus(bu), 1);
750
                AddEq(l, av->Minus(bv), 2);
751
            } else {
752
                EntityBase *plane = SK.GetEntity(entityA);
753
                EntityBase *a = SK.GetEntity(ptA);
754
                EntityBase *b = SK.GetEntity(ptB);
755

756
                Expr *au, *av, *bu, *bv;
757
                a->PointGetExprsInWorkplane(workplane, &au, &av);
758
                b->PointGetExprsInWorkplane(workplane, &bu, &bv);
759
                Expr *mu = Expr::From(0.5)->Times(au->Plus(bu));
760
                Expr *mv = Expr::From(0.5)->Times(av->Plus(bv));
761

762
                ExprVector m = PointInThreeSpace(workplane, mu, mv);
763
                AddEq(l, PointPlaneDistance(m, plane->h), 0);
764

765
                // Construct a vector within the workplane that is normal
766
                // to the symmetry pane's normal (i.e., that lies in the
767
                // plane of symmetry). The line connecting the points is
768
                // perpendicular to that constructed vector.
769
                EntityBase *w = SK.GetEntity(workplane);
770
                ExprVector u = w->Normal()->NormalExprsU();
771
                ExprVector v = w->Normal()->NormalExprsV();
772

773
                ExprVector pa = a->PointGetExprs();
774
                ExprVector pb = b->PointGetExprs();
775
                ExprVector n;
776
                Expr *d;
777
                plane->WorkplaneGetPlaneExprs(&n, &d);
778
                AddEq(l, (n.Cross(u.Cross(v))).Dot(pa.Minus(pb)), 1);
779
            }
780
            return;
781

782
        case Type::SYMMETRIC_HORIZ:
783
        case Type::SYMMETRIC_VERT: {
784
            ssassert(workplane != Entity::FREE_IN_3D,
785
                     "Unexpected horizontal/vertical symmetric constraint in 3d");
786

787
            EntityBase *a = SK.GetEntity(ptA);
788
            EntityBase *b = SK.GetEntity(ptB);
789

790
            Expr *au, *av, *bu, *bv;
791
            a->PointGetExprsInWorkplane(workplane, &au, &av);
792
            b->PointGetExprsInWorkplane(workplane, &bu, &bv);
793

794
            if(type == Type::SYMMETRIC_HORIZ) {
795
                AddEq(l, av->Minus(bv), 0);
796
                AddEq(l, au->Plus(bu), 1);
797
            } else {
798
                AddEq(l, au->Minus(bu), 0);
799
                AddEq(l, av->Plus(bv), 1);
800
            }
801
            return;
802
        }
803

804
        case Type::SYMMETRIC_LINE: {
805
            EntityBase *pa = SK.GetEntity(ptA);
806
            EntityBase *pb = SK.GetEntity(ptB);
807

808
            Expr *pau, *pav, *pbu, *pbv;
809
            pa->PointGetExprsInWorkplane(workplane, &pau, &pav);
810
            pb->PointGetExprsInWorkplane(workplane, &pbu, &pbv);
811

812
            EntityBase *ln = SK.GetEntity(entityA);
813
            EntityBase *la = SK.GetEntity(ln->point[0]);
814
            EntityBase *lb = SK.GetEntity(ln->point[1]);
815
            Expr *lau, *lav, *lbu, *lbv;
816
            la->PointGetExprsInWorkplane(workplane, &lau, &lav);
817
            lb->PointGetExprsInWorkplane(workplane, &lbu, &lbv);
818

819
            Expr *dpu = pbu->Minus(pau), *dpv = pbv->Minus(pav);
820
            Expr *dlu = lbu->Minus(lau), *dlv = lbv->Minus(lav);
821

822
            // The line through the points is perpendicular to the line
823
            // of symmetry.
824
            AddEq(l, (dlu->Times(dpu))->Plus(dlv->Times(dpv)), 0);
825

826
            // And the signed distances of the points to the line are
827
            // equal in magnitude and opposite in sign, so sum to zero
828
            Expr *dista = (dlv->Times(lau->Minus(pau)))->Minus(
829
                          (dlu->Times(lav->Minus(pav))));
830
            Expr *distb = (dlv->Times(lau->Minus(pbu)))->Minus(
831
                          (dlu->Times(lav->Minus(pbv))));
832
            AddEq(l, dista->Plus(distb), 1);
833

834
            return;
835
        }
836

837
        case Type::HORIZONTAL:
838
        case Type::VERTICAL: {
839
            ssassert(workplane != Entity::FREE_IN_3D,
840
                     "Unexpected horizontal/vertical constraint in 3d");
841

842
            hEntity ha, hb;
843
            if(entityA.v) {
844
                EntityBase *e = SK.GetEntity(entityA);
845
                ha = e->point[0];
846
                hb = e->point[1];
847
            } else {
848
                ha = ptA;
849
                hb = ptB;
850
            }
851
            EntityBase *a = SK.GetEntity(ha);
852
            EntityBase *b = SK.GetEntity(hb);
853

854
            Expr *au, *av, *bu, *bv;
855
            a->PointGetExprsInWorkplane(workplane, &au, &av);
856
            b->PointGetExprsInWorkplane(workplane, &bu, &bv);
857

858
            AddEq(l, (type == Type::HORIZONTAL) ? av->Minus(bv) : au->Minus(bu), 0);
859
            return;
860
        }
861

862
        case Type::SAME_ORIENTATION: {
863
            EntityBase *a = SK.GetEntity(entityA);
864
            EntityBase *b = SK.GetEntity(entityB);
865

866
            ExprVector au = a->NormalExprsU(),
867
                       an = a->NormalExprsN();
868
            ExprVector bu = b->NormalExprsU(),
869
                       bv = b->NormalExprsV(),
870
                       bn = b->NormalExprsN();
871

872
            ExprVector eq = VectorsParallel3d(an, bn, valP);
873
            AddEq(l, eq.x, 0);
874
            AddEq(l, eq.y, 1);
875
            AddEq(l, eq.z, 2);
876
            Expr *d1 = au.Dot(bv);
877
            Expr *d2 = au.Dot(bu);
878
            // Allow either orientation for the coordinate system, depending
879
            // on how it was drawn.
880
            if(fabs(d1->Eval()) < fabs(d2->Eval())) {
881
                AddEq(l, d1, 3);
882
            } else {
883
                AddEq(l, d2, 3);
884
            }
885
            return;
886
        }
887

888
        case Type::PERPENDICULAR:
889
        case Type::ANGLE: {
890
            EntityBase *a = SK.GetEntity(entityA);
891
            EntityBase *b = SK.GetEntity(entityB);
892
            ExprVector ae = a->VectorGetExprs();
893
            ExprVector be = b->VectorGetExprs();
894
            if(other) ae = ae.ScaledBy(Expr::From(-1));
895
            Expr *c = DirectionCosine(workplane, ae, be);
896

897
            if(type == Type::ANGLE) {
898
                // The direction cosine is equal to the cosine of the
899
                // specified angle
900
                Expr *rads = exA->Times(Expr::From(PI/180)),
901
                     *rc   = rads->Cos();
902
                double arc = fabs(rc->Eval());
903
                // avoid false detection of inconsistent systems by gaining
904
                // up as the difference in dot products gets small at small
905
                // angles; doubles still have plenty of precision, only
906
                // problem is that rank test
907
                Expr *mult = Expr::From(arc > 0.99 ? 0.01/(1.00001 - arc) : 1);
908
                AddEq(l, (c->Minus(rc))->Times(mult), 0);
909
            } else {
910
                // The dot product (and therefore the direction cosine)
911
                // is equal to zero, perpendicular.
912
                AddEq(l, c, 0);
913
            }
914
            return;
915
        }
916

917
        case Type::EQUAL_ANGLE: {
918
            EntityBase *a = SK.GetEntity(entityA);
919
            EntityBase *b = SK.GetEntity(entityB);
920
            EntityBase *c = SK.GetEntity(entityC);
921
            EntityBase *d = SK.GetEntity(entityD);
922
            ExprVector ae = a->VectorGetExprs();
923
            ExprVector be = b->VectorGetExprs();
924
            ExprVector ce = c->VectorGetExprs();
925
            ExprVector de = d->VectorGetExprs();
926

927
            if(other) ae = ae.ScaledBy(Expr::From(-1));
928

929
            Expr *cab = DirectionCosine(workplane, ae, be);
930
            Expr *ccd = DirectionCosine(workplane, ce, de);
931

932
            AddEq(l, cab->Minus(ccd), 0);
933
            return;
934
        }
935

936
        case Type::ARC_LINE_TANGENT: {
937
            EntityBase *arc  = SK.GetEntity(entityA);
938
            EntityBase *line = SK.GetEntity(entityB);
939

940
            ExprVector ac = SK.GetEntity(arc->point[0])->PointGetExprs();
941
            ExprVector ap =
942
                SK.GetEntity(arc->point[other ? 2 : 1])->PointGetExprs();
943

944
            ExprVector ld = line->VectorGetExprs();
945

946
            // The line is perpendicular to the radius
947
            AddEq(l, ld.Dot(ac.Minus(ap)), 0);
948
            return;
949
        }
950

951
        case Type::CUBIC_LINE_TANGENT: {
952
            EntityBase *cubic = SK.GetEntity(entityA);
953
            EntityBase *line  = SK.GetEntity(entityB);
954

955
            ExprVector a;
956
            if(other) {
957
                a = cubic->CubicGetFinishTangentExprs();
958
            } else {
959
                a = cubic->CubicGetStartTangentExprs();
960
            }
961

962
            ExprVector b = line->VectorGetExprs();
963

964
            if(workplane == EntityBase::FREE_IN_3D) {
965
                ExprVector eq = VectorsParallel3d(a, b, valP);
966
                AddEq(l, eq);
967
            } else {
968
                EntityBase *w = SK.GetEntity(workplane);
969
                ExprVector wn = w->Normal()->NormalExprsN();
970
                AddEq(l, (a.Cross(b)).Dot(wn), 0);
971
            }
972
            return;
973
        }
974

975
        case Type::CURVE_CURVE_TANGENT: {
976
            bool parallel = true;
977
            int i;
978
            ExprVector dir[2];
979
            for(i = 0; i < 2; i++) {
980
                EntityBase *e = SK.GetEntity((i == 0) ? entityA : entityB);
981
                bool oth = (i == 0) ? other : other2;
982

983
                if(e->type == Entity::Type::ARC_OF_CIRCLE) {
984
                    ExprVector center, endpoint;
985
                    center = SK.GetEntity(e->point[0])->PointGetExprs();
986
                    endpoint =
987
                        SK.GetEntity(e->point[oth ? 2 : 1])->PointGetExprs();
988
                    dir[i] = endpoint.Minus(center);
989
                    // We're using the vector from the center of the arc to
990
                    // an endpoint; so that's normal to the tangent, not
991
                    // parallel.
992
                    parallel = !parallel;
993
                } else if(e->type == Entity::Type::CUBIC) { // BRANCH_ALWAYS_TAKEN
994
                    if(oth) {
995
                        dir[i] = e->CubicGetFinishTangentExprs();
996
                    } else {
997
                        dir[i] = e->CubicGetStartTangentExprs();
998
                    }
999
                } else {
1000
                    ssassert(false, "Unexpected entity types for CURVE_CURVE_TANGENT");
1001
                }
1002
            }
1003
            if(parallel) {
1004
                EntityBase *w = SK.GetEntity(workplane);
1005
                ExprVector wn = w->Normal()->NormalExprsN();
1006
                AddEq(l, ((dir[0]).Cross(dir[1])).Dot(wn), 0);
1007
            } else {
1008
                AddEq(l, (dir[0]).Dot(dir[1]), 0);
1009
            }
1010
            return;
1011
        }
1012

1013
        case Type::PARALLEL: {
1014
            EntityBase *ea = SK.GetEntity(entityA), *eb = SK.GetEntity(entityB);
1015
            ExprVector a = ea->VectorGetExprsInWorkplane(workplane);
1016
            ExprVector b = eb->VectorGetExprsInWorkplane(workplane);
1017

1018
            if(workplane == EntityBase::FREE_IN_3D) {
1019
                ExprVector eq = VectorsParallel3d(a, b, valP);
1020
                AddEq(l, eq);
1021
            } else {
1022
                // We use expressions written in workplane csys, so we can assume the workplane
1023
                // normal is (0, 0, 1). We can write the equation as:
1024
                //   Expr *eq = a.Cross(b).Dot(ExprVector::From(0.0, 0.0, 1.0));
1025
                // but this will just result in elimination of x and y terms after dot product.
1026
                // We can only use the z expression:
1027
                //   Expr *eq = a.Cross(b).z;
1028
                // but it's more efficient to write it in the terms of pseudo-scalar product:
1029
                Expr *eq = (a.x->Times(b.y))->Minus(a.y->Times(b.x));
1030
                AddEq(l, eq, 0);
1031
            }
1032

1033
            return;
1034
        }
1035

1036
        case Type::WHERE_DRAGGED: {
1037
            EntityBase *ep = SK.GetEntity(ptA);
1038
            if(workplane == EntityBase::FREE_IN_3D) {
1039
                ExprVector ev = ep->PointGetExprs();
1040
                Vector v = ep->PointGetNum();
1041

1042
                AddEq(l, ev.x->Minus(Expr::From(v.x)), 0);
1043
                AddEq(l, ev.y->Minus(Expr::From(v.y)), 1);
1044
                AddEq(l, ev.z->Minus(Expr::From(v.z)), 2);
1045
            } else {
1046
                Expr *u, *v;
1047
                ep->PointGetExprsInWorkplane(workplane, &u, &v);
1048
                AddEq(l, u->Minus(Expr::From(u->Eval())), 0);
1049
                AddEq(l, v->Minus(Expr::From(v->Eval())), 1);
1050
            }
1051
            return;
1052
        }
1053

1054
        case Type::COMMENT:
1055
            return;
1056
    }
1057
    ssassert(false, "Unexpected constraint ID");
1058
}
1059

1060