Salome HOME
updated copyright message
[modules/shaper.git] / src / SketchSolver / PlaneGCSSolver / PlaneGCSSolver_Solver.cpp
index a49487878b7dbc309abfa6a8e21ef17e2513ceb9..47b03a01cd61570aa52253ced5ed8906ba091475 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2014-2017  CEA/DEN, EDF R&D
+// Copyright (C) 2014-2023  CEA, EDF
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 //
 // You should have received a copy of the GNU Lesser General Public
 // License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-// See http://www.salome-platform.org/ or
-// email : webmaster.salome@opencascade.com<mailto:webmaster.salome@opencascade.com>
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 
 #include <PlaneGCSSolver_Solver.h>
 #include <Events_LongOp.h>
 
 // Multiplier to correlate IDs of SketchPlugin constraint and primitive PlaneGCS constraints
-static const int THE_CONSTRAINT_MULT = 10;
+static const int THE_CONSTRAINT_MULT = 100;
 
 
 PlaneGCSSolver_Solver::PlaneGCSSolver_Solver()
@@ -63,11 +62,11 @@ void PlaneGCSSolver_Solver::addConstraint(const ConstraintID& theMultiConstraint
     GCSConstraintPtr aConstraint = *anIt;
     aConstraint->setTag(anID);
     myEquationSystem->addConstraint(aConstraint.get());
-    myConstraints[theMultiConstraintID].insert(aConstraint);
 
     if (anID > CID_UNKNOWN)
       ++anID;
   }
+  myConstraints[theMultiConstraintID] = theConstraints;
 
   if (theMultiConstraintID >= CID_UNKNOWN)
     myDOF = -1;
@@ -78,7 +77,7 @@ void PlaneGCSSolver_Solver::removeConstraint(const ConstraintID& theID)
 {
   ConstraintMap::iterator aFound = myConstraints.find(theID);
   if (aFound != myConstraints.end()) {
-    for (std::set<GCSConstraintPtr>::iterator anIt = aFound->second.begin();
+    for (std::list<GCSConstraintPtr>::iterator anIt = aFound->second.begin();
          anIt != aFound->second.end(); ++anIt)
       myEquationSystem->clearByTag((*anIt)->getTag());
 
@@ -165,6 +164,17 @@ PlaneGCSSolver_Solver::SolveStatus PlaneGCSSolver_Solver::solve()
     diagnose();
     aResult = (GCS::SolveStatus)myEquationSystem->solve(myParameters);
   }
+
+  if (aResult == GCS::Failed) {
+    // DogLeg solver failed without conflicting constraints, try to use Levenberg-Marquardt solver
+    diagnose(GCS::LevenbergMarquardt);
+    aResult = (GCS::SolveStatus)myEquationSystem->solve(myParameters, true,
+                                                        GCS::LevenbergMarquardt);
+    if (aResult == GCS::Failed) {
+      diagnose(GCS::BFGS);
+      aResult = (GCS::SolveStatus)myEquationSystem->solve(myParameters, true, GCS::BFGS);
+    }
+  }
   Events_LongOp::end(this);
 
   // collect information about conflicting constraints every time,
@@ -173,30 +183,6 @@ PlaneGCSSolver_Solver::SolveStatus PlaneGCSSolver_Solver::solve()
   collectConflicting();
   if (!myConflictingIDs.empty())
     aResult = GCS::Failed;
-  else if (aResult == GCS::Failed) {
-    // DogLeg solver failed without conflicting constraints, try to use Levenberg-Marquardt solver
-    // if there are point-line distance constraints
-    ConstraintMap::iterator aCIt = myConstraints.begin();
-    for (; aCIt != myConstraints.end(); ++aCIt) {
-      if (aCIt->second.size() <= 1)
-        continue;
-      std::set<GCSConstraintPtr>::const_iterator anIt = aCIt->second.begin();
-      for (; anIt != aCIt->second.end(); ++anIt)
-        if ((*anIt)->getTypeId() == GCS::P2LDistance)
-          break;
-      if (anIt != aCIt->second.end())
-        break;
-    }
-
-    if (aCIt != myConstraints.end()) {
-      aResult = (GCS::SolveStatus)myEquationSystem->solve(
-          myParameters, true, GCS::LevenbergMarquardt);
-      myConfCollected = false;
-      collectConflicting();
-      if (!myConflictingIDs.empty())
-        aResult = GCS::Failed;
-    }
-  }
 
   SolveStatus aStatus;
   if (aResult == GCS::Failed)
@@ -225,7 +211,7 @@ bool PlaneGCSSolver_Solver::isConflicting(const ConstraintID& theConstraint) con
   return myConflictingIDs.find((int)theConstraint) != myConflictingIDs.end();
 }
 
-void PlaneGCSSolver_Solver::collectConflicting()
+void PlaneGCSSolver_Solver::collectConflicting(bool withRedundant)
 {
   GCS::VEC_I aConflict;
   myEquationSystem->getConflicting(aConflict);
@@ -233,10 +219,12 @@ void PlaneGCSSolver_Solver::collectConflicting()
   for (GCS::VEC_I::const_iterator anIt = aConflict.begin(); anIt != aConflict.end(); ++anIt)
     myConflictingIDs.insert((*anIt) / THE_CONSTRAINT_MULT);
 
-  myEquationSystem->getRedundant(aConflict);
-  // convert PlaneGCS constraint IDs to SketchPlugin's ID
-  for (GCS::VEC_I::const_iterator anIt = aConflict.begin(); anIt != aConflict.end(); ++anIt)
-    myConflictingIDs.insert((*anIt) / THE_CONSTRAINT_MULT);
+  if (withRedundant) {
+    myEquationSystem->getRedundant(aConflict);
+    // convert PlaneGCS constraint IDs to SketchPlugin's ID
+    for (GCS::VEC_I::const_iterator anIt = aConflict.begin(); anIt != aConflict.end(); ++anIt)
+      myConflictingIDs.insert((*anIt) / THE_CONSTRAINT_MULT);
+  }
 
   myConfCollected = true;
 }
@@ -248,19 +236,115 @@ int PlaneGCSSolver_Solver::dof()
   return myDOF;
 }
 
