From b007889efd4ca591c55faba556a48eae90d848ca Mon Sep 17 00:00:00 2001
From: ifv <ifv@opencascade.com>
Date: Tue, 12 Jan 2021 16:14:35 +0300
Subject: [PATCH] 0032058: Modeling Data - Extrema curve-surface gives wrong
 result for planar surface of revolunion and circle

Extrema_GenExtCS.cxx : added solution refinement, if solution seems to be "bad" according to special criteria.
Extrema_FuncExtCS.cxx : "cosmetic" modifications
BOPAlgo_PaveFiller_3.cxx : adding control of shape index to prevent exception in ShapeInfo
bugs/moddata_3/bug32058 : new test
---
 src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx |   6 ++
 src/Extrema/Extrema_FuncExtCS.cxx    |   5 +-
 src/Extrema/Extrema_GenExtCS.cxx     | 100 ++++++++++++++++++++++++++-
 tests/bugs/moddata_3/bug32058        |  15 ++++
 4 files changed, 121 insertions(+), 5 deletions(-)
 create mode 100644 tests/bugs/moddata_3/bug32058

diff --git a/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx b/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx
index c199a87c8c..0c0835ecbe 100644
--- a/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx
+++ b/src/BOPAlgo/BOPAlgo_PaveFiller_3.cxx
@@ -655,6 +655,12 @@ void BOPAlgo_PaveFiller::FillShrunkData(Handle(BOPDS_PaveBlock)& thePB)
   // Vertices
   Standard_Integer nV1, nV2;
   thePB->Indices(nV1, nV2);
+
+  if (nV1 < 0 || nV2 < 0)
+  {
+    return;
+  }
+
   const TopoDS_Vertex& aV1=(*(TopoDS_Vertex *)(&myDS->Shape(nV1))); 
   const TopoDS_Vertex& aV2=(*(TopoDS_Vertex *)(&myDS->Shape(nV2))); 
   // Get the edge
diff --git a/src/Extrema/Extrema_FuncExtCS.cxx b/src/Extrema/Extrema_FuncExtCS.cxx
index ffacecc19c..230dfeb17a 100644
--- a/src/Extrema/Extrema_FuncExtCS.cxx
+++ b/src/Extrema/Extrema_FuncExtCS.cxx
@@ -208,12 +208,13 @@ Standard_Integer Extrema_FuncExtCS::GetStateNumber()
   std::cout <<"F(1)= "<<Sol(1)<<" F(2)= "<<Sol(2)<<" F(3)= "<<Sol(3)<<std::endl;
 #endif
   //comparison of solution with previous solutions
-  Standard_Real tol2d = Precision::PConfusion() * Precision::PConfusion();
+  Standard_Real tol2d = Precision::SquarePConfusion();
   Standard_Integer i = 1, nbSol = mySqDist.Length();
   for( ; i <=  nbSol; i++)
   {
     Standard_Real aT = myPoint1(i).Parameter();
-    if( (myt - aT) * (myt - aT) <= tol2d )
+    aT -= myt; aT *= aT;
+    if( aT <= tol2d )
       break;
   }
   if (i <= nbSol)
diff --git a/src/Extrema/Extrema_GenExtCS.cxx b/src/Extrema/Extrema_GenExtCS.cxx
index c9aa59c18a..a193ba32b7 100644
--- a/src/Extrema/Extrema_GenExtCS.cxx
+++ b/src/Extrema/Extrema_GenExtCS.cxx
@@ -17,8 +17,6 @@
 
 #include <Extrema_GenExtCS.hxx>
 #include <Adaptor3d_Curve.hxx>
-#include <Adaptor3d_Curve.hxx>
-#include <Adaptor3d_Surface.hxx>
 #include <Adaptor3d_Surface.hxx>
 #include <Geom_OffsetCurve.hxx>
 #include <Extrema_GlobOptFuncCS.hxx>
@@ -38,6 +36,8 @@
 #include <TColgp_Array1OfPnt.hxx>
 #include <Geom_TrimmedCurve.hxx>
 #include <ElCLib.hxx>
+#include <Extrema_GenLocateExtPS.hxx>
+
 
 const Standard_Real MaxParamVal = 1.0e+10;
 const Standard_Real aBorderDivisor = 1.0e+4;
@@ -331,8 +331,8 @@ void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C,
   // Find min approximation
   math_FunctionSetRoot anA(myF, Tol);
   anA.Perform(myF, TUV, TUVinf, TUVsup);
-
   myDone = Standard_True;
+
 }
 //=======================================================================
 //function : GlobMinGenCS
@@ -447,6 +447,7 @@ void Extrema_GenExtCS::GlobMinConicS(const Adaptor3d_Curve& theC,
     anUVsup(i) = theTUVsup(i + 1);
   }
   //
+  //
   math_PSOParticlesPool aParticles(theNbParticles, aNbVar);
 
   math_Vector aMinUV(1, aNbVar);
