diff --git a/sketch.h b/sketch.h index fc5d66f..a1181ad 100644 --- a/sketch.h +++ b/sketch.h @@ -554,6 +554,9 @@ public: class hEquation { public: DWORD v; + + inline bool isFromConstraint(void); + inline hConstraint constraint(void); }; class Equation { @@ -599,5 +602,10 @@ inline hRequest hParam::request(void) inline hEquation hConstraint::equation(int i) { hEquation r; r.v = (v << 16) | i; return r; } +inline bool hEquation::isFromConstraint(void) + { if(v & 0xc0000000) return false; else return true; } +inline hConstraint hEquation::constraint(void) + { hConstraint r; r.v = (v >> 16); return r; } + #endif diff --git a/solvespace.h b/solvespace.h index 8df57e2..468021c 100644 --- a/solvespace.h +++ b/solvespace.h @@ -212,7 +212,7 @@ public: } B; } mat; - static const double RANK_MAG_TOLERANCE; + static const double RANK_MAG_TOLERANCE, CONVERGE_TOLERANCE; int CalculateRank(void); static bool SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS], double B[], int N); diff --git a/system.cpp b/system.cpp index 46d0da6..1e36e32 100644 --- a/system.cpp +++ b/system.cpp @@ -1,6 +1,7 @@ #include "solvespace.h" const double System::RANK_MAG_TOLERANCE = 1e-4; +const double System::CONVERGE_TOLERANCE = 1e-10; void System::WriteJacobian(int tag) { int a, i, j; @@ -328,7 +329,7 @@ bool System::NewtonSolve(int tag) { if(isnan(mat.B.num[i])) { return false; } - if(fabs(mat.B.num[i]) > 1e-10) { + if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) { converged = false; break; } @@ -360,25 +361,39 @@ void System::WriteEquationsExceptFor(hConstraint hc, hGroup hg) { } void System::FindWhichToRemoveToFixJacobian(Group *g) { - int i; + int a, i; (g->solved.remove).Clear(); - for(i = 0; i < SS.constraint.n; i++) { - Constraint *c = &(SS.constraint.elem[i]); - if(c->group.v != g->h.v) continue; + for(a = 0; a < 2; a++) { + for(i = 0; i < SS.constraint.n; i++) { + Constraint *c = &(SS.constraint.elem[i]); + if(c->group.v != g->h.v) continue; + if((c->type == Constraint::POINTS_COINCIDENT && a == 0) || + (c->type != Constraint::POINTS_COINCIDENT && a == 1)) + { + // Do the constraints in two passes: first everything but + // the point-coincident constraints, then only those + // constraints (so they appear last in the list). + continue; + } - param.ClearTags(); - eq.Clear(); - WriteEquationsExceptFor(c->h, g->h); - eq.ClearTags(); + param.ClearTags(); + eq.Clear(); + WriteEquationsExceptFor(c->h, g->h); + eq.ClearTags(); - WriteJacobian(0); - EvalJacobian(); + // It's a major speedup to solve the easy ones by substitution here, + // and that doesn't break anything. + SolveBySubstitution(); - int rank = CalculateRank(); - if(rank == mat.m) { - // We fixed it by removing this constraint - (g->solved.remove).Add(&(c->h)); + WriteJacobian(0); + EvalJacobian(); + + int rank = CalculateRank(); + if(rank == mat.m) { + // We fixed it by removing this constraint + (g->solved.remove).Add(&(c->h)); + } } } } @@ -476,6 +491,25 @@ void System::Solve(Group *g) { didnt_converge: g->solved.how = Group::DIDNT_CONVERGE; + (g->solved.remove).Clear(); + SS.constraint.ClearTags(); + for(i = 0; i < eq.n; i++) { + if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) { + // This constraint is unsatisfied. + if(!mat.eq[i].isFromConstraint()) continue; + + hConstraint hc = mat.eq[i].constraint(); + Constraint *c = SS.constraint.FindByIdNoOops(hc); + if(!c) continue; + // Don't double-show constraints that generated multiple + // unsatisfied equations + if(!c->tag) { + (g->solved.remove).Add(&(c->h)); + c->tag = 1; + } + } + } + TextWindow::ReportHowGroupSolved(g->h); } diff --git a/textscreens.cpp b/textscreens.cpp index 65c38ba..5678ffd 100644 --- a/textscreens.cpp +++ b/textscreens.cpp @@ -509,25 +509,26 @@ void TextWindow::ShowGroupSolveInfo(void) { Printf(true, "%FtGROUP %E%s", g->DescriptionString()); switch(g->solved.how) { case Group::DIDNT_CONVERGE: - Printf(true, " %FxSOLVE FAILED!%Fd no convergence"); + Printf(true, "%FxSOLVE FAILED!%Fd no convergence"); + Printf(true, "the following constraints are unsatisfied"); break; - case Group::SINGULAR_JACOBIAN: { + case Group::SINGULAR_JACOBIAN: Printf(true, "%FxSOLVE FAILED!%Fd inconsistent system"); Printf(true, "remove any one of these to fix it"); - for(int i = 0; i < g->solved.remove.n; i++) { - hConstraint hc = g->solved.remove.elem[i]; - Constraint *c = SS.constraint.FindByIdNoOops(hc); - if(!c) continue; - - Printf(false, "%Bp %Fl%Ll%D%f%h%s%E", - (i & 1) ? 'd' : 'a', - c->h.v, (&TextWindow::ScreenSelectConstraint), - (&TextWindow::ScreenHoverConstraint), - c->DescriptionString()); - } break; - } + } + + for(int i = 0; i < g->solved.remove.n; i++) { + hConstraint hc = g->solved.remove.elem[i]; + Constraint *c = SS.constraint.FindByIdNoOops(hc); + if(!c) continue; + + Printf(false, "%Bp %Fl%Ll%D%f%h%s%E", + (i & 1) ? 'd' : 'a', + c->h.v, (&TextWindow::ScreenSelectConstraint), + (&TextWindow::ScreenHoverConstraint), + c->DescriptionString()); } Printf(true, "It may be possible to fix the problem ");