0023830: BRepExtrema_DistShapeShape does not find intersection of face with edge

Data member myIsDivergent was added to class math_FunctionSetRoot.
The member is initialized in method Perform start by Standard_False. The member is changed in the method to Standard_True if an approximation point is located outside of search set.

Method IsDivergent was added to class math_FunctionSetRoot.
The method returns value of myIsDivergent.

Parameter theStopOnDivergent with default value Standard_False was added to
a constructor and method Perform of class math_FunctionSetRoot.
The parameter shows whether return from method Perform if myIsDivergent
became Standard_True.

Class Extrema_GenExtCS was optimized for dynamic memory consumption:
the class surface values are not stored now.

Method Perform of class Extrema_GenExtCS was improved for local extreme
search in case of iteration algorithm divergence. Size of the algorithm
grid cell for initial approximation is reduced by 2 in the case. The count
of the reduces is limited.
Build fix in *.cdl.
Case of extrusion was carried out of the repeat loop.

Division method of the 3-dimensional grid (curve param, surface first
param, surface second param) by two was replaced by method of moving the
grid.

The moving method is about two times faster than division method.
Now grid moving method is depending on some average linear sizes of the curve and the surface.
In case (theStopOnDivergent && myIsDivergent) method math_FunctionSetRoot::Perform sets Done to Standard_False and discard it's last solution.
Test case for the bug fix.
Test case correction
Initial approximation was corrected for search of maximum in Extrema_GenExtCS::Perform.
Test case bugs moddata_3 bug23830 was formatted.

diff --git a/src/Extrema/Extrema_GenExtCS.cxx b/src/Extrema/Extrema_GenExtCS.cxx
index 33261d1..3b2639c 100755 (executable)
--- a/src/Extrema/Extrema_GenExtCS.cxx
+++ b/src/Extrema/Extrema_GenExtCS.cxx
@@ -117,7 +117,6 @@ void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S,
                                  const Standard_Real Tol2)
 {
   myS = (Adaptor3d_SurfacePtr)&S;
                                  const Standard_Real Tol2)
 {
   myS = (Adaptor3d_SurfacePtr)&S;

-  mypoints2 = new TColgp_HArray2OfPnt(0,NbU+1,0,NbV+1);

   myusample = NbU;
   myvsample = NbV;
   myumin = Umin;
   myusample = NbU;
   myvsample = NbV;
   myumin = Umin;


@@ -125,29 +124,6 @@ void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S,
   myvmin = Vmin;
   myvsup = Vsup;
   mytol2 = Tol2;
   myvmin = Vmin;
   myvsup = Vsup;
   mytol2 = Tol2;

-
-// Parametrage de l echantillon sur S
-
-  Standard_Real PasU = myusup - myumin;
-  Standard_Real PasV = myvsup - myvmin;
-  Standard_Real U0 = PasU / myusample / 100.;
-  Standard_Real V0 = PasV / myvsample / 100.;
-  gp_Pnt P1;
-  PasU = (PasU - U0) / (myusample - 1);
-  PasV = (PasV - V0) / (myvsample - 1);
-  U0 = myumin + U0/2.;
-  V0 = myvmin + V0/2.;
-
-// Calcul des distances
-
-  Standard_Integer NoU, NoV;
-  Standard_Real U, V;
-  for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
-    for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
-      P1 = myS->Value(U, V);
-      mypoints2->SetValue(NoU,NoV,P1);
-    }
-  }

 }
 
 //=======================================================================
 }
 
 //=======================================================================


@@ -181,13 +157,7 @@ void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
   mytsup = tsup;
   mytol1 = Tol1;
   mytsample = NbT;
   mytsup = tsup;
   mytol1 = Tol1;
   mytsample = NbT;