@@ -517,10 +518,103 @@ void Extrema_GenExtCS::GlobMinConicS(const Adaptor3d_Curve& theC,
       aCT = ElCLib::InPeriod(aCT, theTUVinf(1), theTUVinf(1) + 2. * M_PI);
     }
   }
+
   theTUV(1) = aCT;
   theTUV(2) = anUV(1);
   theTUV(3) = anUV(2);
 
+  Standard_Boolean isBadSol = Standard_False;
+  gp_Vec aDU, aDV, aDT;
+  gp_Pnt aPOnS, aPOnC;
+  myS->D1(anUV(1), anUV(2), aPOnS, aDU, aDV);
+  theC.D1(aCT, aPOnC, aDT);
+  Standard_Real aSqDist = aPOnC.SquareDistance(aPOnS);
+  if (aSqDist <= Precision::SquareConfusion())
+    return;
+
+  gp_Vec aN = aDU.Crossed(aDV);
+  if (aN.SquareMagnitude() < Precision::SquareConfusion())
+    return;
+
+  gp_Vec PcPs(aPOnC, aPOnS);
+
+  Standard_Real anAngMin = M_PI_2 - M_PI_2 / 10.;
+  Standard_Real anAngMax = M_PI_2 + M_PI_2 / 10.;
+
+  Standard_Real anAngN = PcPs.Angle(aN);
+  if (anAngN >= anAngMin && anAngN <= anAngMax)
+  {
+    // PcPs is perpendicular to surface normal, it means that
+    // aPOnC can be on surface, but far from aPOnS
+    isBadSol = Standard_True;
+    Standard_Integer iu, iv;
+    for (iu = -1; iu <= 1; ++iu)
+    {
+      Standard_Real u = anUV(1) + iu * aStepSU;
+      u = Max(anUVinf(1), u);
+      u = Min(anUVsup(1), u);
+      for (iv = -1; iv <= 1; ++iv)
+      {
+        Standard_Real v = anUV(2) + iv * aStepSV;
+        v = Max(anUVinf(2), v);
+        v = Min(anUVsup(2), v);
+        myS->D1(u, v, aPOnS, aDU, aDV);
+        if (aPOnC.SquareDistance(aPOnS) < Precision::SquareConfusion())
+        {
+          isBadSol = Standard_False;
+          break;
+        }
+        aN = aDU.Crossed(aDV);
+        if (aN.SquareMagnitude() < Precision::SquareConfusion())
+        {
+          isBadSol = Standard_False;
+          break;
+        }
+        PcPs.SetXYZ(aPOnS.XYZ() - aPOnC.XYZ());
+        anAngN = PcPs.Angle(aN);
+        if (anAngN < anAngMin || anAngN > anAngMax)
+        {
+          isBadSol = Standard_False;
+          break;
+        }
+      }
+      if (!isBadSol)
+      {
+        break;
+      }
+    }
+  }
+
+  if (isBadSol)
+  {
+    //Try to precise solution with help of Extrema PS
+
+    math_Vector aF(1, 3);
+    aF(1) = PcPs.Dot(aDT);
+    aF(2) = PcPs.Dot(aDU);
+    aF(3) = PcPs.Dot(aDV);
+    Standard_Real aFF = aF.Norm2();
+
+    Extrema_GenLocateExtPS anExtPS(*myS, mytol2, mytol2);
+    anExtPS.Perform(aPOnC, anUV(1), anUV(2), Standard_False);
+    if (anExtPS.IsDone())
+    {
+      const Extrema_POnSurf& aPmin = anExtPS.Point();
+      aPmin.Parameter(anUV(1), anUV(2));
+      math_Vector aTUV = theTUV;
+      aTUV(2) = anUV(1);
+      aTUV(3) = anUV(2);
+      myF.Value(aTUV, aF);
+      Standard_Real aFF1 = aF.Norm2();
+
+      if (anExtPS.SquareDistance() < aSqDist && aFF1 <= 1.1 * aFF)
+      {
+        theTUV(2) = aTUV(2);
+        theTUV(3) = aTUV(3);
+      }
+    }
+  }
+ 
 }
 //=======================================================================
 //function : GlobMinCQuadric
diff --git a/tests/bugs/moddata_3/bug32058 b/tests/bugs/moddata_3/bug32058
new file mode 100644
index 0000000000..96e53beb72
--- /dev/null
+++ b/tests/bugs/moddata_3/bug32058
@@ -0,0 +1,15 @@
+puts "========="
+puts "0032058: Modeling Data - Extrema curve-surface gives wrong result for planar surface of revolunion and circle"
+puts "========="
+puts ""
+
+
+restore [locate_data_file bug32058_c] c
+restore [locate_data_file bug32058_s] s
+
+extrema c s
+checklength ext_1 -l 2.3437142008433856e-13
+
+smallview
+fit
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png