0022851: No intersection curve between two surfaces
authorJGV <>
Fri, 13 Jan 2012 16:25:51 +0000 (16:25 +0000)
committerbugmaster <bugmaster@opencascade.com>
Mon, 5 Mar 2012 15:31:57 +0000 (19:31 +0400)
src/math/math_FunctionSetRoot.cxx

index 6866e54..0df7a4c 100755 (executable)
 //     On change de direction
 //============================================================================
 
-static  Standard_Boolean mydebug = Standard_False;
+#define FSR_DEBUG(arg)
+// Uncomment the following code to have debug output to cout 
+/* * /
+static Standard_Boolean mydebug = Standard_False;
+#undef FSR_DEBUG
+#define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }}
+/* */
 
 class MyDirFunction : public math_Function
 {
@@ -202,7 +208,7 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P0,
   F.Initialize(P1, Delta);
 
   // (2) On minimise
-  if (mydebug)cout << "      minimisation dans la direction" << endl;
+  FSR_DEBUG ("      minimisation dans la direction")
   ax = -1; bx = 0;
   cx = (P2-P1).Norm()*invnorme;
   if (cx < 1.e-2) return Standard_False;
@@ -243,7 +249,7 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
     
   // (1) On realise une premiere interpolation quadratique
   Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis;
-  if (mydebug) { cout << "     essai d interpolation" << endl;}
+  FSR_DEBUG("     essai d interpolation")
 
   df1 = Gradient*Dir;
   df2 = DGradient*Dir;
@@ -285,7 +291,7 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
   if (fsol<PValue) { 
     good = Standard_True;
     Result = fsol;
-    if (mydebug) cout << "t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue <<endl;
+    FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue)
   }
 
   // (2) Si l'on a pas assez progresser on realise une recherche 
@@ -298,14 +304,13 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
     else {
       ax = 0.0; bx = tsol; cx = 1.0;
     }
-    if (mydebug) cout << " minimisation dans la direction" << endl;
+    FSR_DEBUG(" minimisation dans la direction")
     math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
     if(Sol.IsDone()) {
       if (Sol.Minimum() <= Result) {
-       tsol = Sol.Location();
-       good = Standard_True;
-       if (mydebug) cout << "t= "<<tsol<<" F ="<< Sol.Minimum() 
-                         << " OldF = "<<Result <<endl;
+        tsol = Sol.Location();
+        good = Standard_True;
+       FSR_DEBUG("t= "<<tsol<<" F ="<< Sol.Minimum() << " OldF = "<<Result)
       }
     }
   }