-  mypoints1 = new TColgp_HArray1OfPnt(0,NbT);
-
-  Standard_Real t1, U, V;
-  Standard_Integer NoT, NoU, NoV;
-  gp_Pnt P1, P2;
-
-// Modif de lvt pour trimer la surface non pas aux infinis mais  a +/- 10000
+  // Modif de lvt pour trimer la surface non pas aux infinis mais  a +/- 10000

 
   Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin;
   if (Precision::IsInfinite(trimusup)){
 
   Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin;
   if (Precision::IsInfinite(trimusup)){


@@ -202,40 +172,7 @@ void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
   if (Precision::IsInfinite(trimvmin)){
     trimvmin = -1.0e+10;
   }
   if (Precision::IsInfinite(trimvmin)){
     trimvmin = -1.0e+10;
   }

-  
-
-// Parametrage de l echantillon sur C
-
-  Standard_Real PasT = mytsup - mytmin;
-  Standard_Real t0 = PasT / mytsample / 100.;
-  PasT = (PasT - t0) / (mytsample - 1);
-  t0 = mytmin + t0/2.;
-
-// Parametrage de l echantillon sur S
-
-  Standard_Real PasU = trimusup - trimumin;
-  Standard_Real PasV = trimvsup - trimvmin;
-  Standard_Real U0 = PasU / myusample / 100.;
-  Standard_Real V0 = PasV / myvsample / 100.;
-  PasU = (PasU - U0) / (myusample - 1);
-  PasV = (PasV - V0) / (myvsample - 1);
-  U0 = trimumin + U0/2.;
-  V0 = trimvmin + V0/2.;
-
-// Calcul des distances
-
-  for ( NoT = 1, t1 = t0; NoT <= mytsample; NoT++, t1 += PasT) {
-    P1 = C.Value(t1);
-    mypoints1->SetValue(NoT,P1);
-  }
-  
-
-/*
-b- Calcul des minima:
-   -----------------
-   b.a) Initialisations:
-*/
-
+  //

   math_Vector Tol(1, 3);
   Tol(1) = mytol1;
   Tol(2) = mytol2;
   math_Vector Tol(1, 3);
   Tol(1) = mytol1;
   Tol(2) = mytol2;


@@ -248,7 +185,6 @@ b- Calcul des minima:
   UVsup(1) = mytsup;
   UVsup(2) = trimusup;
   UVsup(3) = trimvsup;
   UVsup(1) = mytsup;
   UVsup(2) = trimusup;
   UVsup(3) = trimvsup;

-

   // 18/02/02 akm vvv : (OCC163) bad extremas - extrusion surface versus the line.
   //                    Try to preset the initial solution as extrema between
   //                    extrusion direction and the curve.
   // 18/02/02 akm vvv : (OCC163) bad extremas - extrusion surface versus the line.
   //                    Try to preset the initial solution as extrema between
   //                    extrusion direction and the curve.


@@ -268,80 +204,168 @@ b- Calcul des minima:
       Extrema_POnCurv aP1, aP2;
       for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
       {
       Extrema_POnCurv aP1, aP2;
       for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
       {

-       aLocator.Points (iExt, aP1, aP2);
-       // Parameter on curve
-       UV(1) = aP1.Parameter();
-       // To find parameters on surf, try ExtPS
-       Extrema_ExtPS aPreciser (aP1.Value(), *myS, mytol2, mytol2);
-       if (aPreciser.IsDone())
-       {
-         // Managed to find extremas between point and surface
-         Standard_Integer iPExt;
-         for (iPExt=1; iPExt<=aPreciser.NbExt(); iPExt++)
-         {
-           aPreciser.Point(iPExt).Parameter(UV(2),UV(3));
-           math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
-         }
-       }
-       else
-       {
-         // Failed... try the point on iso line
-         UV(2) = dfUFirst;
-         UV(3) = aP2.Parameter();
-         math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
-       }
+  aLocator.Points (iExt, aP1, aP2);
+  // Parameter on curve
+  UV(1) = aP1.Parameter();
+  // To find parameters on surf, try ExtPS
+  Extrema_ExtPS aPreciser (aP1.Value(), *myS, mytol2, mytol2);
+  if (aPreciser.IsDone())
+  {
+    // Managed to find extremas between point and surface
+    Standard_Integer iPExt;
+    for (iPExt=1; iPExt<=aPreciser.NbExt(); iPExt++)
+    {
+      aPreciser.Point(iPExt).Parameter(UV(2),UV(3));
+      math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
+    }
+  }
+  else
+  {
+    // Failed... try the point on iso line
+    UV(2) = dfUFirst;
+    UV(3) = aP2.Parameter();
+    math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
+  }

       } // for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
     } // if (aLocator.IsDone() && aLocator.NbExt()>0)
   } // if (myS.Type() == GeomAbs_ExtrusionSurface)
       } // for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
     } // if (aLocator.IsDone() && aLocator.NbExt()>0)
   } // if (myS.Type() == GeomAbs_ExtrusionSurface)