-void PlaneGCSSolver_Solver::diagnose()
+void PlaneGCSSolver_Solver::diagnose(const GCS::Algorithm& theAlgo)
 {
   myEquationSystem->declareUnknowns(myParameters);
-  myDOF = myEquationSystem->diagnose();
+  myDOF = myEquationSystem->diagnose(theAlgo);
   myDiagnoseBeforeSolve = false;
 }
 
-void PlaneGCSSolver_Solver::addFictiveConstraintIfNecessary()
+void PlaneGCSSolver_Solver::getFreeParameters(GCS::SET_pD& theFreeParams)
 {
-  if (!myConstraints.empty() &&
-      myConstraints.find(CID_MOVEMENT) == myConstraints.end())
+  if (myConstraints.empty())
+    theFreeParams.insert(myParameters.begin(), myParameters.end());
+  else {
+    GCS::VEC_pD aParametersCopy = myParameters;
+    ConstraintMap aConstraintCopy = myConstraints;
+
+    // clear the set of equations
+    clear();
+    // reset constraints
+    myParameters = aParametersCopy;
+    for (ConstraintMap::iterator anIt = aConstraintCopy.begin();
+         anIt != aConstraintCopy.end(); ++anIt)
+      addConstraint(anIt->first, anIt->second);
+
+    // parameters detection works for Dense QR only
+    GCS::QRAlgorithm aQRAlgo = myEquationSystem->qrAlgorithm;
+    myEquationSystem->qrAlgorithm = GCS::EigenDenseQR;
+    diagnose();
+    GCS::VEC_pD aFreeParams;
+    myEquationSystem->getDependentParams(aFreeParams);
+    theFreeParams.insert(aFreeParams.begin(), aFreeParams.end());
+    // revert QR decomposition algorithm
+    myEquationSystem->qrAlgorithm = aQRAlgo;
+  }
+
+  if (theFreeParams.empty())
     return;
 
+  // find all equal parameters too
+  struct EqualParameters
+  {
+    typedef std::map<double*, std::list<GCS::SET_pD>::iterator> MapParamGroup;
+
+    std::list<GCS::SET_pD> myEqualParams;
+    MapParamGroup myGroups;
+
+    void add(double* theParam1, double* theParam2)
+    {
+      MapParamGroup::iterator aFound1 = myGroups.find(theParam1);
+      MapParamGroup::iterator aFound2 = myGroups.find(theParam2);
+
+      if (aFound1 == myGroups.end()) {
+        if (aFound2 == myGroups.end()) {
+          // create new group
+          myEqualParams.push_back(GCS::SET_pD());
+          std::list<GCS::SET_pD>::iterator aGroup = --myEqualParams.end();
+          aGroup->insert(theParam1);
+          aGroup->insert(theParam2);
+          myGroups[theParam1] = aGroup;
+          myGroups[theParam2] = aGroup;
+        }
+        else {
+          // add first parameter to the second group
+          aFound2->second->insert(theParam1);
+          myGroups[theParam1] = aFound2->second;
+        }
+      }
+      else {
+        if (aFound2 == myGroups.end()) {
+          // add second parameter to the first group
+          aFound1->second->insert(theParam2);
+          myGroups[theParam2] = aFound1->second;
+        }
+        else if (aFound1 != aFound2) {
+          // merge two groups
+          GCS::SET_pD aCopy = *(aFound2->second);
+          myEqualParams.erase(aFound2->second);
+          for (GCS::SET_pD::iterator anIt = aCopy.begin(); anIt != aCopy.end(); ++anIt)
+            myGroups[*anIt] = aFound1->second;
+          aFound1->second->insert(aCopy.begin(), aCopy.end());
+        }
+      }
+    }
+  } anEqualParams;
+
+  for (ConstraintMap::iterator anIt = myConstraints.begin(); anIt != myConstraints.end(); ++anIt)
+    for (std::list<GCSConstraintPtr>::iterator aCIt = anIt->second.begin();
+         aCIt != anIt->second.end(); ++aCIt) {
+      if ((*aCIt)->getTypeId() == GCS::Equal)
+        anEqualParams.add((*aCIt)->params()[0], (*aCIt)->params()[1]);
+    }
+
+  GCS::SET_pD aFreeParamsCopy = theFreeParams;
+  for (GCS::SET_pD::iterator anIt = aFreeParamsCopy.begin();
+       anIt != aFreeParamsCopy.end(); ++anIt) {
+    EqualParameters::MapParamGroup::iterator aFound = anEqualParams.myGroups.find(*anIt);
+    if (aFound != anEqualParams.myGroups.end())
+      theFreeParams.insert(aFound->second->begin(), aFound->second->end());
+  }
+}
+
+void PlaneGCSSolver_Solver::addFictiveConstraintIfNecessary()
+{
+  bool hasOnlyMovement = true;
+  for (ConstraintMap::iterator anIt = myConstraints.begin();
+       anIt != myConstraints.end() && hasOnlyMovement; ++anIt)
+    hasOnlyMovement = anIt->first == CID_MOVEMENT;
+  if (!hasOnlyMovement)
+    return; // regular constraints are available too
+
   if (myFictiveConstraint)
     return; // no need several fictive constraints
 
@@ -278,12 +362,14 @@ void PlaneGCSSolver_Solver::addFictiveConstraintIfNecessary()
 void PlaneGCSSolver_Solver::removeFictiveConstraint()
 {
   if (myFictiveConstraint) {
-    myEquationSystem->removeConstraint(myFictiveConstraint);
+    myEquationSystem->clearByTag(myFictiveConstraint->getTag());
     myParameters.pop_back();
 
     GCS::VEC_pD aParams = myFictiveConstraint->params();
-    for (GCS::VEC_pD::iterator anIt = aParams.begin(); anIt != aParams.end(); ++ anIt)
-      delete *anIt;
+    for (GCS::VEC_pD::iterator anIt = aParams.begin(); anIt != aParams.end(); ++anIt) {
+      double* aPar = *anIt;
+      delete aPar;
+    }
     delete myFictiveConstraint;
     myFictiveConstraint = 0;
   }