0027310: Huge tolerance obtained in the result of intersection of two cylindrical...

Sometimes start point of the intersection line is in the surface boundary strictly. I.e. the parameter of this point in the surface can be equal to both 0 or 2*PI equivalently. It is important to chose correct parameter value.

The algorithm of prediction is based on monotonicity property (see CylCylMonotonicity(...) function in IntPatch_ImpImpIntersection_4.gxx). Now, this function is used wrongly. The fix improves this situation.

Small optimization in the code.
Creation of test cases .

The logic of returning value by the method BoundariesComputing() has been corrected.

diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
index 8673a17..cee3358 100644 (file)
--- a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
+++ b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
@@ -1607,6 +1607,8 @@ static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
 //=======================================================================
 static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
                                             const Standard_Integer theNOfMembers,
+                                            const Standard_Real theUSurf1f,
+                                            const Standard_Real theUSurf1l,
                                             const Standard_Real theTol)
 {
   Standard_Boolean aRetVal = Standard_False;
@@ -1622,7 +1624,10 @@ static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
 
     if((anA-anB) < theTol)
     {
-      anA = (anA + anB)/2.0;
+      if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l)) 
+        anA = (anA + anB)/2.0;
+      else
+        anA = anB;
 
       //Make this element infinite an forget it
       //(we will not use it in next iterations).
@@ -2092,9 +2097,170 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
   }
 }
 