-  else 
-  // 18/02/02 akm ^^^
+  else

   {
   {

-    Standard_Real distmin = RealLast(), distmax = 0.0, TheDist;
-
-    Standard_Integer NTmin=0,NUmin=0,NVmin=0;
-    gp_Pnt PP1min, PP2min;
-    Standard_Integer NTmax=0,NUmax=0,NVmax=0;
-    gp_Pnt PP1max, PP2max;
-    for ( NoT = 1, t1 = t0; NoT <= mytsample; NoT++, t1 += PasT) {
-      P1 = mypoints1->Value(NoT);
-      for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
-       for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
-         P2 = mypoints2->Value(NoU, NoV);
-         TheDist = P1.SquareDistance(P2);
-         if (TheDist < distmin) {
-           distmin = TheDist;
-           NTmin = NoT;
-           NUmin = NoU;
-           NVmin = NoV;
-           PP1min = P1;
-           PP2min = P2;
-         }
-         if (TheDist > distmax) {
-           distmax = TheDist;
-           NTmax = NoT;
-           NUmax = NoU;
-           NVmax = NoV;
-           PP1max = P1;
-           PP2max = P2;
-         }
-       }
+    Standard_Real aCUAdd = (mytsup - mytmin) / mytsample;
+    Standard_Real aSUAdd = (myusup - myumin) / myusample;
+    Standard_Real aSVAdd = (myvsup - myvmin) / myvsample;
+    TColgp_HArray1OfPnt aCPs(1, mytsample);
+    TColgp_HArray2OfPnt aSPs(1, myusample, 1, myvsample);
+    Standard_Integer aRestIterCount = 3;
+      // The value is calculated by the bug CR23830.
+    Standard_Integer aCUDen = 2, aSUDen = 2, aSVDen = 2;
+    Standard_Boolean anAreAvSqsInited = Standard_False;
+    Standard_Real aCUSq = 0, aSUSq = 0, aSVSq = 0;
+    while (aRestIterCount--)
+    {
+      Standard_Real aMinCU, aMinSU, aMinSV, aMaxCU, aMaxSU, aMaxSV;
+      Standard_Real aMinSqDist = DBL_MAX, aMaxSqDist = -DBL_MAX;
+      for (Standard_Integer aSUNom = 1; aSUNom < aSUDen; aSUNom += 2)
+      {
+        Standard_Real aSU0 = myumin + (aSUNom * aSUAdd) / aSUDen;
+        for (Standard_Integer aSVNom = 1; aSVNom < aSVDen; aSVNom += 2)
+        {
+          Standard_Real aSV0 = myvmin + (aSVNom * aSVAdd) / aSVDen;
+          for (Standard_Integer aCUNom = 1; aCUNom < aCUDen; aCUNom += 2)
+          {
+            Standard_Real aCU0 = mytmin + (aCUNom * aCUAdd) / aCUDen;
+            Standard_Real aCU = aCU0;
+            for (Standard_Integer aCUI = 1; aCUI <= mytsample;
+              aCUI++, aCU += aCUAdd)
+            {
+              aCPs.SetValue(aCUI, C.Value(aCU));
+            }
+            //
+            aCU = aCU0;
+            Standard_Real aSU = aSU0;
+            for (Standard_Integer aSUI = 1; aSUI <= myusample;
+              aSUI++, aSU += aSUAdd)
+            {
+              Standard_Real aSV = aSV0;
+              for (Standard_Integer aSVI = 1; aSVI <= myvsample;
+                aSVI++, aSV += aSVAdd)
+              {
+                gp_Pnt aSP = myS->Value(aSU, aSV);
+                aSPs.ChangeValue(aSUI, aSVI) = aSP;
+                Standard_Real aCU = aCU0;
+                for (Standard_Integer aCUI = 1; aCUI <= mytsample;
+                  aCUI++, aCU += aCUAdd)
+                {
+                  Standard_Real aSqDist = aSP.SquareDistance(aCPs.Value(aCUI));
+                  if (aSqDist < aMinSqDist)
+                  {
+                    aMinSqDist = aSqDist;
+                    aMinCU = aCU;
+                    aMinSU = aSU;
+                    aMinSV = aSV;
+                  }
+                  if (aMaxSqDist < aSqDist)
+                  {
+                    aMaxSqDist = aSqDist;
+                    aMaxCU = aCU;
+                    aMaxSU = aSU;
+                    aMaxSV = aSV;
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+      // Find min approximation.
+      UV(1) = aMinCU;
+      UV(2) = aMinSU;
+      UV(3) = aMinSV;
+      math_FunctionSetRoot anA(myF, UV, Tol, UVinf, UVsup, 100, aRestIterCount);
+      // Find max approximation.
+      if (!anA.IsDivergent() || !aRestIterCount)
+      {
+        UV(1) = aMaxCU;
+        UV(2) = aMaxSU;
+        UV(3) = aMaxSV;
+        math_FunctionSetRoot(myF, UV, Tol, UVinf, UVsup);
+        break;
+      }
+      //
+      if (!anAreAvSqsInited)
+      {
+        anAreAvSqsInited = Standard_True;
+        //
+        const gp_Pnt * aCP1 = &aCPs.Value(1);
+        for (Standard_Integer aCUI = 2; aCUI <= mytsample; aCUI++)
+        {
+          const gp_Pnt & aCP2 = aCPs.Value(aCUI);
+          aCUSq += aCP1->SquareDistance(aCP2);
+          aCP1 = &aCP2;
+        }
+        aCUSq /= mytsample - 1;
+        //
+        for (Standard_Integer aSVI = 1; aSVI <= myvsample; aSVI++)
+        {
+          const gp_Pnt * aSP1 = &aSPs.Value(1, aSVI);
+          for (Standard_Integer aSUI = 2; aSUI <= myusample; aSUI++)
+          {
+            const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
+            aSUSq += aSP1->SquareDistance(aSP2);
+            aSP1 = &aSP2;
+          }
+        }
+        aSUSq /= myvsample * (myusample - 1);
+        //
+        for (Standard_Integer aSUI = 1; aSUI <= myusample; aSUI++)
+        {
+          const gp_Pnt * aSP1 = &aSPs.Value(aSUI, 1);
+          for (Standard_Integer aSVI = 2; aSVI <= myvsample; aSVI++)
+          {
+            const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
+            aSVSq += aSP1->SquareDistance(aSP2);
+            aSP1 = &aSP2;
+          }
+        }
+        aSVSq /= myusample * (myvsample - 1);
+      }
+      //
+      if ((aSUSq <= aCUSq) && (aSVSq <= aCUSq))
+      {
+        aCUDen += aCUDen;
+        aCUSq /= 4;
+      }
+      else if ((aCUSq <= aSUSq) && (aSVSq <= aSUSq))
+      {
+        aSUDen += aSUDen;
+        aSUSq /= 4;
+      }
+      else
+      {
+        aSVDen += aSVDen;
+        aSVSq /= 4;

       }
     }
       }
     }

-  
-    // minimum:
-    UV(1) = t0 + (NTmin - 1) * PasT;
-    UV(2) = U0 + (NUmin - 1) * PasU;
-    UV(3) = V0 + (NVmin - 1) * PasV;
-
-    math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
-
-    // maximum:
-    UV(1) = t0 + (NTmax - 1) * PasT;
-    UV(2) = U0 + (NUmax - 1) * PasU;
-    UV(3) = V0 + (NVmax - 1) * PasV;
-
-    math_FunctionSetRoot S2 (myF,UV,Tol,UVinf,UVsup);
-

   }
 
   myDone = Standard_True;
   }
 
   myDone = Standard_True;

diff --git a/src/math/math_FunctionSetRoot.cdl b/src/math/math_FunctionSetRoot.cdl
index 721c212..e842d6b 100755 (executable)
--- a/src/math/math_FunctionSetRoot.cdl
+++ b/src/math/math_FunctionSetRoot.cdl
@@ -82,7 +82,7 @@ is
 
     Create(F: in out FunctionSetWithDerivatives; StartingPoint: Vector;
           Tolerance: Vector; infBound, supBound: Vector;
 
     Create(F: in out FunctionSetWithDerivatives; StartingPoint: Vector;
           Tolerance: Vector; infBound, supBound: Vector;

-          NbIterations: Integer = 100)
+          NbIterations: Integer = 100; theStopOnDivergent : Boolean from Standard = Standard_False)

        ---Purpose: is used to improve the root of the function F
        --          from the initial guess StartingPoint.
        --          The maximum number of iterations allowed is given 
        ---Purpose: is used to improve the root of the function F
        --          from the initial guess StartingPoint.
        --          The maximum number of iterations allowed is given 


@@ -105,8 +105,8 @@ is
     
     
     Perform(me: in out; F: in out FunctionSetWithDerivatives;
     
     
     Perform(me: in out; F: in out FunctionSetWithDerivatives;

-           StartingPoint: Vector;
-           infBound, supBound: Vector)
+        StartingPoint: Vector;
+        infBound, supBound: Vector; theStopOnDivergent : Boolean from Standard = Standard_False)

        ---Purpose: Improves the root of function F from the initial guess
        -- StartingPoint. infBound and supBound may be given to constrain the solution.
        -- Warning
        ---Purpose: Improves the root of function F from the initial guess
        -- StartingPoint. infBound and supBound may be given to constrain the solution.
        -- Warning


@@ -231,6 +231,10 @@ is
 
     is static;
 
 
     is static;
 

+    IsDivergent(me)
+      returns Boolean
+      is static;
+

     
 fields
 
     
 fields
 


@@ -260,5 +264,6 @@ Temp1            : Vector;
 Temp2            : Vector;
 Temp3            : Vector;
 Temp4            : Vector;
 Temp2            : Vector;
 Temp3            : Vector;
 Temp4            : Vector;

+myIsDivergent    : Boolean from Standard;

 
 end FunctionSetRoot;
 
 end FunctionSetRoot;

diff --git a/src/math/math_FunctionSetRoot.cxx b/src/math/math_FunctionSetRoot.cxx
index 9f338b3..e6da643 100755 (executable)
--- a/src/math/math_FunctionSetRoot.cxx
+++ b/src/math/math_FunctionSetRoot.cxx
@@ -616,7 +616,8 @@ math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
                                           const math_Vector& Tolerance,
                                           const math_Vector& infBound,
                                           const math_Vector& supBound,
                                           const math_Vector& Tolerance,
                                           const math_Vector& infBound,
                                           const math_Vector& supBound,

-                                          const Standard_Integer NbIterations) :
+                                          const Standard_Integer NbIterations,
+                                          Standard_Boolean theStopOnDivergent) :

                                            Delta(1, F.NbVariables()),
                                            Sol(1, F.NbVariables()),
                                            DF(1, F.NbEquations(),
                                            Delta(1, F.NbVariables()),
                                            Sol(1, F.NbVariables()),
                                            DF(1, F.NbEquations(),


@@ -646,7 +647,7 @@ math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
     Tol(i) =Tolerance(i);
   }
   Itermax = NbIterations;
     Tol(i) =Tolerance(i);
   }
   Itermax = NbIterations;

-  Perform(F, StartingPoint, infBound, supBound);
+  Perform(F, StartingPoint, infBound, supBound, theStopOnDivergent);

 }
 
 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
 }
 
 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,


@@ -698,7 +699,8 @@ void math_FunctionSetRoot::SetTolerance(const math_Vector& Tolerance)
 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
                                   const math_Vector& StartingPoint,
                                   const math_Vector& InfBound,
 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
                                   const math_Vector& StartingPoint,
                                   const math_Vector& InfBound,

-                                  const math_Vector& SupBound) 
+                                  const math_Vector& SupBound,
+                                  Standard_Boolean theStopOnDivergent)

 {
   Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
   
 {
   Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
   


@@ -729,6 +731,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
   Sol = StartingPoint;
   Kount = 0;
 
   Sol = StartingPoint;
   Kount = 0;
 

+  //
+  myIsDivergent = Standard_False;
+  for (i = 1; i <= Ninc; i++)
+  {
+    myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+  }
+  if (theStopOnDivergent & myIsDivergent)
+  {
+    return;
+  }

   // Verification de la validite des inconnues par rapport aux bornes.
   // Recentrage sur les bornes si pas valide.
   for ( i = 1; i <= Ninc; i++) {
   // Verification de la validite des inconnues par rapport aux bornes.
   // Recentrage sur les bornes si pas valide.
   for ( i = 1; i <= Ninc; i++) {


@@ -739,7 +751,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
   // Calcul de la premiere valeur de F et de son gradient
   if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
     Done = Standard_False;
   // Calcul de la premiere valeur de F et de son gradient
   if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
     Done = Standard_False;

-    State = F.GetStateNumber();
+    if (!theStopOnDivergent || !myIsDivergent)
+    {
+      State = F.GetStateNumber();
+    }

     return;
   }
   Ambda2 = Gnr1;
     return;
   }
   Ambda2 = Gnr1;


@@ -752,8 +767,12 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
   FSR_DEBUG("    F2 Initial = " << F2)
 
   if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
   FSR_DEBUG("    F2 Initial = " << F2)
 
   if ((F2 <= Eps) || (Gnr1 <= Eps2)) {

-    Done = Standard_True;
-    State = F.GetStateNumber();
+    Done = Standard_False;
+    if (!theStopOnDivergent || !myIsDivergent)
+    {
+      Done = Standard_True;
+      State = F.GetStateNumber();
+    }

     return;
   }
 
     return;
   }
 


@@ -765,11 +784,15 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 
     SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, DH, Dy);
     if (Abs(Dy) <= Eps) {
 
     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();
+      Done = Standard_False;
+      if (!theStopOnDivergent || !myIsDivergent)
+      {
+        Done = Standard_True;
+        ////modified by jgv, 31.08.2011////
+        F.Value(Sol, FF); //update F before GetStateNumber
+        ///////////////////////////////////
+        State = F.GetStateNumber();
+      }

       return;
     }
     if (ChangeDirection) {
       return;
     }
     if (ChangeDirection) {


@@ -784,6 +807,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
     for( i = Sol.Lower(); i <= Sol.Upper(); i++) { 
       Sol(i) = Sol(i) + Ambda * DH(i);
     }
     for( i = Sol.Lower(); i <= Sol.Upper(); i++) { 
       Sol(i) = Sol(i) + Ambda * DH(i);
     }

+    //
+    for (i = 1; i <= Ninc; i++)
+    {
+      myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+    }
+    if (theStopOnDivergent & myIsDivergent)
+    {
+      return;
+    }
+    //

     Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                  Constraints, Delta, isNewSol);
 
     Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                  Constraints, Delta, isNewSol);
 


