diff --git a/src/Mod/Sketcher/App/Sketch.cpp b/src/Mod/Sketcher/App/Sketch.cpp index bfdc0aa3a..59b6410a7 100644 --- a/src/Mod/Sketcher/App/Sketch.cpp +++ b/src/Mod/Sketcher/App/Sketch.cpp @@ -120,8 +120,11 @@ int Sketch::setUpSketch(const std::vector &GeoList, addConstraints(ConstraintList); GCSsys.clearByTag(-1); - GCSsys.initSolution(Parameters); - return diagnose(); + GCSsys.declareUnknowns(Parameters); + GCSsys.initSolution(); + GCSsys.getConflicting(Conflicting); + GCSsys.getRedundant(Redundant); + return GCSsys.dofsNumber(); } const char* nameByType(Sketch::GeoType type) @@ -1594,7 +1597,7 @@ int Sketch::solve(void) InitParameters[i] = **it; GCSsys.addConstraintEqual(*it, &InitParameters[i], -1); } - GCSsys.initSolution(Parameters); + GCSsys.initSolution(); ret = GCSsys.solve(isFine); break; } @@ -1603,8 +1606,11 @@ int Sketch::solve(void) if (ret == GCS::Success) { GCSsys.applySolution(); valid_solution = updateGeometry(); - if (!valid_solution) + if (!valid_solution) { + GCSsys.undoSolution(); + updateGeometry(); Base::Console().Warning("Invalid solution from %s solver.\n", solvername.c_str()); + } } else { valid_solution = false; //Base::Console().Log("NotSolved "); @@ -1630,11 +1636,6 @@ int Sketch::solve(void) } } // soltype - if (!valid_solution) { // undo any changes - GCSsys.undoSolution(); - updateGeometry(); - } - Base::TimeInfo end_time; //Base::Console().Log("T:%s\n",Base::TimeInfo::diffTime(start_time,end_time).c_str()); SolveTime = Base::TimeInfo::diffTimeF(start_time,end_time); @@ -1752,7 +1753,7 @@ int Sketch::initMove(int geoId, PointPos pos, bool fine) } InitParameters = MoveParameters; - GCSsys.initSolution(Parameters); + GCSsys.initSolution(); isInitMove = true; return 0; } @@ -1830,18 +1831,6 @@ Base::Vector3d Sketch::getPoint(int geoId, PointPos pos) return Base::Vector3d(); } -int Sketch::diagnose(void) -{ - Conflicting.clear(); - if (GCSsys.isInit()) { - int dofs = GCSsys.diagnose(Parameters, Conflicting); - return dofs; - } - else { - return -1; - } -} - TopoShape Sketch::toShape(void) const diff --git a/src/Mod/Sketcher/App/Sketch.h b/src/Mod/Sketcher/App/Sketch.h index 3f22f7fdd..3ea95bc4f 100644 --- a/src/Mod/Sketcher/App/Sketch.h +++ b/src/Mod/Sketcher/App/Sketch.h @@ -83,20 +83,10 @@ public: /// retrieves a point Base::Vector3d getPoint(int geoId, PointPos pos); - /** returns the degree of freedom of a sketch and calculates a list of - * conflicting constraints - * - * 0 degrees of freedom correspond to a fully constrained sketch - * -1 degrees of freedom correspond to an over-constrained sketch - * positive degrees of freedom correspond to an under-constrained sketch - * - * an over-constrained sketch will always contain conflicting constraints - * a fully constrained or under-constrained sketch may contain conflicting - * constraints or may not - */ - int diagnose(void); - bool hasConflicts(void) const { return (Conflicting.size() > 0); }; - const std::vector &getConflicting(void) const { return Conflicting; }; + bool hasConflicts(void) const { return (Conflicting.size() > 0); } + const std::vector &getConflicting(void) const { return Conflicting; } + bool hasRedundancies(void) const { return (Redundant.size() > 0); } + const std::vector &getRedundant(void) const { return Redundant; } /** set the datum of a distance or angle constraint to a certain value and solve * This can cause the solving to fail! @@ -213,6 +203,7 @@ protected: GCS::System GCSsys; int ConstraintsCounter; std::vector Conflicting; + std::vector Redundant; // solving parameters std::vector Parameters; // with memory allocation diff --git a/src/Mod/Sketcher/App/SketchObject.cpp b/src/Mod/Sketcher/App/SketchObject.cpp index c9e5b3845..8d75c5f24 100644 --- a/src/Mod/Sketcher/App/SketchObject.cpp +++ b/src/Mod/Sketcher/App/SketchObject.cpp @@ -108,6 +108,11 @@ App::DocumentObjectExecReturn *SketchObject::execute(void) appendConflictMsg(sketch.getConflicting(), msg); return new App::DocumentObjectExecReturn(msg.c_str(),this); } + if (sketch.hasRedundancies()) { // redundant constraints + std::string msg="Sketch with redundant constraints\n"; + appendRedundantMsg(sketch.getRedundant(), msg); + return new App::DocumentObjectExecReturn(msg.c_str(),this); + } // solve the sketch if (sketch.solve() != 0) @@ -1396,6 +1401,20 @@ void SketchObject::appendConflictMsg(const std::vector &conflicting, std::s msg = ss.str(); } +void SketchObject::appendRedundantMsg(const std::vector &redundant, std::string &msg) +{ + std::stringstream ss; + if (msg.length() > 0) + ss << msg; + if (redundant.size() > 0) { + ss << "The following constraints were identified as redundant and should be removed:\n" << redundant[0]; + for (unsigned int i=1; i < redundant.size(); i++) + ss << ", " << redundant[i]; + ss << "\n"; + } + msg = ss.str(); +} + PyObject *SketchObject::getPyObject(void) { if (PythonObject.is(Py::_None())) { diff --git a/src/Mod/Sketcher/App/SketchObject.h b/src/Mod/Sketcher/App/SketchObject.h index 83802d83a..cf865d405 100644 --- a/src/Mod/Sketcher/App/SketchObject.h +++ b/src/Mod/Sketcher/App/SketchObject.h @@ -136,6 +136,8 @@ public: /// generates a warning message about constraint conflicts and appends it to the given message static void appendConflictMsg(const std::vector &conflicting, std::string &msg); + /// generates a warning message about redundant constraints and appends it to the given message + static void appendRedundantMsg(const std::vector &redundant, std::string &msg); // from base class virtual PyObject *getPyObject(void); diff --git a/src/Mod/Sketcher/App/freegcs/GCS.cpp b/src/Mod/Sketcher/App/freegcs/GCS.cpp index c608c9bd0..aa850c9a2 100644 --- a/src/Mod/Sketcher/App/freegcs/GCS.cpp +++ b/src/Mod/Sketcher/App/freegcs/GCS.cpp @@ -41,19 +41,20 @@ typedef boost::adjacency_list Gra // System System::System() -: clist(0), +: plist(0), clist(0), c2p(), p2c(), - subsyslist(0), - reference(), - init(false) + subSystems(0), subSystemsAux(0), + reference(0), + hasUnknowns(false), hasDiagnosis(false), isInit(false) { } System::System(std::vector clist_) -: c2p(), p2c(), - subsyslist(0), - reference(), - init(false) +: plist(0), + c2p(), p2c(), + subSystems(0), subSystemsAux(0), + reference(0), + hasUnknowns(false), hasDiagnosis(false), isInit(false) { // create own (shallow) copy of constraints for (std::vector::iterator constr=clist_.begin(); @@ -125,7 +126,16 @@ System::~System() void System::clear() { - clearReference(); + plist.clear(); + pIndex.clear(); + hasUnknowns = false; + hasDiagnosis = false; + + redundant.clear(); + conflictingTags.clear(); + redundantTags.clear(); + + reference.clear(); clearSubSystems(); free(clist); c2p.clear(); @@ -148,7 +158,9 @@ void System::clearByTag(int tagId) int System::addConstraint(Constraint *constr) { - clearReference(); + isInit = false; + if (constr->getTag() >= 0) // negatively tagged constraints have no impact + hasDiagnosis = false; // on the diagnosis clist.push_back(constr); VEC_pD constr_params = constr->params(); @@ -163,12 +175,15 @@ int System::addConstraint(Constraint *constr) void System::removeConstraint(Constraint *constr) { - clearReference(); - clearSubSystems(); - std::vector::iterator it; it = std::find(clist.begin(), clist.end(), constr); + if (it == clist.end()) + return; + clist.erase(it); + if (constr->getTag() >= 0) + hasDiagnosis = false; + clearSubSystems(); VEC_pD constr_params = c2p[constr]; for (VEC_pD::const_iterator param=constr_params.begin(); @@ -545,9 +560,18 @@ void System::rescaleConstraint(int id, double coeff) clist[id]->rescale(coeff); } - -void System::initSolution(VEC_pD ¶ms) +void System::declareUnknowns(VEC_pD ¶ms) { + plist = params; + pIndex.clear(); + for (int i=0; i < int(plist.size()); ++i) + pIndex[plist[i]] = i; + hasUnknowns = true; +} + +void System::initSolution() +{ + // - Stores the current parameters values in the vector "reference" // - identifies any decoupled subsystems and partitions the original // system into corresponding components // - Stores the current parameters in the vector "reference" @@ -556,18 +580,38 @@ void System::initSolution(VEC_pD ¶ms) // - Organizes the rest of constraints into two subsystems for // tag ids >=0 and < 0 respectively and applies the // system reduction specified in the previous step - MAP_pD_I pIndex; - for (int i=0; i < int(params.size()); ++i) - pIndex[params[i]] = i; + + isInit = false; + if (!hasUnknowns) + return; + + // storing reference configuration + setReference(); + + // diagnose conflicting or redundant constraints + if (!hasDiagnosis) { + diagnose(); + if (!hasDiagnosis) + return; + } + std::vector clistR; + if (redundant.size()) { + for (std::vector::const_iterator constr=clist.begin(); + constr != clist.end(); ++constr) + if (redundant.count(*constr) == 0) + clistR.push_back(*constr); + } + else + clistR = clist; // partitioning into decoupled components Graph g; - for (int i=0; i < int(params.size() + clist.size()); i++) + for (int i=0; i < int(plist.size() + clistR.size()); i++) boost::add_vertex(g); - int cvtid = int(params.size()); - for (std::vector::const_iterator constr=clist.begin(); - constr != clist.end(); ++constr, cvtid++) { + int cvtid = int(plist.size()); + for (std::vector::const_iterator constr=clistR.begin(); + constr != clistR.end(); ++constr, cvtid++) { VEC_pD &cparams = c2p[*constr]; for (VEC_pD::const_iterator param=cparams.begin(); param != cparams.end(); ++param) { @@ -580,113 +624,131 @@ void System::initSolution(VEC_pD ¶ms) VEC_I components(boost::num_vertices(g)); int componentsSize = boost::connected_components(g, &components[0]); - clearReference(); - for (VEC_pD::const_iterator param=params.begin(); - param != params.end(); ++param) - reference[*param] = **param; - // identification of equality constraints and parameter reduction - std::set eliminated; // constraints that will be eliminated through reduction - reductionmap.clear(); + std::set reducedConstrs; // constraints that will be eliminated through reduction + reductionmaps.clear(); // destroy any maps + reductionmaps.resize(componentsSize); // create empty maps to be filled in { - VEC_pD reduced_params=params; + VEC_pD reducedParams=plist; - for (std::vector::const_iterator constr=clist.begin(); - constr != clist.end(); ++constr) { + for (std::vector::const_iterator constr=clistR.begin(); + constr != clistR.end(); ++constr) { if ((*constr)->getTag() >= 0 && (*constr)->getTypeId() == Equal) { MAP_pD_I::const_iterator it1,it2; it1 = pIndex.find((*constr)->params()[0]); it2 = pIndex.find((*constr)->params()[1]); if (it1 != pIndex.end() && it2 != pIndex.end()) { - eliminated.insert(*constr); - double *p_kept = reduced_params[it1->second]; - double *p_replaced = reduced_params[it2->second]; - for (int i=0; i < int(params.size()); ++i) - if (reduced_params[i] == p_replaced) - reduced_params[i] = p_kept; + reducedConstrs.insert(*constr); + double *p_kept = reducedParams[it1->second]; + double *p_replaced = reducedParams[it2->second]; + for (int i=0; i < int(plist.size()); ++i) + if (reducedParams[i] == p_replaced) + reducedParams[i] = p_kept; } } } - for (int i=0; i < int(params.size()); ++i) - if (params[i] != reduced_params[i]) - reductionmap[params[i]] = reduced_params[i]; + for (int i=0; i < int(plist.size()); ++i) + if (plist[i] != reducedParams[i]) { + int cid = components[i]; + reductionmaps[cid][plist[i]] = reducedParams[i]; + } } - std::vector< std::vector > clists0(componentsSize), - clists1(componentsSize), - clists2(componentsSize); - int i = int(params.size()); - for (std::vector::const_iterator constr=clist.begin(); - constr != clist.end(); ++constr, i++) { - if (eliminated.count(*constr) == 0) { - int id = components[i]; - if ((*constr)->getTag() >= 0) - clists0[id].push_back(*constr); - else if ((*constr)->getTag() == -1) // move constraints - clists1[id].push_back(*constr); - else // distance from reference constraints - clists2[id].push_back(*constr); + clists.clear(); // destroy any lists + clists.resize(componentsSize); // create empty lists to be filled in + int i = int(plist.size()); + for (std::vector::const_iterator constr=clistR.begin(); + constr != clistR.end(); ++constr, i++) { + if (reducedConstrs.count(*constr) == 0) { + int cid = components[i]; + clists[cid].push_back(*constr); } } - std::vector< std::vector > plists(componentsSize); - for (int i=0; i < int(params.size()); ++i) { - int id = components[i]; - plists[id].push_back(params[i]); + plists.clear(); // destroy any lists + plists.resize(componentsSize); // create empty lists to be filled in + for (int i=0; i < int(plist.size()); ++i) { + int cid = components[i]; + plists[cid].push_back(plist[i]); } + // calculates subSystems and subSystemsAux from clists, plists and reductionmaps clearSubSystems(); - for (int cid=0; cid < componentsSize; cid++) { - subsyslist.push_back(std::vector(0)); - if (clists0[cid].size() > 0) - subsyslist[cid].push_back(new SubSystem(clists0[cid], plists[cid], reductionmap)); - if (clists1[cid].size() > 0) - subsyslist[cid].push_back(new SubSystem(clists1[cid], plists[cid], reductionmap)); - if (clists2[cid].size() > 0) - subsyslist[cid].push_back(new SubSystem(clists2[cid], plists[cid], reductionmap)); + for (int cid=0; cid < clists.size(); cid++) { + std::vector clist0, clist1; + for (std::vector::const_iterator constr=clists[cid].begin(); + constr != clists[cid].end(); ++constr) { + if ((*constr)->getTag() >= 0) + clist0.push_back(*constr); + else // move or distance from reference constraints + clist1.push_back(*constr); + } + + subSystems.push_back(NULL); + subSystemsAux.push_back(NULL); + if (clist0.size() > 0) + subSystems[cid] = new SubSystem(clist0, plists[cid], reductionmaps[cid]); + if (clist1.size() > 0) + subSystemsAux[cid] = new SubSystem(clist1, plists[cid], reductionmaps[cid]); } - init = true; + + isInit = true; } -void System::clearReference() +void System::setReference() { - init = false; reference.clear(); + reference.reserve(plist.size()); + for (VEC_pD::const_iterator param=plist.begin(); + param != plist.end(); ++param) + reference.push_back(**param); } void System::resetToReference() { - for (MAP_pD_D::const_iterator it=reference.begin(); - it != reference.end(); ++it) - *(it->first) = it->second; + if (reference.size() == plist.size()) { + VEC_D::const_iterator ref=reference.begin(); + VEC_pD::iterator param=plist.begin(); + for (; ref != reference.end(); ++ref, ++param) + **param = *ref; + } } int System::solve(VEC_pD ¶ms, bool isFine, Algorithm alg) { - initSolution(params); + declareUnknowns(params); + initSolution(); return solve(isFine, alg); } int System::solve(bool isFine, Algorithm alg) { + if (!isInit) + return Failed; + bool isReset = false; // return success by default in order to permit coincidence constraints to be applied // even if no other system has to be solved int res = Success; - for (int cid=0; cid < int(subsyslist.size()); cid++) { - if (subsyslist[cid].size() > 0 && !isReset) { + for (int cid=0; cid < int(subSystems.size()); cid++) { + if ((subSystems[cid] || subSystemsAux[cid]) && !isReset) { resetToReference(); isReset = true; } - if (subsyslist[cid].size() == 1) - res = std::max(res, solve(subsyslist[cid][0], isFine, alg)); - else if (subsyslist[cid].size() == 2) - res = std::max(res, solve(subsyslist[cid][0], subsyslist[cid][1], isFine)); - else if (subsyslist[cid].size() > 2) - // subsystem 1 has higher priority than subsystems 2,3,... - // these subsystems act like a preconditioner - for (int i=subsyslist[cid].size()-1; i > 0; i--) - res = std::max(res, solve(subsyslist[cid][0], subsyslist[cid][i], isFine)); + if (subSystems[cid] && subSystemsAux[cid]) + res = std::max(res, solve(subSystems[cid], subSystemsAux[cid], isFine)); + else if (subSystems[cid]) + res = std::max(res, solve(subSystems[cid], isFine, alg)); + else if (subSystemsAux[cid]) + res = std::max(res, solve(subSystemsAux[cid], isFine, alg)); + } + if (res == Success) { + for (std::set::const_iterator constr=redundant.begin(); + constr != redundant.end(); constr++) + if ((*constr)->error() > XconvergenceFine) { + res = Converged; + return res; + } } return res; } @@ -1068,7 +1130,7 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) int xsizeB = subsysB->pSize(); int csizeA = subsysA->cSize(); - VEC_pD plist(xsizeA+xsizeB); + VEC_pD plistAB(xsizeA+xsizeB); { VEC_pD plistA, plistB; subsysA->getParamList(plistA); @@ -1079,10 +1141,10 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) VEC_pD::const_iterator it; it = std::set_union(plistA.begin(),plistA.end(), - plistB.begin(),plistB.end(),plist.begin()); - plist.resize(it-plist.begin()); + plistB.begin(),plistB.end(),plistAB.begin()); + plistAB.resize(it-plistAB.begin()); } - int xsize = plist.size(); + int xsize = plistAB.size(); Eigen::MatrixXd B = Eigen::MatrixXd::Identity(xsize, xsize); Eigen::MatrixXd JA(csizeA, xsize); @@ -1100,12 +1162,12 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) subsysA->redirectParams(); subsysB->redirectParams(); - subsysB->getParams(plist,x); - subsysA->getParams(plist,x); - subsysB->setParams(plist,x); // just to ensure that A and B are synchronized + subsysB->getParams(plistAB,x); + subsysA->getParams(plistAB,x); + subsysB->setParams(plistAB,x); // just to ensure that A and B are synchronized - subsysB->calcGrad(plist,grad); - subsysA->calcJacobi(plist,JA); + subsysB->calcGrad(plistAB,grad); + subsysA->calcJacobi(plistAB,JA); subsysA->calcResidual(resA); double convergence = isFine ? XconvergenceFine : XconvergenceRough; @@ -1130,7 +1192,7 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) double tau=0.5; double rho=0.5; double alpha=1; - alpha = std::min(alpha, subsysA->maxStep(plist,xdir)); + alpha = std::min(alpha, subsysA->maxStep(plistAB,xdir)); // Eq. 18.32 // double mu = lambda.lpNorm() + 0.01; @@ -1148,8 +1210,8 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) double deriv = grad.dot(xdir) - mu * resA.lpNorm<1>(); x = x0 + alpha * xdir; - subsysA->setParams(plist,x); - subsysB->setParams(plist,x); + subsysA->setParams(plistAB,x); + subsysB->setParams(plistAB,x); subsysA->calcResidual(resA); double f = subsysB->error() + mu * resA.lpNorm<1>(); @@ -1161,8 +1223,8 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) // Eigen::ComputeThinV).solve(-resA); xdir1 = -Y*resA; x += xdir1; // = x0 + alpha * xdir + xdir1 - subsysA->setParams(plist,x); - subsysB->setParams(plist,x); + subsysA->setParams(plistAB,x); + subsysB->setParams(plistAB,x); subsysA->calcResidual(resA); f = subsysB->error() + mu * resA.lpNorm<1>(); if (f < f0 + eta * alpha * deriv) @@ -1172,8 +1234,8 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) if (alpha < 1e-8) // let the linesearch fail alpha = 0.; x = x0 + alpha * xdir; - subsysA->setParams(plist,x); - subsysB->setParams(plist,x); + subsysA->setParams(plistAB,x); + subsysB->setParams(plistAB,x); subsysA->calcResidual(resA); f = subsysB->error() + mu * resA.lpNorm<1>(); if (alpha < 1e-8) // let the linesearch fail @@ -1186,8 +1248,8 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) y = grad - JA.transpose() * lambda; { - subsysB->calcGrad(plist,grad); - subsysA->calcJacobi(plist,JA); + subsysB->calcGrad(plistAB,grad); + subsysA->calcJacobi(plistAB,JA); subsysA->calcResidual(resA); } y = grad - JA.transpose() * lambda - y; // Eq. 18.13 @@ -1225,13 +1287,15 @@ int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine) void System::applySolution() { - for (int cid=0; cid < int(subsyslist.size()); cid++) - for (int i=subsyslist[cid].size()-1; i >= 0; i--) - subsyslist[cid][i]->applySolution(); - - for (MAP_pD_pD::const_iterator it=reductionmap.begin(); - it != reductionmap.end(); ++it) - *(it->first) = *(it->second); + for (int cid=0; cid < int(subSystems.size()); cid++) { + if (subSystemsAux[cid]) + subSystemsAux[cid]->applySolution(); + if (subSystems[cid]) + subSystems[cid]->applySolution(); + for (MAP_pD_pD::const_iterator it=reductionmaps[cid].begin(); + it != reductionmaps[cid].end(); ++it) + *(it->first) = *(it->second); + } } void System::undoSolution() @@ -1239,7 +1303,7 @@ void System::undoSolution() resetToReference(); } -int System::diagnose(VEC_pD ¶ms, VEC_I &conflictingTags) +int System::diagnose() { // Analyses the constrainess grad of the system and provides feedback // The vector "conflictingTags" will hold a group of conflicting constraints @@ -1251,22 +1315,24 @@ int System::diagnose(VEC_pD ¶ms, VEC_I &conflictingTags) // will provide no feedback about possible conflicts between // two high priority constraints. For this reason, tagging // constraints with 0 should be used carefully. - if (!isInit()) - return -1; + hasDiagnosis = false; + if (!hasUnknowns) { + dofs = -1; + return dofs; + } + redundant.clear(); conflictingTags.clear(); - std::vector conflictingIndex; - VEC_I tags; - Eigen::MatrixXd J(clist.size(), params.size()); + redundantTags.clear(); + Eigen::MatrixXd J(clist.size(), plist.size()); int count=0; for (std::vector::iterator constr=clist.begin(); constr != clist.end(); ++constr) { (*constr)->revertParams(); if ((*constr)->getTag() >= 0) { count++; - tags.push_back((*constr)->getTag()); - for (int j=0; j < int(params.size()); j++) - J(count-1,j) = (*constr)->grad(params[j]); + for (int j=0; j < int(plist.size()); j++) + J(count-1,j) = (*constr)->grad(plist[j]); } } @@ -1284,7 +1350,7 @@ int System::diagnose(VEC_pD ¶ms, VEC_I &conflictingTags) R = qrJT.matrixQR().topRows(constrNum) .triangularView(); - if (constrNum > rank) { // conflicting constraints + if (constrNum > rank) { // conflicting or redundant constraints for (int i=1; i < rank; i++) { // eliminate non zeros above pivot assert(R(i,i) != 0); @@ -1296,43 +1362,139 @@ int System::diagnose(VEC_pD ¶ms, VEC_I &conflictingTags) } } } - conflictingIndex.resize(constrNum-rank); + std::vector< std::vector > conflictGroups(constrNum-rank); for (int j=rank; j < constrNum; j++) { for (int row=0; row < rank; row++) { if (fabs(R(row,j)) > 1e-10) { int origCol = qrJT.colsPermutation().indices()[row]; - conflictingIndex[j-rank].push_back(origCol); + conflictGroups[j-rank].push_back(clist[origCol]); } } int origCol = qrJT.colsPermutation().indices()[j]; - conflictingIndex[j-rank].push_back(origCol); + conflictGroups[j-rank].push_back(clist[origCol]); } - SET_I tags_set; - for (int i=0; i < conflictingIndex.size(); i++) { - for (int j=0; j < conflictingIndex[i].size(); j++) { - tags_set.insert(tags[conflictingIndex[i][j]]); + // try to remove the conflicting constraints and solve the + // system in order to check if the removed constraints were + // just redundant but not really conflicting + std::set skipped; + SET_I satisfiedGroups; + while (1) { + std::map< Constraint *, SET_I > conflictingMap; + for (int i=0; i < conflictGroups.size(); i++) { + if (satisfiedGroups.count(i) == 0) { + for (int j=0; j < conflictGroups[i].size(); j++) { + Constraint *constr = conflictGroups[i][j]; + if (constr->getTag() != 0) // exclude constraints tagged with zero + conflictingMap[constr].insert(i); + } + } + } + if (conflictingMap.empty()) + break; + + int maxPopularity = 0; + Constraint *mostPopular = NULL; + for (std::map< Constraint *, SET_I >::const_iterator it=conflictingMap.begin(); + it != conflictingMap.end(); it++) { + if (it->second.size() >= maxPopularity) { + mostPopular = it->first; + maxPopularity = it->second.size(); + } + } + if (maxPopularity > 0) { + skipped.insert(mostPopular); + for (SET_I::const_iterator it=conflictingMap[mostPopular].begin(); + it != conflictingMap[mostPopular].end(); it++) + satisfiedGroups.insert(*it); } } - tags_set.erase(0); // exclude constraints tagged with zero - conflictingTags.resize(tags_set.size()); - std::copy(tags_set.begin(), tags_set.end(), conflictingTags.begin()); - if (paramsNum == rank) // over-constrained - return paramsNum - constrNum; + std::vector clistTmp; + clistTmp.reserve(clist.size()); + for (std::vector::iterator constr=clist.begin(); + constr != clist.end(); ++constr) + if (skipped.count(*constr) == 0) + clistTmp.push_back(*constr); + + SubSystem *subSysTmp = new SubSystem(clistTmp, plist); + int res = solve(subSysTmp); + if (res == Success) { + subSysTmp->applySolution(); + for (std::set::const_iterator constr=skipped.begin(); + constr != skipped.end(); constr++) { + if ((*constr)->error() < XconvergenceFine) + redundant.insert(*constr); + } + resetToReference(); + + std::vector< std::vector > conflictGroupsOrig=conflictGroups; + conflictGroups.clear(); + for (int i=conflictGroupsOrig.size()-1; i >= 0; i--) { + bool isRedundant = false; + for (int j=0; j < conflictGroupsOrig[i].size(); j++) { + if (redundant.count(conflictGroupsOrig[i][j]) > 0) { + isRedundant = true; + break; + } + } + if (!isRedundant) + conflictGroups.push_back(conflictGroupsOrig[i]); + else + constrNum--; + } + } + delete subSysTmp; + + // simplified output of conflicting tags + SET_I conflictingTagsSet; + for (int i=0; i < conflictGroups.size(); i++) { + for (int j=0; j < conflictGroups[i].size(); j++) { + conflictingTagsSet.insert(conflictGroups[i][j]->getTag()); + } + } + conflictingTagsSet.erase(0); // exclude constraints tagged with zero + conflictingTags.resize(conflictingTagsSet.size()); + std::copy(conflictingTagsSet.begin(), conflictingTagsSet.end(), + conflictingTags.begin()); + + // output of redundant tags + SET_I redundantTagsSet; + for (std::set::iterator constr=redundant.begin(); + constr != redundant.end(); ++constr) + redundantTagsSet.insert((*constr)->getTag()); + // remove tags represented at least in one non-redundant constraint + for (std::vector::iterator constr=clist.begin(); + constr != clist.end(); ++constr) + if (redundant.count(*constr) == 0) + redundantTagsSet.erase((*constr)->getTag()); + redundantTags.resize(redundantTagsSet.size()); + std::copy(redundantTagsSet.begin(), redundantTagsSet.end(), + redundantTags.begin()); + + if (paramsNum == rank && constrNum > rank) { // over-constrained + hasDiagnosis = true; + dofs = paramsNum - constrNum; + return dofs; + } } - return paramsNum - rank; + hasDiagnosis = true; + dofs = paramsNum - rank; + return dofs; } - return params.size(); + hasDiagnosis = true; + dofs = plist.size(); + return dofs; } void System::clearSubSystems() { - init = false; - for (int i=0; i < int(subsyslist.size()); i++) - free(subsyslist[i]); - subsyslist.clear(); + isInit = false; + free(subSystems); + free(subSystemsAux); + subSystems.clear(); + subSystemsAux.clear(); } double lineSearch(SubSystem *subsys, Eigen::VectorXd &xdir) diff --git a/src/Mod/Sketcher/App/freegcs/GCS.h b/src/Mod/Sketcher/App/freegcs/GCS.h index 5cd1466ea..f0d87906f 100644 --- a/src/Mod/Sketcher/App/freegcs/GCS.h +++ b/src/Mod/Sketcher/App/freegcs/GCS.h @@ -49,25 +49,31 @@ namespace GCS // This is the main class. It holds all constraints and information // about partitioning into subsystems and solution strategies private: - std::vector clist; + VEC_pD plist; // list of the unknown parameters + MAP_pD_I pIndex; + std::vector clist; std::map c2p; // constraint to parameter adjacency list std::map > p2c; // parameter to constraint adjacency list - // each row of subsyslist contains up to 3 subsystems. - // the first one has the highest priority, always used as the primary subsystem - // the second one is used as secondary subsystem - // the third one is used as secondary system and serves as a preconditioner - std::vector< std::vector > subsyslist; + std::vector subSystems, subSystemsAux; void clearSubSystems(); - MAP_pD_D reference; - void clearReference(); - void resetToReference(); + VEC_D reference; + void setReference(); // copies the current parameter values to reference + void resetToReference(); // reverts all parameter values to the stored reference - MAP_pD_pD reductionmap; // for simplification of equality constraints + std::vector< VEC_pD > plists; // partitioned plist except equality constraints + std::vector< std::vector > clists; // partitioned clist except equality constraints + std::vector< MAP_pD_pD > reductionmaps; // for simplification of equality constraints - bool init; + int dofs; + std::set redundant; + VEC_I conflictingTags, redundantTags; + + bool hasUnknowns; // if plist is filled with the unknown parameters + bool hasDiagnosis; // if dofs, conflictingTags, redundantTags are up to date + bool isInit; // if plists, clists, reductionmaps are up to date int solve_BFGS(SubSystem *subsys, bool isFine); int solve_LM(SubSystem *subsys); @@ -147,7 +153,8 @@ namespace GCS int addConstraintP2PSymmetric(Point &p1, Point &p2, Line &l, int tagId=0); void rescaleConstraint(int id, double coeff); - void initSolution(VEC_pD ¶ms); + void declareUnknowns(VEC_pD ¶ms); + void initSolution(); int solve(bool isFine=true, Algorithm alg=DogLeg); int solve(VEC_pD ¶ms, bool isFine=true, Algorithm alg=DogLeg); @@ -157,9 +164,12 @@ namespace GCS void applySolution(); void undoSolution(); - bool isInit() const { return init; } - - int diagnose(VEC_pD ¶ms, VEC_I &conflictingTags); + int diagnose(); + int dofsNumber() { return hasDiagnosis ? dofs : -1; } + void getConflicting(VEC_I &conflictingOut) const + { conflictingOut = hasDiagnosis ? conflictingTags : VEC_I(0); } + void getRedundant(VEC_I &redundantOut) const + { redundantOut = hasDiagnosis ? redundantTags : VEC_I(0); } }; /////////////////////////////////////// diff --git a/src/Mod/Sketcher/Gui/TaskSketcherMessages.cpp b/src/Mod/Sketcher/Gui/TaskSketcherMessages.cpp index 7061f3a97..5ef6c65fb 100644 --- a/src/Mod/Sketcher/Gui/TaskSketcherMessages.cpp +++ b/src/Mod/Sketcher/Gui/TaskSketcherMessages.cpp @@ -86,6 +86,9 @@ void TaskSketcherMessages::slotSetUp(int type, int dofs, const std::string &msg) case 3: ui->labelConstrainStatus->setText(QString::fromLatin1("Over-constrained sketch
%1").arg(QString::fromStdString(msg))); break; + case 4: + ui->labelConstrainStatus->setText(QString::fromLatin1("Sketch contains redundant constraints
%1").arg(QString::fromStdString(msg))); + break; } } diff --git a/src/Mod/Sketcher/Gui/ViewProviderSketch.cpp b/src/Mod/Sketcher/Gui/ViewProviderSketch.cpp index 78841d18c..2239f7445 100644 --- a/src/Mod/Sketcher/Gui/ViewProviderSketch.cpp +++ b/src/Mod/Sketcher/Gui/ViewProviderSketch.cpp @@ -2825,21 +2825,30 @@ void ViewProviderSketch::solveSketch(void) signalSetUp(2, dofs, msg); signalSolved(-1, 0); } - else if (edit->ActSketch.solve() == 0) { // solving the sketch - if (dofs == 0) { - // color the sketch as fully constrained - edit->FullyConstrained = true; - //Base::Console().Message("Fully constrained sketch\n"); - signalSetUp(0, 0, msg); + else { + if (edit->ActSketch.hasRedundancies()) { // redundant constraints + SketchObject::appendRedundantMsg(edit->ActSketch.getRedundant(), msg); + //Base::Console().Warning("Sketch with redundant constraints\n%s",msg.c_str()); + signalSetUp(4, dofs, msg); + } + if (edit->ActSketch.solve() == 0) { // solving the sketch + if (dofs == 0) { + // color the sketch as fully constrained + edit->FullyConstrained = true; + if (!edit->ActSketch.hasRedundancies()) { + //Base::Console().Message("Fully constrained sketch\n"); + signalSetUp(0, 0, msg); + } + } + else if (!edit->ActSketch.hasRedundancies()) { + //Base::Console().Message("Under-constrained sketch with %d degrees of freedom\n", dofs); + signalSetUp(1, dofs, msg); + } + signalSolved(0, edit->ActSketch.SolveTime); } else { - //Base::Console().Message("Under-constrained sketch with %d degrees of freedom\n", dofs); - signalSetUp(1, dofs, msg); + signalSolved(1, edit->ActSketch.SolveTime); } - signalSolved(0, edit->ActSketch.SolveTime); - } - else { - signalSolved(1, edit->ActSketch.SolveTime); } }