0025757: distmini returns wrong solution for ellipse/vertex

Analytical handling of degenerated cases added.

Test case for issue CR25757

Correction of test case for issue CR25757

diff --git a/src/math/math_TrigonometricFunctionRoots.cxx b/src/math/math_TrigonometricFunctionRoots.cxx
index 6f61628a66c4e7af98734387e04b058982c79d3a..4d5f4f99a40e16c3d62248e4ac5002c13e8ccc81 100644 (file)
--- a/src/math/math_TrigonometricFunctionRoots.cxx
+++ b/src/math/math_TrigonometricFunctionRoots.cxx
@@ -29,6 +29,7 @@
 #include <Standard_OutOfRange.hxx>
 #include <math_FunctionWithDerivative.hxx>
 #include <math_NewtonFunctionRoot.hxx>
+#include <Precision.hxx>
 
 
 class MyTrigoFunction: public math_FunctionWithDerivative
@@ -178,97 +179,207 @@ void math_TrigonometricFunctionRoots::Perform(const Standard_Real A,
   if ((Abs(A) <= Eps) && (Abs(B) <= Eps)) {
     if (Abs(C) <= Eps) {
       if (Abs(D) <= Eps) {
-       if (Abs(E) <= Eps) {
-         InfiniteStatus = Standard_True;   // infinite de solutions.
-         return;
-       }
-       else {
-         NbSol = 0;
-         return;
-       }
+        if (Abs(E) <= Eps) {
+          InfiniteStatus = Standard_True;   // infinite de solutions.
+          return;
+        }
+        else {
+          NbSol = 0;
+          return;
+        }
       }
       else { 
-// Equation du type d*sin(x) + e = 0
-// =================================   
-       NbSol = 0;
-       AA = -E/D;
-       if (Abs(AA) > 1.) {
-         return;
-       }
-
-       Zer(1) = ASin(AA);
-       Zer(2) = M_PI - Zer(1);
-       NZer = 2;
-       for (i = 1; i <= NZer; i++) {
-         if (Zer(i) <= -Eps) {
-           Zer(i) = Depi - Abs(Zer(i));
-         }
-         // On rend les solutions entre InfBound et SupBound:
-         // =================================================
-         Zer(i) += IntegerPart(Mod)*Depi;
-         X = Zer(i)-MyBorneInf;
-         if ((X > (-Epsilon(Delta))) && (X < Delta+ Epsilon(Delta))) {
-           NbSol++;
-           Sol(NbSol) = Zer(i);
-         }
-       }
+        // Equation du type d*sin(x) + e = 0
+        // =================================   
+        NbSol = 0;
+        AA = -E/D;
+        if (Abs(AA) > 1.) {
+          return;
+        }
+
+        Zer(1) = ASin(AA);
+        Zer(2) = M_PI - Zer(1);
+        NZer = 2;
+        for (i = 1; i <= NZer; i++) {
+          if (Zer(i) <= -Eps) {
+            Zer(i) = Depi - Abs(Zer(i));
+          }
+          // On rend les solutions entre InfBound et SupBound:
+          // =================================================
+          Zer(i) += IntegerPart(Mod)*Depi;
+          X = Zer(i)-MyBorneInf;
+          if ((X > (-Epsilon(Delta))) && (X < Delta+ Epsilon(Delta))) {
+            NbSol++;
+            Sol(NbSol) = Zer(i);
+          }
+        }
       }
       return;
     }
     else if (Abs(D) <= Eps)  {
 
-// Equation du premier degre de la forme c*cos(x) + e = 0
-// ======================================================        
+      // Equation du premier degre de la forme c*cos(x) + e = 0
+      // ======================================================          
       NbSol = 0;
       AA = -E/C;
       if (Abs(AA) >1.) {
-       return;
+        return;
       }
       Zer(1) = ACos(AA);
       Zer(2) = -Zer(1);
       NZer = 2;
 
       for (i = 1; i <= NZer; i++) {
-       if (Zer(i) <= -Eps) {
-         Zer(i) = Depi-Abs(Zer(i));
-       }
-       // On rend les solutions entre InfBound et SupBound:
-       // =================================================
-       Zer(i) += IntegerPart(Mod)*2.*M_PI;
-       X = Zer(i)-MyBorneInf;
-       if ((X >= (-Epsilon(Delta))) && (X <= Delta+ Epsilon(Delta))) {
-         NbSol++;
-         Sol(NbSol) = Zer(i);
-       }
+        if (Zer(i) <= -Eps) {
+          Zer(i) = Depi - Abs(Zer(i));
+        }
+        // On rend les solutions entre InfBound et SupBound:
+        // =================================================
+        Zer(i) += IntegerPart(Mod)*2.*M_PI;
+        X = Zer(i)-MyBorneInf;
+        if ((X >= (-Epsilon(Delta))) && (X <= Delta+ Epsilon(Delta))) {
+          NbSol++;
+          Sol(NbSol) = Zer(i);
+        }
       }
       return;
     }
     else {
 
-// Equation du second degre:
-// =========================
+      // Equation du second degre:
+      // =========================
       AA = E - C;
       BB = 2.0*D;
       CC = E + C;
 
       math_DirectPolynomialRoots Resol(AA, BB, CC);
       if (!Resol.IsDone()) {
-       Done = Standard_False;
-       return;
+        Done = Standard_False;
+        return;
       }
       else if(!Resol.InfiniteRoots()) {
-       NZer = Resol.NbSolutions();
-       for (i = 1; i <= NZer; i++) {
-         Zer(i) = Resol.Value(i);
-       }
+        NZer = Resol.NbSolutions();
+        for (i = 1; i <= NZer; i++) {
+          Zer(i) = Resol.Value(i);
+        }
       }
       else if (Resol.InfiniteRoots()) {
-       InfiniteStatus = Standard_True;
-       return;
+        InfiniteStatus = Standard_True;
+        return;
       }
     }
   }
   else {
+    // Two additional analytical cases.
+    if ((Abs(A) <= Eps) && 
+        (Abs(E) <= Eps))
+    {
+      if (Abs(C) <= Eps)
+      {
+        // 2 * B * sin * cos + D * sin = 0
+        NZer = 2;
+        Zer(1) = 0.0;
+        Zer(2) = M_PI;
+
+        AA = -D/(B * 2);
+        if (Abs(AA) <= 1.0 + Precision::PConfusion())
+        {
+          NZer = 4;
+          if (AA >= 1.0)
+          {
+            Zer(3)= 0.0;
+            Zer(4)= 0.0;
+          }
+          else if (AA <= -1.0)
+          {
+            Zer(3)= M_PI;
+            Zer(4)= M_PI;
+          }
+          else
+          {
+            Zer(3)= ACos(AA);
+            Zer(4) = Depi - Zer(3);
+          }
+        }
+
+        NbSol = 0;
+        for (i = 1; i <= NZer; i++) 
+        {
+          if (Zer(i) <= MyBorneInf - Eps)
+          {
+            Zer(i) += Depi;
+          }
+          // On rend les solutions entre InfBound et SupBound:
+          // =================================================
+          Zer(i) += IntegerPart(Mod)*2.*M_PI;
+          X = Zer(i)-MyBorneInf;
+          if ((X >= (-Precision::PConfusion())) && 
+              (X <= Delta + Precision::PConfusion()))
+          {
+            if (Zer(i) < InfBound)
+              Zer(i) = InfBound;
+            if (Zer(i) > SupBound)
+              Zer(i) = SupBound;
+            NbSol++;
+            Sol(NbSol) = Zer(i);
+          }
+        }
+        return;
+      }
+      if (Abs(D) <= Eps)
+      {
+        // 2 * B * sin * cos + C * cos = 0
+        NZer = 2;
+        Zer(1) = M_PI / 2.0;
+        Zer(2) = M_PI * 3.0 / 2.0;
+
+        AA = -C/(B * 2);
+        if (Abs(AA) <= 1.0 + Precision::PConfusion())
+        {
+          NZer = 4;
+          if (AA >= 1.0)
+          {
+            Zer(3) = M_PI / 2.0;
+            Zer(4) = M_PI / 2.0;
+          }
+          else if (AA <= -1.0)
+          {
+            
+            Zer(3) = M_PI * 3.0 / 2.0;
+            Zer(4) = M_PI * 3.0 / 2.0;
+          }
+          else
+          {
+            Zer(3)= ASin(AA);
+            Zer(4) = M_PI - Zer(3);
+          }
+        }
+
+        NbSol = 0;
+        for (i = 1; i <= NZer; i++) 
+        {
+          if (Zer(i) <= MyBorneInf - Eps)
+          {
+            Zer(i) += Depi;
+          }
+          // On rend les solutions entre InfBound et SupBound:
+          // =================================================
+          Zer(i) += IntegerPart(Mod)*2.*M_PI;
+          X = Zer(i)-MyBorneInf;
+          if ((X >= (-Precision::PConfusion())) && 
+              (X <= Delta + Precision::PConfusion()))
+          {
+            if (Zer(i) < InfBound)
+              Zer(i) = InfBound;
+            if (Zer(i) > SupBound)
+              Zer(i) = SupBound;
+            NbSol++;
+            Sol(NbSol) = Zer(i);
+          }
+        }
+        return;
+      }
+    }
 
 // Equation du 4 ieme degre
 // ========================

diff --git a/tests/bugs/fclasses/bug25757 b/tests/bugs/fclasses/bug25757
new file mode 100755 (executable)
index 0000000..f16c875
--- /dev/null
+++ b/tests/bugs/fclasses/bug25757
@@ -0,0 +1,29 @@
+puts "========"
+puts "OCC25757"
+puts "========"
+puts ""
+##############################################
+# distmini returns wrong solution for ellipse/vertex
+##############################################
+
+restore [locate_data_file bug25757_ellipse.brep] ellipse
+
+vertex vertex1 7.32050807568877 3.999999999999999 10.0
+vertex vertex2 7.32050807568877 3.99999999999999 10.0
+
+distmini dv1 vertex1 ellipse
+set dist1 [dval dv1_val]
+puts "vertex1 distance = ${dist1}"
+
+distmini dv2 vertex2 ellipse
+set dist2 [dval dv2_val]
+puts "vertex2 distance = ${dist2}"
+
+set tol_abs_dist 1.0e-12
+set tol_rel_dist 1.0e-2
+
+set expected_dist1 0.0
+set expected_dist2 0.0
+
+checkreal "Distance 1" ${dist1} ${expected_dist1} ${tol_abs_dist} ${tol_rel_dist}
+checkreal "Distance 2" ${dist2} ${expected_dist2} ${tol_abs_dist} ${tol_rel_dist}