@@ -793,7 +826,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 //      F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
       if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
         Done = Standard_False;
 //      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();
+        if (!theStopOnDivergent || !myIsDivergent)
+        {
+          State = F.GetStateNumber();
+        }

         return;
       }
     }
         return;
       }
     }


@@ -803,11 +839,15 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
     FSR_DEBUG("Direction     = " << ChangeDirection)
 
     if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
     FSR_DEBUG("Direction     = " << ChangeDirection)
 
     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();
+      Done = Standard_False;
+      if (!theStopOnDivergent || !myIsDivergent)
+      {
+        Done = Standard_True;
+        ////modified by jgv, 31.08.2011////
+        F.Value(Sol, FF); //update F before GetStateNumber
+        ///////////////////////////////////
+        State = F.GetStateNumber();
+      }

       return;
     }
 
       return;
     }
 


@@ -833,6 +873,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
            for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
              Sol(i) = Sol(i) + Ambda * DH(i);
            }
            for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
              Sol(i) = Sol(i) + Ambda * DH(i);
            }

+      //
+      for (i = 1; i <= Ninc; i++)
+      {
+        myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+      }
+      if (theStopOnDivergent & myIsDivergent)
+      {
+        return;
+      }
+      //

            Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave, 
                          Constraints, Delta, isNewSol);
            FSR_DEBUG(" Augmentation de lambda")
            Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave, 
                          Constraints, Delta, isNewSol);
            FSR_DEBUG(" Augmentation de lambda")