@@ -336,11 +341,11 @@ static void SearchDirection(const math_Matrix& DF,
       math_Gauss Solut(DF, 1.e-9);
       if (Solut.IsDone()) Solut.Solve(Direction);
       else { // we have to "forget" singular directions.
-       if (mydebug) cout << " Matrice singuliere : On prend SVD" << endl;
-       math_SVD SolvebySVD(DF);
+       FSR_DEBUG(" Matrice singuliere : On prend SVD")
+        math_SVD SolvebySVD(DF);
         if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
-       else ChangeDirection = Standard_True;
-      } 
+        else ChangeDirection = Standard_True;
+      }
     }
     else if (Ninc > Neq) {
       math_SVD Solut(DF);
@@ -457,22 +462,28 @@ static void SearchDirection(const math_Matrix& DF,
 
 
 //====================================================
-Standard_Boolean Bounds(const math_Vector& InfBound,
-                       const math_Vector& SupBound,
-                       const math_Vector& Tol,
-                       math_Vector&       Sol,
-                       const math_Vector& SolSave,
-                       math_IntegerVector& Constraints,
-                       math_Vector&        Delta)
+Standard_Boolean Bounds(const math_Vector&  InfBound,
+                        const math_Vector&  SupBound,
+                        const math_Vector&  Tol,
+                        math_Vector&        Sol,
+                        const math_Vector&  SolSave,
+                        math_IntegerVector& Constraints,
+                        math_Vector&        Delta,
+                        Standard_Boolean&   theIsNewSol)
 //
-// Purpose : Troncate un pas d'optimisation pour rester 
-//           dans le domaine, Delta donne le pas final
+// Purpose: Trims an initial solution Sol to be within a domain defined by
+//   InfBound and SupBound. Delta will contain a distance between final Sol and
+//   SolSave.
+//   IsNewSol returns False, if final Sol fully coincides with SolSave, i.e.
+//   if SolSave already lied on a boundary and initial Sol was fully beyond it
 //======================================================
 {
   Standard_Boolean Out = Standard_False;
   Standard_Integer i, Ninc = Sol.Length();
   Standard_Real    monratio = 1;
   
+  theIsNewSol = Standard_True;
+
   // Calcul du ratio de recadrage
   for (i = 1; i <= Ninc; i++) {
     Constraints(i) = 0;
@@ -484,28 +495,36 @@ Standard_Boolean Bounds(const math_Vector& InfBound,
     else if(Sol(i) < InfBound(i)) {
       Constraints(i) = 1;
       Out = Standard_True;
-      if (Abs(Delta(i)) > Tol(i)) // Afin d'eviter des ratio nulles pour rien
-       monratio = Min(monratio, Abs( (SolSave(i)-InfBound(i))/Delta(i)) );
+      // Delta(i) is negative
+      if (-Delta(i) > Tol(i)) // Afin d'eviter des ratio nulles pour rien
+        monratio = Min(monratio, (InfBound(i) - SolSave(i))/Delta(i) );
     }
     else if (Sol(i) > SupBound(i)) {
       Constraints(i) = 1;
       Out = Standard_True;
-      if (Abs(Delta(i)) > Tol(i))
-       monratio = Min(monratio, Abs( (SolSave(i)-SupBound(i))/Delta(i)) );
+      // Delta(i) is positive
+      if (Delta(i) > Tol(i))
+        monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) );
     }
   }
 
   if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques)
-    Delta *= monratio;
-    Sol = SolSave+Delta;
-    for (i = 1; i <= Ninc; i++) {
-      if(Sol(i) < InfBound(i))  {
-       Sol(i) = InfBound(i);
-       Delta(i) = Sol(i) - SolSave(i);
-      }
-      else if (Sol(i) > SupBound(i)) {
-       Sol(i) = SupBound(i);
-       Delta(i) = Sol(i) - SolSave(i);
+    if (monratio == 0.0) {
+      theIsNewSol = Standard_False;
+      Sol = SolSave;
+      Delta.Init (0.0);
+    } else {
+      Delta *= monratio;
+      Sol = SolSave+Delta;
+      for (i = 1; i <= Ninc; i++) {
+        if(Sol(i) < InfBound(i))  {
+          Sol(i) = InfBound(i);
+          Delta(i) = Sol(i) - SolSave(i);
+        }
+        else if (Sol(i) > SupBound(i)) {
+          Sol(i) = SupBound(i);
+          Delta(i) = Sol(i) - SolSave(i);
+        }
       }
     }
   }
@@ -667,7 +686,6 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
                                   const math_Vector& InfBound,
                                   const math_Vector& SupBound) 
 {
-                                      
   Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
   
   if ((Neq <= 0)                      ||
@@ -676,10 +694,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
       (SupBound.Length() != Ninc))  { Standard_DimensionError:: Raise(); }
 
   Standard_Integer i;
-  Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False;
+  Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False;
   Standard_Boolean Good, Verif;
   Standard_Boolean Stop;
-  Standard_Real Eps = 1.e-32, Progres = 0.005;
+  const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005;
   Standard_Real F2, PreviousMinimum, Dy, OldF;
   Standard_Real Ambda, Ambda2, Gnr1, Oldgr;
   math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine
@@ -714,14 +732,12 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
   // Le rang 0 de Save ne doit servir q'au test accelarteur en fin de boucle
   // s'il on est dejas sur la solution, il faut leurer ce test pour eviter
   // de faire une seconde iteration...
-  Save(0) = Max (F2, Sqrt(Eps));
+  Save(0) = Max (F2, EpsSqrt);
 
-  if (mydebug) {
-    cout << "=== Mode Debug de Function Set Root"<<endl;
-    cout << "    F2 Initial = " << F2 << endl;
-  }
+  FSR_DEBUG("=== Mode Debug de Function Set Root" << endl)
+  FSR_DEBUG("    F2 Initial = " << F2)
 
-  if ((F2 <= Eps) || (Sqrt(Gnr1) <= Eps)) {
+  if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
     Done = Standard_True;
     State = F.GetStateNumber();
     return;
@@ -736,11 +752,14 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
     SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, DH, Dy);
     if (Abs(Dy) <= Eps) {
       Done = Standard_True;
+      ////modified by jgv, 31.08.2011////
+      F.Value(Sol, FF); //update F before GetStateNumber
+      ///////////////////////////////////
       State = F.GetStateNumber();
       return;
     }
     if (ChangeDirection) {
-      Ambda = Ambda2/Sqrt(Abs(Dy));
+      Ambda = Ambda2 / Sqrt(Abs(Dy));
       if (Ambda > 1.0) Ambda = 1.0;
     }
     else {
@@ -752,27 +771,28 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
       Sol(i) = Sol(i) + Ambda * DH(i);
     }
     Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-                 Constraints, Delta);
+                 Constraints, Delta, isNewSol);
 
       
     DHSave = GH;
