0029535: Modeling Algorithms - BRepExtrema_DistShapeShape returns only one solution but the wires have two intersections (for 6.9.1)

src/Extrema/Extrema_GenExtCC.gxx

@@ -25,6 +25,9 @@
 #include <Precision.hxx>
 #include <NCollection_Vector.hxx>
 #include <NCollection_CellFilter.hxx>
+#include <TColStd_HArray1OfReal.hxx>
+#include <TColStd_Array1OfInteger.hxx>
+#include <GCPnts_AbscissaPoint.hxx>
 
 // Comparator, used in std::sort.
 static Standard_Boolean comp(const gp_XY& theA,
@@ -45,6 +48,57 @@ static Standard_Boolean comp(const gp_XY& theA,
 
   return Standard_False;
 }
+static void ChangeIntervals(Handle(TColStd_HArray1OfReal)& theInts, const Standard_Integer theNbInts)
+{
+  Standard_Integer aNbInts = theInts->Length() - 1;
+  Standard_Integer aNbAdd = theNbInts - aNbInts;
+  Handle(TColStd_HArray1OfReal) aNewInts = new TColStd_HArray1OfReal(1, theNbInts + 1);
+  Standard_Integer aNbLast = theInts->Length();
+  Standard_Integer i;
+  if (aNbInts == 1)
+  {
+    aNewInts->SetValue(1, theInts->Value(1));
+    aNewInts->SetValue(theNbInts + 1, theInts->Value(aNbLast));
+    Standard_Real dt = (theInts->Value(aNbLast) - theInts->Value(1)) / theNbInts;
+    Standard_Real t = theInts->Value(1) + dt;
+    for (i = 2; i <= theNbInts; ++i, t += dt)
+    {
+      aNewInts->SetValue(i, t);
+    }
+    theInts = aNewInts;
+    return;
+  }
+  //
+  for (i = 1; i <= aNbLast; ++i)
+  {
+    aNewInts->SetValue(i, theInts->Value(i));
+  }
+  //
+  while (aNbAdd > 0)
+  {
+    Standard_Real anLIntMax = -1.;
+    Standard_Integer aMaxInd = -1;
+    for (i = 1; i < aNbLast; ++i)
+    {
+      Standard_Real anL = aNewInts->Value(i + 1) - aNewInts->Value(i);
+      if (anL > anLIntMax)
+      {
+        anLIntMax = anL;
+        aMaxInd = i;
+      }
+    }
+
+    Standard_Real t = (aNewInts->Value(aMaxInd + 1) + aNewInts->Value(aMaxInd)) / 2.;
+    for (i = aNbLast; i > aMaxInd; --i)
+    {
+      aNewInts->SetValue(i + 1, aNewInts->Value(i));
+    }
+    aNbLast++;
+    aNbAdd--;
+    aNewInts->SetValue(aMaxInd + 1, t);
+  }
+  theInts = aNewInts;
+}
 
 class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY
 {
@@ -190,15 +244,79 @@ void Extrema_GenExtCC::Perform()
   Curve1 &C1 = *(Curve1*)myC[0];
   Curve2 &C2 = *(Curve2*)myC[1];
 
+  GeomAbs_Shape aContinuity = GeomAbs_C2;
   Standard_Integer aNbInter[2];
-  aNbInter[0] = C1.NbIntervals(GeomAbs_C2);
+  aNbInter[0] = C1.NbIntervals(aContinuity);
-  aNbInter[1] = C2.NbIntervals(GeomAbs_C2);
+  aNbInter[1] = C2.NbIntervals(aContinuity);
-  TColStd_Array1OfReal anIntervals1(1, aNbInter[0] + 1);
+  if (aNbInter[0] * aNbInter[1] > 100)
-  TColStd_Array1OfReal anIntervals2(1, aNbInter[1] + 1);
+  {
-  C1.Intervals(anIntervals1, GeomAbs_C2);
+    aContinuity = GeomAbs_C1;
-  C2.Intervals(anIntervals2, GeomAbs_C2);
+    aNbInter[0] = C1.NbIntervals(aContinuity);
+    aNbInter[1] = C2.NbIntervals(aContinuity);
+  }
+
+  Standard_Real anL[2];
+  Standard_Integer indmax = -1, indmin = -1;
+  const Standard_Real mult = 20.;
+  if (!(Precision::IsInfinite(C1.FirstParameter()) || Precision::IsInfinite(C1.LastParameter()) ||
+        Precision::IsInfinite(C2.FirstParameter()) || Precision::IsInfinite(C2.LastParameter())))
+  {
+    anL[0] = GCPnts_AbscissaPoint::Length(C1);
+    anL[1] = GCPnts_AbscissaPoint::Length(C2);
+    if (anL[0] / aNbInter[0] > mult * anL[1] / aNbInter[1])
+    {
+      indmax = 0;
+      indmin = 1;
+    }
+    else if (anL[1] / aNbInter[1] > mult * anL[0] / aNbInter[0])
+    {
+      indmax = 1;
+      indmin = 0;
+    }
+  }
+  Standard_Integer aNbIntOpt = 0;
+  if (indmax >= 0)
+  {
+    aNbIntOpt = RealToInt(anL[indmax] * aNbInter[indmin] / anL[indmin] / (mult / 4.)) + 1;
+    if (aNbIntOpt > 100 || aNbIntOpt < aNbInter[indmax])
+    {
+      indmax = -1;
+    }
+    else
+    {
+      if (aNbIntOpt * aNbInter[indmin] > 100)
+      {
+        aNbIntOpt = 100 / aNbInter[indmin];
+        if (aNbIntOpt < aNbInter[indmax])
+        {
+          indmax = -1;
+        }
+      }
+    }
+  }
+
+  Handle(TColStd_HArray1OfReal) anIntervals1 = new TColStd_HArray1OfReal(1, aNbInter[0] + 1);
+  Handle(TColStd_HArray1OfReal) anIntervals2 = new TColStd_HArray1OfReal(1, aNbInter[1] + 1);
+  C1.Intervals(anIntervals1->ChangeArray1(), aContinuity);
+  C2.Intervals(anIntervals2->ChangeArray1(), aContinuity);
+  if (indmax >= 0)
+  {
+    if (indmax == 0)
+    {
+      //Change anIntervals1
+      ChangeIntervals(anIntervals1, aNbIntOpt);
+      aNbInter[0] = anIntervals1->Length() - 1;
+    }
+    else
+    {
+      //Change anIntervals2;
+      ChangeIntervals(anIntervals2, aNbIntOpt);
+      aNbInter[1] = anIntervals2->Length() - 1;
+    }
+  }
 
   // Lipchitz constant computation.
+  const Standard_Real aMaxLC = 10000.;
   Standard_Real aLC = 9.0; // Default value.
   const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0);
   const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0);
