From 984956cbc7152bc872b97b195ce61f2c1ca4215f Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Thu, 8 Jan 2009 09:22:59 -0800 Subject: [PATCH] Add a constraint to make line length equal arc length. That's quite tricky; can't just use the dot product, since that blows up when you cross pi radians. A gear shift approach, use either sin or cos, same kind of thing as the 3d-parallel constraint. And report a NaN constraint as unconverged, of course. [git-p4: depot-paths = "//depot/solvespace/": change = 1890] --- constraint.cpp | 14 +++++++++++- constrainteq.cpp | 43 +++++++++++++++++++++++++++++++++++ drawconstraint.cpp | 56 +++++++++++++++++++++++++++++----------------- expr.cpp | 15 +++++++++++++ expr.h | 4 ++++ sketch.h | 2 ++ system.cpp | 2 +- textwin.cpp | 3 +++ 8 files changed, 117 insertions(+), 22 deletions(-) diff --git a/constraint.cpp b/constraint.cpp index 3350857..2fb8836 100644 --- a/constraint.cpp +++ b/constraint.cpp @@ -36,6 +36,7 @@ char *Constraint::DescriptionString(void) { case PERPENDICULAR: s = "perpendicular"; break; case EQUAL_RADIUS: s = "eq-radius"; break; case EQUAL_ANGLE: s = "eq-angle"; break; + case EQUAL_LINE_ARC_LEN:s = "eq-line-len-arc-len"; break; case COMMENT: s = "comment"; break; default: s = "???"; break; } @@ -239,6 +240,15 @@ void Constraint::MenuConstrain(int id) { c.type = EQUAL_RADIUS; c.entityA = gs.entity[0]; c.entityB = gs.entity[1]; + } else if(gs.arcs == 1 && gs.lineSegments == 1 && gs.n == 2) { + c.type = EQUAL_LINE_ARC_LEN; + if(SS.GetEntity(gs.entity[0])->type == Entity::ARC_OF_CIRCLE) { + c.entityA = gs.entity[1]; + c.entityB = gs.entity[0]; + } else { + c.entityA = gs.entity[0]; + c.entityB = gs.entity[1]; + } } else { Error("Bad selection for equal length / radius constraint. " "This constraint can apply to:\r\n\r\n" @@ -253,7 +263,9 @@ void Constraint::MenuConstrain(int id) { "(equal angle between A,B and C,D)\r\n" " * three line segments or normals " "(equal angle between A,B and B,C)\r\n" - " * two circles or arcs (equal radius)\r\n"); + " * two circles or arcs (equal radius)\r\n" + " * a line segment and an arc " + "(line segment length equals arc length)\r\n"); return; } if(c.type == EQUAL_ANGLE) { diff --git a/constrainteq.cpp b/constrainteq.cpp index 65f75e6..f03277f 100644 --- a/constrainteq.cpp +++ b/constrainteq.cpp @@ -249,6 +249,49 @@ void Constraint::GenerateReal(IdList *l) { break; } + case EQUAL_LINE_ARC_LEN: { + Entity *line = SS.GetEntity(entityA), + *arc = SS.GetEntity(entityB); + + // Get the line length + ExprVector l0 = SS.GetEntity(line->point[0])->PointGetExprs(), + l1 = SS.GetEntity(line->point[1])->PointGetExprs(); + Expr *ll = (l1.Minus(l0)).Magnitude(); + + // And get the arc radius, and the cosine of its angle + Entity *ao = SS.GetEntity(arc->point[0]), + *as = SS.GetEntity(arc->point[1]), + *af = SS.GetEntity(arc->point[2]); + + ExprVector aos = (as->PointGetExprs()).Minus(ao->PointGetExprs()), + aof = (af->PointGetExprs()).Minus(ao->PointGetExprs()); + Expr *r = aof.Magnitude(); + + ExprVector n = arc->Normal()->NormalExprsN(); + ExprVector u = aos.WithMagnitude(Expr::From(1.0)); + ExprVector v = n.Cross(u); + // so in our new csys, we start at (1, 0, 0) + Expr *costheta = aof.Dot(u)->Div(r); + Expr *sintheta = aof.Dot(v)->Div(r); + + double thetas, thetaf, dtheta; + arc->ArcGetAngles(&thetas, &thetaf, &dtheta); + Expr *theta; + if(dtheta < 3*PI/4) { + theta = costheta->ACos(); + } else if(dtheta < 5*PI/4) { + // As the angle crosses pi, cos theta is not invertible; + // so use the sine to stop blowing up + theta = Expr::From(PI)->Minus(sintheta->ASin()); + } else { + theta = (Expr::From(2*PI))->Minus(costheta->ACos()); + } + + // And write the equation; r*theta = L + AddEq(l, (r->Times(theta))->Minus(ll), 0); + break; + } + case POINTS_COINCIDENT: { Entity *a = SS.GetEntity(ptA); Entity *b = SS.GetEntity(ptB); diff --git a/drawconstraint.cpp b/drawconstraint.cpp index 8eb4e67..b6017d6 100644 --- a/drawconstraint.cpp +++ b/drawconstraint.cpp @@ -104,6 +104,30 @@ void Constraint::DoEqualLenTicks(Vector a, Vector b, Vector gn) { LineDrawOrGetDistance(m.Minus(n), m.Plus(n)); } +void Constraint::DoEqualRadiusTicks(hEntity he) { + Entity *circ = SS.GetEntity(he); + + Vector center = SS.GetEntity(circ->point[0])->PointGetNum(); + double r = circ->CircleGetRadiusNum(); + Quaternion q = circ->Normal()->NormalGetNum(); + Vector u = q.RotationU(), v = q.RotationV(); + + double theta; + if(circ->type == Entity::CIRCLE) { + theta = PI/2; + } else if(circ->type == Entity::ARC_OF_CIRCLE) { + double thetaa, thetab, dtheta; + circ->ArcGetAngles(&thetaa, &thetab, &dtheta); + theta = thetaa + dtheta/2; + } else oops(); + + Vector d = u.ScaledBy(cos(theta)).Plus(v.ScaledBy(sin(theta))); + d = d.ScaledBy(r); + Vector p = center.Plus(d); + Vector tick = d.WithMagnitude(10/SS.GW.scale); + LineDrawOrGetDistance(p.Plus(tick), p.Minus(tick)); +} + void Constraint::DoArcForAngle(Vector a0, Vector da, Vector b0, Vector db, Vector offset, Vector *ref) { @@ -505,30 +529,22 @@ void Constraint::DrawOrGetDistance(Vector *labelPos) { case EQUAL_RADIUS: { for(int i = 0; i < 2; i++) { - Entity *circ = SS.GetEntity(i == 0 ? entityA : entityB); - Vector center = SS.GetEntity(circ->point[0])->PointGetNum(); - double r = circ->CircleGetRadiusNum(); - Quaternion q = circ->Normal()->NormalGetNum(); - Vector u = q.RotationU(), v = q.RotationV(); - - double theta; - if(circ->type == Entity::CIRCLE) { - theta = PI/2; - } else if(circ->type == Entity::ARC_OF_CIRCLE) { - double thetaa, thetab, dtheta; - circ->ArcGetAngles(&thetaa, &thetab, &dtheta); - theta = thetaa + dtheta/2; - } else oops(); - - Vector d = u.ScaledBy(cos(theta)).Plus(v.ScaledBy(sin(theta))); - d = d.ScaledBy(r); - Vector p = center.Plus(d); - Vector tick = d.WithMagnitude(10/SS.GW.scale); - LineDrawOrGetDistance(p.Plus(tick), p.Minus(tick)); + DoEqualRadiusTicks(i == 0 ? entityA : entityB); } break; } + case EQUAL_LINE_ARC_LEN: { + Entity *line = SS.GetEntity(entityA); + DoEqualLenTicks( + SS.GetEntity(line->point[0])->PointGetNum(), + SS.GetEntity(line->point[1])->PointGetNum(), + gn); + + DoEqualRadiusTicks(entityB); + break; + } + case LENGTH_RATIO: case EQUAL_LENGTH_LINES: { Vector a, b; diff --git a/expr.cpp b/expr.cpp index 817aac3..5254426 100644 --- a/expr.cpp +++ b/expr.cpp @@ -243,6 +243,8 @@ int Expr::Children(void) { case SQUARE: case SIN: case COS: + case ASIN: + case ACOS: return 1; default: oops(); @@ -312,6 +314,8 @@ double Expr::Eval(void) { case SQUARE: { double r = a->Eval(); return r*r; } case SIN: return sin(a->Eval()); case COS: return cos(a->Eval()); + case ACOS: return acos(a->Eval()); + case ASIN: return asin(a->Eval()); default: oops(); } @@ -349,6 +353,13 @@ Expr *Expr::PartialWrt(hParam p) { case SIN: return (a->Cos())->Times(a->PartialWrt(p)); case COS: return ((a->Sin())->Times(a->PartialWrt(p)))->Negate(); + case ASIN: + return (From(1)->Div((From(1)->Minus(a->Square()))->Sqrt())) + ->Times(a->PartialWrt(p)); + case ACOS: + return (From(-1)->Div((From(1)->Minus(a->Square()))->Sqrt())) + ->Times(a->PartialWrt(p)); + default: oops(); } } @@ -421,6 +432,8 @@ Expr *Expr::FoldConstants(void) { case NEGATE: case SIN: case COS: + case ASIN: + case ACOS: if(n->a->op == CONSTANT) { double nv = n->Eval(); n->op = CONSTANT; @@ -526,6 +539,8 @@ p: case SQUARE: App("(square "); a->PrintW(); App(")"); break; case SIN: App("(sin "); a->PrintW(); App(")"); break; case COS: App("(cos "); a->PrintW(); App(")"); break; + case ASIN: App("(asin "); a->PrintW(); App(")"); break; + case ACOS: App("(acos "); a->PrintW(); App(")"); break; default: oops(); } diff --git a/expr.h b/expr.h index b4775b2..da0babe 100644 --- a/expr.h +++ b/expr.h @@ -29,6 +29,8 @@ public: static const int SQUARE = 106; static const int SIN = 107; static const int COS = 108; + static const int ASIN = 109; + static const int ACOS = 110; // Special helpers for when we're parsing an expression from text. // Initially, literals (like a constant number) appear in the same @@ -70,6 +72,8 @@ public: inline Expr *Square(void) { return AnyOp(SQUARE, NULL); } inline Expr *Sin (void) { return AnyOp(SIN, NULL); } inline Expr *Cos (void) { return AnyOp(COS, NULL); } + inline Expr *ASin (void) { return AnyOp(ASIN, NULL); } + inline Expr *ACos (void) { return AnyOp(ACOS, NULL); } Expr *PartialWrt(hParam p); double Eval(void); diff --git a/sketch.h b/sketch.h index bd77362..622f2d8 100644 --- a/sketch.h +++ b/sketch.h @@ -458,6 +458,7 @@ public: static const int EQ_LEN_PT_LINE_D = 52; static const int EQ_PT_LN_DISTANCES = 53; static const int EQUAL_ANGLE = 54; + static const int EQUAL_LINE_ARC_LEN = 55; static const int SYMMETRIC = 60; static const int SYMMETRIC_HORIZ = 61; static const int SYMMETRIC_VERT = 62; @@ -527,6 +528,7 @@ public: void DoLabel(Vector ref, Vector *labelPos, Vector gr, Vector gu); void DoProjectedPoint(Vector *p); void DoEqualLenTicks(Vector a, Vector b, Vector gn); + void DoEqualRadiusTicks(hEntity he); double GetDistance(Point2d mp); Vector GetLabelPos(void); diff --git a/system.cpp b/system.cpp index 5757a56..006ecb8 100644 --- a/system.cpp +++ b/system.cpp @@ -516,7 +516,7 @@ didnt_converge: (g->solved.remove).Clear(); SS.constraint.ClearTags(); for(i = 0; i < eq.n; i++) { - if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) { + if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE || isnan(mat.B.num[i])) { // This constraint is unsatisfied. if(!mat.eq[i].isFromConstraint()) continue; diff --git a/textwin.cpp b/textwin.cpp index 9db08fc..0b8cec0 100644 --- a/textwin.cpp +++ b/textwin.cpp @@ -343,6 +343,9 @@ void TextWindow::DescribeSelection(void) { double r = e->CircleGetRadiusNum(); Printf(true, " diameter = %Fi%s", SS.MmToString(r*2)); Printf(false, " radius = %Fi%s", SS.MmToString(r)); + double thetas, thetaf, dtheta; + e->ArcGetAngles(&thetas, &thetaf, &dtheta); + Printf(false, " arc len = %Fi%s", SS.MmToString(dtheta*r)); break; } case Entity::CIRCLE: {