+//=======================================================================
+//function : BoundariesComputing
+//purpose  : Computes true domain of future intersection curve (allows
+//            avoiding excess iterations)
+//=======================================================================
+//=======================================================================
+static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
+                                            const Standard_Real thePeriod,
+                                            Standard_Real theU1f[],
+                                            Standard_Real theU1l[])
+{
+  if(theCoeffs.mB > 0.0)
+  {
+    if(theCoeffs.mB + Abs(theCoeffs.mC) < -1.0)
+    {//There is NOT intersection
+      return Standard_False;
+    }
+    else if(theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+    {//U=[0;2*PI]+aFI1
+      theU1f[0] = theCoeffs.mFI1;
+      theU1l[0] = thePeriod + theCoeffs.mFI1;
+    }
+    else if((1 + theCoeffs.mC <= theCoeffs.mB) &&
+            (theCoeffs.mB <= 1 - theCoeffs.mC))
+    {
+      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
+      if(anArg > 1.0)
+        anArg = 1.0;
+      if(anArg < -1.0)
+        anArg = -1.0;
+
+      const Standard_Real aDAngle = acos(anArg);
+      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
+      theU1f[0] = theCoeffs.mFI1;
+      theU1l[0] = aDAngle + theCoeffs.mFI1;
+      theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1;
+      theU1l[1] = thePeriod + theCoeffs.mFI1;
+    }
+    else if((1 - theCoeffs.mC <= theCoeffs.mB) &&
+            (theCoeffs.mB <= 1 + theCoeffs.mC))
+    {
+      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
+      if(anArg > 1.0)
+        anArg = 1.0;
+      if(anArg < -1.0)
+        anArg = -1.0;
+
+      const Standard_Real aDAngle = acos(anArg);
+      //U=[aDAngle;2*PI-aDAngle]+aFI1
+
+      theU1f[0] = aDAngle + theCoeffs.mFI1;
+      theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1;
+    }
+    else if(theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
+    {
+      Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB,
+                    anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB;
+      if(anArg1 > 1.0)
+        anArg1 = 1.0;
+      if(anArg1 < -1.0)
+        anArg1 = -1.0;
+
+      if(anArg2 > 1.0)
+        anArg2 = 1.0;
+      if(anArg2 < -1.0)
+        anArg2 = -1.0;
+
+      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
+      //(U=[aDAngle1;aDAngle2]+aFI1) ||
+      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
+
+      theU1f[0] = aDAngle1 + theCoeffs.mFI1;
+      theU1l[0] = aDAngle2 + theCoeffs.mFI1;
+      theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1;
+      theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1;
+    }
+    else
+    {
+      return Standard_False;
+    }
+  }
+  else if(theCoeffs.mB < 0.0)
+  {
+    if(theCoeffs.mB + Abs(theCoeffs.mC) > 1.0)
+    {//There is NOT intersection
+      return Standard_False;
+    }
+    else if(-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+    {//U=[0;2*PI]+aFI1
+      theU1f[0] = theCoeffs.mFI1;
+      theU1l[0] = thePeriod + theCoeffs.mFI1;
+    }
+    else if((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
+            ( theCoeffs.mB <= theCoeffs.mC - 1))
+    {
+      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
+      if(anArg > 1.0)
+        anArg = 1.0;
+      if(anArg < -1.0)
+        anArg = -1.0;
+
+      const Standard_Real aDAngle = acos(anArg);
+      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
+
+      theU1f[0] = theCoeffs.mFI1;
+      theU1l[0] = aDAngle + theCoeffs.mFI1;
+      theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1;
+      theU1l[1] = thePeriod + theCoeffs.mFI1;
+    }
+    else if((theCoeffs.mC - 1 <= theCoeffs.mB) &&
+            (theCoeffs.mB <= -theCoeffs.mB - 1))
+    {
+      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
+      if(anArg > 1.0)
+        anArg = 1.0;
+      if(anArg < -1.0)
+        anArg = -1.0;
+
+      const Standard_Real aDAngle = acos(anArg);
+      //U=[aDAngle;2*PI-aDAngle]+aFI1
+
+      theU1f[0] = aDAngle + theCoeffs.mFI1;
+      theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1;
+    }
+    else if(-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
+    {
+      Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB,
+                    anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB;
+      if(anArg1 > 1.0)
+        anArg1 = 1.0;
+      if(anArg1 < -1.0)
+        anArg1 = -1.0;
+
+      if(anArg2 > 1.0)
+        anArg2 = 1.0;
+      if(anArg2 < -1.0)
+        anArg2 = -1.0;
+
+      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
+      //(U=[aDAngle1;aDAngle2]+aFI1) ||
+      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
+
+      theU1f[0] = aDAngle1 + theCoeffs.mFI1;
+      theU1l[0] = aDAngle2 + theCoeffs.mFI1;
+      theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1;
+      theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1;
+    }
+    else
+    {
+      return Standard_False;
+    }
+  }
+  else
+  {
+    return Standard_False;
+  }
+
+  return Standard_True;
+}
+
 //=======================================================================
 //function : CriticalPointsComputing
-//purpose  : theNbCritPointsMax contains true number of critical points
+//purpose  : theNbCritPointsMax contains true number of critical points.
+//            It must be initialized correctly before function calling
 //=======================================================================
 static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
                                     const Standard_Real theUSurf1f,
@@ -2106,15 +2272,15 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
                                     Standard_Integer& theNbCritPointsMax,
                                     Standard_Real theU1crit[])
 {
-  //[0...1] - in these points parameter U1 goes through
-  //          the seam-edge of the first cylinder.
-  //[2...3] - First and last U1 parameter.
-  //[4...5] - in these points parameter U2 goes through
-  //          the seam-edge of the second cylinder.
-  //[6...9] - in these points an intersection line goes through
-  //          U-boundaries of the second surface.
-
-  theNbCritPointsMax = 10;
+  //[0...1]   - in these points parameter U1 goes through
+  //            the seam-edge of the first cylinder.
+  //[2...3]   - First and last U1 parameter.
+  //[4...5]   - in these points parameter U2 goes through
+  //            the seam-edge of the second cylinder.
+  //[6...9]   - in these points an intersection line goes through
+  //            U-boundaries of the second surface.
+  //[10...11] - Boundary of monotonicity interval of U2(U1) function
+  //            (see CylCylMonotonicity() function)
 
   theU1crit[0] = 0.0;
   theU1crit[1] = thePeriod;
@@ -2154,6 +2320,9 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
                   acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
                   Precision::Infinite();
 
+  theU1crit[10] = theCoeffs.mFI1;
+  theU1crit[11] = M_PI+theCoeffs.mFI1;
+
   //preparative treatment of array. This array must have faled to contain negative
   //infinity number
 
@@ -2175,7 +2344,8 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
   {
     std::sort(theU1crit, theU1crit + theNbCritPointsMax);
   }
-  while(ExcludeNearElements(theU1crit, theNbCritPointsMax, theTol2D));
+  while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
+                            theUSurf1f, theUSurf1l, theTol2D));
 
   //Here all not infinite elements in theU1crit are different and sorted.
 
@@ -2281,151 +2451,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
   Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
   Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
 
-  if(anEquationCoeffs.mB > 0.0)
-  {
-    if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0)
-    {//There is NOT intersection
-      return Standard_True;
-    }
-    else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
-    {//U=[0;2*PI]+aFI1
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
-    }
-    else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
-            (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC))
-    {
-      Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
-      if(anArg > 1.0)
-        anArg = 1.0;
-      if(anArg < -1.0)
-        anArg = -1.0;
-
-      const Standard_Real aDAngle = acos(anArg);
-      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
-    }
-    else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
-            (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC))
-    {
-      Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
-      if(anArg > 1.0)
-        anArg = 1.0;
-      if(anArg < -1.0)
-        anArg = -1.0;
-
-      const Standard_Real aDAngle = acos(anArg);
-      //U=[aDAngle;2*PI-aDAngle]+aFI1
-
-      aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
-    }
-    else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
-    {
-      Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB,
-                    anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
-      if(anArg1 > 1.0)
-        anArg1 = 1.0;
-      if(anArg1 < -1.0)
-        anArg1 = -1.0;
-
-      if(anArg2 > 1.0)
-        anArg2 = 1.0;
-      if(anArg2 < -1.0)
-        anArg2 = -1.0;
-
-      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
-      //(U=[aDAngle1;aDAngle2]+aFI1) ||
-      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
-
-      aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
-    }
-    else
-    {
-      Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
-    }
-  }
-  else if(anEquationCoeffs.mB < 0.0)
-  {
-    if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0)
-    {//There is NOT intersection
-      return Standard_True;
-    }
-    else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
-    {//U=[0;2*PI]+aFI1
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
-    }
-    else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
-            ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1))
-    {
-      Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
-      if(anArg > 1.0)
-        anArg = 1.0;
-      if(anArg < -1.0)
-        anArg = -1.0;
-
-      const Standard_Real aDAngle = acos(anArg);
-      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
-
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
-    }
-    else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
-            (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1))
-    {
-      Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
-      if(anArg > 1.0)
-        anArg = 1.0;
-      if(anArg < -1.0)
-        anArg = -1.0;
-
-      const Standard_Real aDAngle = acos(anArg);
-      //U=[aDAngle;2*PI-aDAngle]+aFI1
-
-      aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
-    }
-    else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
-    {
-      Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB,
-                    anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
-      if(anArg1 > 1.0)
-        anArg1 = 1.0;
-      if(anArg1 < -1.0)
-        anArg1 = -1.0;
-
-      if(anArg2 > 1.0)
-        anArg2 = 1.0;
-      if(anArg2 < -1.0)
-        anArg2 = -1.0;
-
-      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
-      //(U=[aDAngle1;aDAngle2]+aFI1) ||
-      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
-
-      aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
-    }
-    else
-    {
-      Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
-    }
-  }
-  else
-  {
-    Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!");
-  }
+  if(!BoundariesComputing(anEquationCoeffs, aPeriod, aU1f, aU1l))
+    return Standard_True;
 
   for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
   {
@@ -2450,7 +2477,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
   }
 
   //Critical points