-//    F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-    if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-      Done = Standard_False;
-      State = F.GetStateNumber();
-      return;
-    }
-
-#if DEB
-    if (mydebug) {
-      cout << "Kount         = " << Kount << endl;
-      cout << "Le premier F2 = " << F2 << endl;
-      cout << "Direction     = " << ChangeDirection << endl;
+    if (isNewSol) {
+//      F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+      if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+        Done = Standard_False;
+        State = F.GetStateNumber();
+        return;
+      }
     }
-#endif
+    
+    FSR_DEBUG("Kount         = " << Kount)
+    FSR_DEBUG("Le premier F2 = " << F2)
+    FSR_DEBUG("Direction     = " << ChangeDirection)
 
-    if ((F2 <= Eps) || (Sqrt(Gnr1) <= Eps)) {
+    if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
       Done = Standard_True;
+      ////modified by jgv, 31.08.2011////
+      F.Value(Sol, FF); //update F before GetStateNumber
+      ///////////////////////////////////
       State = F.GetStateNumber();
       return;
     }
@@ -792,16 +812,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
        while((F2/PreviousMinimum > Progres) && !Stop) {
          if (F2 < OldF && (Dy < 0.0)) {
            // On essaye de progresser dans cette direction.
-           if (mydebug)  cout << " iteration de descente = " << DescenteIter<<endl;
+           FSR_DEBUG(" iteration de descente = " << DescenteIter)
            DescenteIter++;
            SolSave = Sol;
            OldF = F2;
            for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
-           Sol(i) = Sol(i) + Ambda * DH(i);
-         }
+             Sol(i) = Sol(i) + Ambda * DH(i);
+           }
            Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave, 
-                         Constraints, Delta);
-           if (mydebug) { cout << " Augmentation de lambda" << endl;}
+                         Constraints, Delta, isNewSol);
+           FSR_DEBUG(" Augmentation de lambda")
            Ambda *= 1.7;
          }
          else {
@@ -826,36 +846,40 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
                F2 = OldF;
              }
              else {
-               Sol = SolSave+Delta;
+               Sol = SolSave+Delta;
+               Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
+                 Constraints, Delta, isNewSol);
              }
-             Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-                           Constraints, Delta);
              Sort = Standard_False; // On a rejete le point sur la frontiere
            }
            Stop = Standard_True; // et on sort dans tous les cas...
          }
          DHSave = GH;
-//       F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-         if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-           Done = Standard_False;
-           State = F.GetStateNumber();
-           return;
-         }
+          if (isNewSol) {
+//            F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+            if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+              Done = Standard_False;
+              State = F.GetStateNumber();
+              return;
+            }
+          }
          Dy = GH*DH;
          if (Abs(Dy) <= Eps) {
-           State = F.GetStateNumber();
            Done = Standard_True;
-           if (F2 > OldF) Sol = SolSave;
+           if (F2 > OldF)
+              Sol = SolSave;
+            ////modified by jgv, 31.08.2011////
+            F.Value(Sol, FF); //update F before GetStateNumber
+            ///////////////////////////////////
+           State = F.GetStateNumber();
            return;
          }
          if (DescenteIter >= 10) {
            Stop = Standard_True;
          }
        }
-       if (mydebug) {
-         cout << "--- Sortie du Traitement Standard"<<endl;
-         cout << "    DescenteIter = "<<DescenteIter << " F2 = " << F2 << endl;
-       }
+       FSR_DEBUG("--- Sortie du Traitement Standard")
+       FSR_DEBUG("    DescenteIter = "<<DescenteIter << " F2 = " << F2)
       }
       // ------------------------------------
       //  on passe au traitement des bords
@@ -872,12 +896,12 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
          OldF = F2;
          SearchDirection(DF, GH, FF,  Constraints, Sol, 
                          ChangeDirection, InvLengthMax, DH, Dy);