@@ -861,6 +911,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
              }
              else {
                Sol = SolSave+Delta;
              }
              else {
                Sol = SolSave+Delta;

+          //
+          for (i = 1; i <= Ninc; i++)
+          {
+            myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+          }
+          if (theStopOnDivergent & myIsDivergent)
+          {
+            return;
+          }
+          //

                Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                  Constraints, Delta, isNewSol);
              }
                Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                  Constraints, Delta, isNewSol);
              }


@@ -873,19 +933,26 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 //            F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
             if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
               Done = Standard_False;
 //            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();
+              if (!theStopOnDivergent || !myIsDivergent)
+              {
+                State = F.GetStateNumber();
+              }

               return;
             }
           }
          Dy = GH*DH;
          if (Abs(Dy) <= Eps) {
               return;
             }
           }
          Dy = GH*DH;
          if (Abs(Dy) <= Eps) {

-           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();
+      Done = Standard_False;
+      if (!theStopOnDivergent || !myIsDivergent)
+      {
+        Done = Standard_True;
+              ////modified by jgv, 31.08.2011////
+              F.Value(Sol, FF); //update F before GetStateNumber
+              ///////////////////////////////////
+        State = F.GetStateNumber();
+      }

            return;
          }
          if (DescenteIter >= 10) {
            return;
          }
          if (DescenteIter >= 10) {


@@ -926,6 +993,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
            for( i = Sol.Lower(); i <= Sol.Upper(); i++) { 
              Sol(i) = Sol(i) + Ambda * DH(i);
            }
            for( i = Sol.Lower(); i <= Sol.Upper(); i++) { 
              Sol(i) = Sol(i) + Ambda * DH(i);
            }

+      //
+      for (i = 1; i <= Ninc; i++)
+      {
+        myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+      }
+      if (theStopOnDivergent & myIsDivergent)
+      {
+        return;
+      }
+      //

            Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                             Constraints, Delta, isNewSol);
 
            Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                             Constraints, Delta, isNewSol);
 