-  const Standard_Integer aNbCritPointsMax = 10;
+  const Standard_Integer aNbCritPointsMax = 12;
   Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
@@ -2460,6 +2487,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
+                                              Precision::Infinite(),
+                                              Precision::Infinite(),
                                               Precision::Infinite()};
 
   Standard_Integer aNbCritPoints = aNbCritPointsMax;
@@ -2622,7 +2651,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
               //correct value.
 
               Standard_Boolean isIncreasing = Standard_True;
-              CylCylMonotonicity(anU1, i, anEquationCoeffs, aPeriod, isIncreasing);
+              CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs, aPeriod, isIncreasing);
 
               //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
               //then after the next step (when U1 will be increased) U2 will be

diff --git a/tests/bugs/modalg_6/bug27310_1 b/tests/bugs/modalg_6/bug27310_1
new file mode 100644 (file)
index 0000000..3aff9cf
--- /dev/null
+++ b/tests/bugs/modalg_6/bug27310_1
@@ -0,0 +1,33 @@
+puts "========"
+puts "OCC27310"
+puts "========"
+puts ""
+#################################################
+# Huge tolerance obtained in the result of intersection of two cylindrical faces
+#################################################
+
+set ExpTol 1.0e-7
+set GoodNbCurv 2
+
+restore [locate_data_file OCC496a.brep] a 
+restore [locate_data_file OCC496b.brep] b
+
+explode a f
+explode b f
+
+set log [bopcurves a_8 b_2 -2d]
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} ${log} full Toler NbCurv
+
+if {${NbCurv} != ${GoodNbCurv}} {
+  puts "Error: Number of curves is bad!"
+}
+
+checkreal TolReached $Toler $ExpTol 0.0 0.1
+
+smallview
+don c_*
+fit
+disp a_8 b_2
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

diff --git a/tests/bugs/modalg_6/bug27310_2 b/tests/bugs/modalg_6/bug27310_2
new file mode 100644 (file)
index 0000000..30e5754
--- /dev/null
+++ b/tests/bugs/modalg_6/bug27310_2
@@ -0,0 +1,33 @@
+puts "========"
+puts "OCC27310"
+puts "========"
+puts ""
+#################################################
+# Huge tolerance obtained in the result of intersection of two cylindrical faces
+#################################################
+
+set ExpTol 7.7015195151142059e-006
+set GoodNbCurv 2
+
+restore [locate_data_file OCC496a.brep] a 
+restore [locate_data_file OCC496b.brep] b
+
+explode a f
+explode b f
+
+set log [bopcurves a_10 b_4 -2d]
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} ${log} full Toler NbCurv
+
+if {${NbCurv} != ${GoodNbCurv}} {
+  puts "Error: Number of curves is bad!"
+}
+
+checkreal TolReached $Toler $ExpTol 0.0 0.1
+
+smallview
+don c_*
+fit
+disp a_10 b_4
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png