@@ -207,6 +325,17 @@ void Extrema_GenExtCC::Perform()
     aLC = aMaxDer;
 
   Standard_Boolean isConstLockedFlag = Standard_False;
+  //To prevent LipConst to became too small
+  const Standard_Real aCR = 0.001;
+  if (aMaxDer1 / aMaxDer < aCR || aMaxDer2 / aMaxDer < aCR)
+  {
+    isConstLockedFlag = Standard_True;
+  }
+  if (aMaxDer > aMaxLC)
+  {
+    aLC = aMaxLC;
+    isConstLockedFlag = Standard_True;
+  }
   if (C1.GetType() == GeomAbs_Line)
   {
     Standard_Real aMaxDer = 1.0 / C2.Resolution(1.0);
@@ -236,10 +365,11 @@ void Extrema_GenExtCC::Perform()
   aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0.
 
   // Size computed to have cell index inside of int32 value.
-  const Standard_Real aCellSize = Max(anIntervals1.Upper() - anIntervals1.Lower(),
+  const Standard_Real aCellSize = Max(Max(anIntervals1->Value(anIntervals1->Upper()) - anIntervals1->Value(1),
-                                      anIntervals2.Upper() - anIntervals2.Lower())
+                                      anIntervals2->Value(anIntervals2->Upper()) - anIntervals2->Value(1))
-                                  * Precision::PConfusion() / (2.0 * Sqrt(2.0));
+                                      * Precision::PConfusion() / (2.0 * Sqrt(2.0)),
-  Extrema_CCPointsInspector anInspector(Precision::PConfusion());
+                                      Precision::PConfusion());
+  Extrema_CCPointsInspector anInspector(aCellSize);
   NCollection_CellFilter<Extrema_CCPointsInspector> aFilter(aCellSize);
   NCollection_Vector<gp_XY> aPnts;
 
@@ -252,10 +382,10 @@ void Extrema_GenExtCC::Perform()
   {
     for(j = 1; j <= aNbInter[1]; j++)
     {
-      aFirstBorderInterval(1) = anIntervals1(i);
+      aFirstBorderInterval(1) = anIntervals1->Value(i);
-      aFirstBorderInterval(2) = anIntervals2(j); 
+      aFirstBorderInterval(2) = anIntervals2->Value(j); 
-      aSecondBorderInterval(1) = anIntervals1(i + 1);
+      aSecondBorderInterval(1) = anIntervals1->Value(i + 1);
-      aSecondBorderInterval(2) = anIntervals2(j + 1);
+      aSecondBorderInterval(2) = anIntervals2->Value(j + 1);
 
       aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
       aFinder.Perform();