@@ -934,7 +1011,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
               if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
                 Done = Standard_False;
 //              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();
+                if (!theStopOnDivergent || !myIsDivergent)
+                {
+                  State = F.GetStateNumber();
+                }

                 return;
               }
             }
                 return;
               }
             }


@@ -956,6 +1036,16 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
              for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
                Sol(i) = Sol(i) + Ambda * DH(i);
              }
              for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
                Sol(i) = Sol(i) + Ambda * DH(i);
              }

+        //
+        for (i = 1; i <= Ninc; i++)
+        {
+          myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+        }
+        if (theStopOnDivergent & myIsDivergent)
+        {
+          return;
+        }
+        //

              Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                               Constraints, Delta, isNewSol);     
            }
              Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
                               Constraints, Delta, isNewSol);     
            }


@@ -964,7 +1054,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
               if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
                 Done = Standard_False;
 //              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();
+                if (!theStopOnDivergent || !myIsDivergent)
+                {
+                  State = F.GetStateNumber();
+                }

                 return;
               }
             }
                 return;
               }
             }


@@ -986,13 +1079,26 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
          }
          else {
            Sol = SolSave + Delta;
          }
          else {
            Sol = SolSave + Delta;

+      //
+      for (i = 1; i <= Ninc; i++)
+      {
+        myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+      }
+      if (theStopOnDivergent & myIsDivergent)
+      {
+        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;
            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();
+                if (!theStopOnDivergent || !myIsDivergent)
+                {
+                  State = F.GetStateNumber();
+                }

                 return;
               }
            }
                 return;
               }
            }