-         if (mydebug) { cout << " Conditional Direction = " << ChangeDirection << endl;}
+         FSR_DEBUG(" Conditional Direction = " << ChangeDirection)
          if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
            if (ChangeDirection) {
 
-//           Ambda = Ambda2/Sqrt(Abs(Dy));
-             Ambda = Ambda2/Sqrt(-Dy);
+//           Ambda = Ambda2 / Sqrt(Abs(Dy));
+             Ambda = Ambda2 / Sqrt(-Dy);
              if (Ambda > 1.0) Ambda = 1.0;
            }
            else {
@@ -889,20 +913,20 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
              Sol(i) = Sol(i) + Ambda * DH(i);
            }
            Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-                            Constraints, Delta);
+                            Constraints, Delta, isNewSol);
 
            DHSave = GH;
-//         F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-           if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-             Done = Standard_False;
-             State = F.GetStateNumber();
-             return;
-           }
+            if (isNewSol) {
+//              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+              if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+                Done = Standard_False;
+                State = F.GetStateNumber();
+                return;
+              }
+            }
            Ambda2 = Gnr1;
-            if (mydebug) {
-              cout << "---  Iteration au bords : " << DescenteIter << endl; 
-              cout << "---  F2 = " << F2 << endl;
-            }
+           FSR_DEBUG("---  Iteration au bords : " << DescenteIter)
+            FSR_DEBUG("---  F2 = " << F2)
          }
          else {
            Stop = Standard_True;
@@ -910,8 +934,7 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 
          while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
            DescenteIter++;
-           if (mydebug)  cout << "--- Iteration de descente conditionnel = " 
-             << DescenteIter<<endl;
+           FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter)
            if (F2 < OldF && Dy < 0.0) {
              // On essaye de progresser dans cette direction.
              SolSave = Sol;
@@ -920,15 +943,17 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
                Sol(i) = Sol(i) + Ambda * DH(i);
              }
              Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-                              Constraints, Delta );     
+                              Constraints, Delta, isNewSol);     
            }
            DHSave = GH;
-//         F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-           if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-             Done = Standard_False;
-             State = F.GetStateNumber();
-             return;
-           }
+            if (isNewSol) {
+//              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+              if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+                Done = Standard_False;
+                State = F.GetStateNumber();
+                return;
+              }
+            }
            Ambda2 = Gnr1;
            Dy = GH*DH;
            Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
@@ -943,24 +968,25 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
                                   DHSave, GH, Tol, F_Dir);
          if (!Good) {
            Sol = SolSave;
+            Sort = Standard_False;
          }
          else {
            Sol = SolSave + Delta;
-         }
-         Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-                       Constraints, Delta);
-//       F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-         if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-           Done = Standard_False;
-           State = F.GetStateNumber();
-           return;
+           Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
+             Constraints, Delta, isNewSol);
+            if (isNewSol) {
+//              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+              if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+                Done = Standard_False;
+                State = F.GetStateNumber();
+                return;
+              }
+           }
          }
          Dy = GH*DH;
        }       
-       if (mydebug) {
-         cout << "--- Sortie du Traitement des Bords"<<endl;
-         cout << "--- DescenteIter = "<<DescenteIter << " F2 = " << F2 << endl;
-       }      
+       FSR_DEBUG("--- Sortie du Traitement des Bords")
+       FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2)
       }
     }
 
@@ -982,12 +1008,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
       for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
        Delta(i) = PreviousSolution(i) - Sol(i);
       }
+      
       if (IsSolutionReached(F)) {
        if (PreviousMinimum < F2) {
          Sol = SolSave;
        }
-       State = F.GetStateNumber();
        Done = Standard_True;
+        ////modified by jgv, 31.08.2011////
+        F.Value(Sol, FF); //update F before GetStateNumber
+        ///////////////////////////////////
+       State = F.GetStateNumber();
        return;
       }
     }
@@ -1021,6 +1051,9 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
        }
        if (IsSolutionReached(F)) {
          Done = Standard_True;
+          ////modified by jgv, 31.08.2011////
+          F.Value(Sol, FF); //update F before GetStateNumber
+          ///////////////////////////////////
          State = F.GetStateNumber();
          return;
        }
@@ -1078,10 +1111,3 @@ void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const{
   Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " ");
   Err = Delta;
 }
-
-
-
-
-
-
-