@@ -1027,11 +1133,15 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
        if (PreviousMinimum < F2) {
          Sol = SolSave;
        }
        if (PreviousMinimum < F2) {
          Sol = SolSave;
        }

-       Done = Standard_True;
-        ////modified by jgv, 31.08.2011////
-        F.Value(Sol, FF); //update F before GetStateNumber
-        ///////////////////////////////////
-       State = F.GetStateNumber();
+  Done = Standard_False;
+  if (!theStopOnDivergent || !myIsDivergent)
+  {
+    Done = Standard_True;
+          ////modified by jgv, 31.08.2011////
+          F.Value(Sol, FF); //update F before GetStateNumber
+          ///////////////////////////////////
+    State = F.GetStateNumber();
+  }

        return;
       }
     }
        return;
       }
     }


@@ -1046,8 +1156,12 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
          if (!ChangeDirection) ChangeDirection = Standard_True;
          else 
     {
          if (!ChangeDirection) ChangeDirection = Standard_True;
          else 
     {

-      Done = Standard_True;
-      State = F.GetStateNumber();
+      Done = Standard_False;
+      if (!theStopOnDivergent || !myIsDivergent)
+      {
+        Done = Standard_True;
+        State = F.GetStateNumber();
+      }

       return; //  si un gain inf a 5% on sort
     }
        }
       return; //  si un gain inf a 5% on sort
     }
        }


@@ -1070,11 +1184,15 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
          Delta(i) = PreviousSolution(i) - Sol(i);
        }
        if (IsSolutionReached(F)) {
          Delta(i) = PreviousSolution(i) - Sol(i);
        }
        if (IsSolutionReached(F)) {

-         Done = Standard_True;
-          ////modified by jgv, 31.08.2011////
-          F.Value(Sol, FF); //update F before GetStateNumber
-          ///////////////////////////////////
-         State = F.GetStateNumber();
+    Done = Standard_False;
+    if (!theStopOnDivergent || !myIsDivergent)
+    {
+      Done = Standard_True;
+            ////modified by jgv, 31.08.2011////
+            F.Value(Sol, FF); //update F before GetStateNumber
+            ///////////////////////////////////
+      State = F.GetStateNumber();
+    }

          return;
        }
       }
          return;
        }
       }


@@ -1086,19 +1204,28 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 //     F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
        if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
          Done = Standard_False;
 //     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();
+    if (!theStopOnDivergent || !myIsDivergent)
+    {
+      State = F.GetStateNumber();
+    }

          return;
        }
       }
       else 
       {
         
          return;
        }
       }
       else 
       {
         

-        State = F.GetStateNumber();
+        if (!theStopOnDivergent || !myIsDivergent)
+        {
+          State = F.GetStateNumber();
+        }

         return; // y a plus d'issues
       }
     }
   }
         return; // y a plus d'issues
       }
     }
   }

-  State = F.GetStateNumber();
+  if (!theStopOnDivergent || !myIsDivergent)
+  {
+    State = F.GetStateNumber();
+  }

 }
 
 
 }
 
 

diff --git a/src/math/math_FunctionSetRoot.lxx b/src/math/math_FunctionSetRoot.lxx
index 2cae557..c76f56c 100755 (executable)
--- a/src/math/math_FunctionSetRoot.lxx
+++ b/src/math/math_FunctionSetRoot.lxx
@@ -66,3 +66,7 @@ inline Standard_Integer math_FunctionSetRoot::NbIterations() const{
   return Kount;
 }    
 
   return Kount;
 }    
 

+inline Standard_Boolean math_FunctionSetRoot::IsDivergent() const
+{
+  return myIsDivergent;
+}

diff --git a/tests/bugs/moddata_3/bug23830 b/tests/bugs/moddata_3/bug23830
new file mode 100755 (executable)
index 0000000..4f67c88
--- /dev/null
+++ b/tests/bugs/moddata_3/bug23830
@@ -0,0 +1,41 @@
+puts "================"
+puts "OCC23830"
+puts "================"
+puts ""
+#######################################################################
+# BRepExtrema_DistShapeShape does not find intersection of face with edge
+#######################################################################
+
+restore [locate_data_file bug23830_face.brep] s
+mksurface s s
+
+set Indices {1 2 6 7 11 16 21 36 41 42}
+
+foreach i ${Indices} {
+  restore [locate_data_file bug23830_circle$i.brep] c_$i
+  mkcurve c_$i c_$i
+  extrema c_$i s
+  if { [isdraw ext_1] } {
+    mkedge e ext_1
+    regexp {Mass +: +([-0-9.+eE]+)} [lprops e] full l
+    if {$l > 1e-12} {
+      puts "Error: invalid result"
+    }
+    renamevar ext_1 r_$i
+  }
+}
+
+smallview
+l
+l
+l
+l
+l
+l
+l
+l
+erase
+display s
+foreach i ${Indices} { if {[isdraw r_$i]} {display r_$i} }
+fit
+set only_screen_axo 1