0029289: Wrong derivatives in math_TrigonometricFunctionRoots.cxx file
authorifv <ifv@opencascade.com>
Thu, 9 Nov 2017 14:20:10 +0000 (17:20 +0300)
committerbugmaster <bugmaster@opencascade.com>
Mon, 27 Nov 2017 08:01:21 +0000 (11:01 +0300)
Class MyTrigoFunction is removed from file math_TrigonometricFunctionRoots.cxx.
New class math_TrigonometricEquationFunction with the same functionality is created to provide possibilities
for individual testing.
Expressions for derivatives are corrected.
New Draw command "intconcon" for intersection 2d conic curves is created.
Test command OCC29289 (file QABugs_20.cxx) is created for individual testing math_TrigonometricEquationFunction.
It is used in tests/bugs/modalg_7/bug29289

dox/user_guides/draw_test_harness/draw_test_harness.md
src/GeomliteTest/GeomliteTest_API2dCommands.cxx
src/QABugs/QABugs_20.cxx
src/math/FILES
src/math/math_TrigonometricEquationFunction.hxx [new file with mode: 0644]
src/math/math_TrigonometricFunctionRoots.cxx
src/math/math_TrigonometricFunctionRoots.hxx
tests/bugs/modalg_7/bug29289 [new file with mode: 0644]

index a39b506..d49c821 100644 (file)
@@ -5529,6 +5529,7 @@ surface_radius c pi 3 c1 c2
 
 * **intersect** computes intersections of surfaces; 
 * **2dintersect** computes intersections of 2d curves.
+* **intconcon** computes intersections of 2d conic curves.
 
 @subsubsection occt_draw_6_7_1  intersect
 
@@ -5547,7 +5548,7 @@ plane p 0 0 40 0 1 5
 intersect e c p 
 ~~~~~
 
-@subsubsection occt_draw_6_7_2  dintersect
+@subsubsection occt_draw_6_7_2  2dintersect
 
 Syntax:      
 ~~~~~
@@ -5564,6 +5565,25 @@ ellipse e2 0 0 0 1 5 2
 2dintersect e1 e2 
 ~~~~~
 
+@subsubsection occt_draw_6_7_3 intconcon
+
+Syntax:      
+~~~~~
+intconcon curve1 curve2 
+~~~~~
+
+Displays the intersection points between two 2d curves. 
+Curves must be only conic sections: 2d lines, circles, ellipses,
+hyperbolas, parabolas. Algorithm from IntAna2d_AnaIntersection is used.
+
+**Example:** 
+~~~~~
+# intersect two 2d ellipses 
+ellipse e1 0 0 5 2 
+ellipse e2 0 0 0 1 5 2 
+intconcon e1 e2 
+~~~~~
+
 @subsection occt_draw_6_8  Approximations
 
 Draw provides command to create curves and surfaces by approximation. 
index 28aa55c..11f4daa 100644 (file)
 #include <Geom2d_Circle.hxx>
 #include <IntAna2d_AnaIntersection.hxx>
 #include <IntAna2d_IntPoint.hxx>
+#include <IntAna2d_Conic.hxx>
 #include <IntRes2d_IntersectionPoint.hxx>
+#include <Geom2dAdaptor_GHCurve.hxx>
+#include <memory>
 
 #include <stdio.h>
 #ifdef _WIN32
@@ -359,24 +362,24 @@ static Standard_Integer intersect(Draw_Interpretor& di, Standard_Integer n, cons
 }
 
 //=======================================================================
-//function : intersect
+//function : intersect_ana
 //purpose  : 
 //=======================================================================
 
 static Standard_Integer intersect_ana(Draw_Interpretor& di, Standard_Integer n, const char** a)
 {
-  if( n < 2) 
+  if (n < 2)
   {
-    cout<< "2dintana circle circle "<<endl;
+    cout << "2dintana circle circle " << endl;
     return 1;
   }
-  
+
   Handle(Geom2d_Curve) C1 = DrawTrSurf::GetCurve2d(a[1]);
-  if ( C1.IsNull() && !C1->IsKind(STANDARD_TYPE(Geom2d_Circle))) 
+  if (C1.IsNull() && !C1->IsKind(STANDARD_TYPE(Geom2d_Circle)))
     return 1;
 
   Handle(Geom2d_Curve) C2 = DrawTrSurf::GetCurve2d(a[2]);
-  if ( C2.IsNull() && !C2->IsKind(STANDARD_TYPE(Geom2d_Circle)))
+  if (C2.IsNull() && !C2->IsKind(STANDARD_TYPE(Geom2d_Circle)))
     return 1;
 
   Handle(Geom2d_Circle) aCir1 = Handle(Geom2d_Circle)::DownCast(C1);
@@ -386,11 +389,121 @@ static Standard_Integer intersect_ana(Draw_Interpretor& di, Standard_Integer n,
 
   Standard_Integer i;
 
+  for (i = 1; i <= Intersector.NbPoints(); i++) {
+    gp_Pnt2d P = Intersector.Point(i).Value();
+    di << "Intersection point " << i << " : " << P.X() << " " << P.Y() << "\n";
+    di << "parameter on the fist: " << Intersector.Point(i).ParamOnFirst();
+    di << " parameter on the second: " << Intersector.Point(i).ParamOnSecond() << "\n";
+    Handle(Draw_Marker2D) mark = new Draw_Marker2D(P, Draw_X, Draw_vert);
+    dout << mark;
+  }
+  dout.Flush();
+
+  return 0;
+}
+
+//=======================================================================
+//function : intconcon
+//purpose  : 
+//=======================================================================
+
+static Standard_Integer intconcon(Draw_Interpretor& di, Standard_Integer n, const char** a)
+{
+  if( n < 2) 
+  {
+    cout<< "intconcon con1 con2 "<<endl;
+    return 1;
+  }
+  
+  Handle(Geom2d_Curve) C1 = DrawTrSurf::GetCurve2d(a[1]);
+  if (C1.IsNull())
+  {
+    cout << a[1] << " is Null " << endl;
+    return 1;
+  }
+
+  Handle(Geom2d_Curve) C2 = DrawTrSurf::GetCurve2d(a[2]);
+  if (C2.IsNull())
+  {
+    cout << a[2] << " is Null " << endl;
+    return 1;
+  }
+
+  Geom2dAdaptor_Curve AC1(C1), AC2(C2);
+  GeomAbs_CurveType T1 = AC1.GetType(), T2 = AC2.GetType();
+#if (defined(_MSC_VER) && (_MSC_VER < 1600))
+  std::auto_ptr<IntAna2d_Conic> pCon;
+#else
+  std::unique_ptr<IntAna2d_Conic> pCon;
+#endif  
+  switch (T2)
+  {
+  case GeomAbs_Line:
+  {
+    pCon.reset(new IntAna2d_Conic(AC2.Line()));
+    break;
+  }
+  case GeomAbs_Circle:
+  {
+    pCon.reset(new IntAna2d_Conic(AC2.Circle()));
+    break;
+  }
+  case GeomAbs_Ellipse:
+  {
+    pCon.reset(new IntAna2d_Conic(AC2.Ellipse()));
+    break;
+  }
+  case GeomAbs_Hyperbola:
+  {
+    pCon.reset(new IntAna2d_Conic(AC2.Hyperbola()));
+    break;
+  }
+  case GeomAbs_Parabola:
+  {
+    pCon.reset(new IntAna2d_Conic(AC2.Parabola()));
+    break;
+  }
+  default:
+    cout << a[2] << " is not conic " << endl;
+    return 1;
+  }
+
+  IntAna2d_AnaIntersection Intersector;
+  switch (T1)
+  {
+  case GeomAbs_Line:
+    Intersector.Perform(AC1.Line(), *pCon);
+    break;
+  case GeomAbs_Circle:
+    Intersector.Perform(AC1.Circle(), *pCon);
+    break;
+  case GeomAbs_Ellipse:
+    Intersector.Perform(AC1.Ellipse(), *pCon);
+    break;
+  case GeomAbs_Hyperbola:
+    Intersector.Perform(AC1.Hyperbola(), *pCon);
+    break;
+  case GeomAbs_Parabola:
+    Intersector.Perform(AC1.Parabola(), *pCon);
+    break;
+  default:
+    cout << a[1] << " is not conic " << endl;
+    return 1;
+  }
+  
+  Standard_Integer i;
   for ( i = 1; i <= Intersector.NbPoints(); i++) {
     gp_Pnt2d P = Intersector.Point(i).Value();
     di<<"Intersection point "<<i<<" : "<<P.X()<<" "<<P.Y()<<"\n";
-       di<<"parameter on the fist: "<<Intersector.Point(i).ParamOnFirst();
-       di<<" parameter on the second: "<<Intersector.Point(i).ParamOnSecond()<<"\n";
+    di << "parameter on the fist: " << Intersector.Point(i).ParamOnFirst();
+    if (!Intersector.Point(i).SecondIsImplicit())
+    {
+      di << " parameter on the second: " << Intersector.Point(i).ParamOnSecond() << "\n";
+    }
+    else
+    {
+      di << "\n";
+    }
     Handle(Draw_Marker2D) mark = new Draw_Marker2D( P, Draw_X, Draw_vert); 
     dout << mark;
   }
@@ -430,6 +543,8 @@ void GeomliteTest::API2dCommands(Draw_Interpretor& theCommands)
   theCommands.Add("2dintersect", "intersect curve curve [Tol]",__FILE__,
                  intersect,g);
 
-  theCommands.Add("2dintanalytical", "intersect curve curve using IntAna",__FILE__,
+  theCommands.Add("2dintanalytical", "intersect circle1 and circle2 using IntAna",__FILE__,
                  intersect_ana,g);
+  theCommands.Add("intconcon", "intersect conic curve1 and conic curve2 using IntAna", __FILE__,
+    intconcon, g);
 }
index 4b7a116..df972f2 100644 (file)
@@ -2530,6 +2530,65 @@ static Standard_Integer OCC28131 (Draw_Interpretor&, Standard_Integer theNbArgs,
 */
   return 0;
 }
+#include <math_NewtonFunctionRoot.hxx>
+#include <math_TrigonometricFunctionRoots.hxx>
+#include <math_TrigonometricEquationFunction.hxx>
+#include <gp_Elips2d.hxx>
+#include <Geom2d_Ellipse.hxx>
+#include <Geom2dAPI_InterCurveCurve.hxx>
+static Standard_Integer OCC29289(Draw_Interpretor&, Standard_Integer , const char** )
+{
+  gp_Elips2d e1(gp_Ax2d(gp_Pnt2d(0., 0.), gp_Dir2d(1., 0)), 2., 1.);
+  Handle(Geom2d_Ellipse) Ge1 = new Geom2d_Ellipse(e1);
+  gp_Elips2d e2(gp_Ax2d(gp_Pnt2d(0.5, 0.5), gp_Dir2d(1., 1.)), 2., 1.);
+  Handle(Geom2d_Ellipse) Ge2 = new Geom2d_Ellipse(e2);
+
+  Standard_Integer err = 0;
+  Geom2dAPI_InterCurveCurve Intersector;
+  Intersector.Init(Ge1, Ge2, 1.e-7);
+  if (Intersector.NbPoints() == 0)
+  {
+    std::cout << "Error: intersector is not done  \n";
+    err = 1;
+  }
+
+
+  Standard_Real A, B, C, D, E;
+  A = 1.875;
+  B = -.75;
+  C = -.5;
+  D = -.25;
+  E = -.25;
+  math_TrigonometricEquationFunction MyF(A, B, C, D, E);
+  Standard_Real X, Tol1, Eps, Teta, TetaNewton;
+  Tol1 = 1.e-15;
+  Eps = 1.5e-12;
+  Standard_Integer Nit[] = { 5, 6, 7, 6 };
+
+  Standard_Real TetaPrev = 0.;
+  Standard_Integer i;
+  for (i = 1; i <= Intersector.NbPoints(); i++) {
+    Teta = Intersector.Intersector().Point(i).ParamOnFirst();
+    X = Teta - 0.1 * (Teta - TetaPrev);
+    TetaPrev = Teta;
+    math_NewtonFunctionRoot Resol(MyF, X, Tol1, Eps, Nit[i-1]);
+    if (Resol.IsDone()) {
+      TetaNewton = Resol.Root();
+      if (Abs(Teta - TetaNewton) > 1.e-7)
+      {
+        std::cout << "Error: Newton root is wrong for " << Teta << " \n";
+        err = 1;
+      }
+    }
+    else
+    {
+      std::cout << "Error: Newton is not done for " << Teta << " \n";
+      err = 1;
+    }
+  }
+
+  return err;
+}
 
 void QABugs::Commands_20(Draw_Interpretor& theCommands) {
   const char *group = "QABugs";
@@ -2559,6 +2618,7 @@ void QABugs::Commands_20(Draw_Interpretor& theCommands) {
                   "\n\t\t: Check interface for reading BRep from memory.",
                   __FILE__, OCC28887, group);
   theCommands.Add("OCC28131", "OCC28131 name: creates face problematic for offset", __FILE__, OCC28131, group);
+  theCommands.Add("OCC29289", "OCC29289 : searching trigonometric root by Newton iterations", __FILE__, OCC29289, group);
 
   return;
 }
index dc40792..2be57d8 100755 (executable)
@@ -121,6 +121,7 @@ math_SVD.lxx
 math_TrigonometricFunctionRoots.cxx
 math_TrigonometricFunctionRoots.hxx
 math_TrigonometricFunctionRoots.lxx
+math_TrigonometricEquationFunction.hxx
 math_Uzawa.cxx
 math_Uzawa.hxx
 math_Uzawa.lxx
diff --git a/src/math/math_TrigonometricEquationFunction.hxx b/src/math/math_TrigonometricEquationFunction.hxx
new file mode 100644 (file)
index 0000000..e7bb96d
--- /dev/null
@@ -0,0 +1,75 @@
+// Created on: 1991-05-13
+// Created by: Laurent Painnot
+// Copyright (c) 1991-1999 Matra Datavision
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _math_TrigonometricEquationFunction_HeaderFile
+#define _math_TrigonometricEquationFunction_HeaderFile
+
+#include <math_FunctionWithDerivative.hxx>
+
+//! This is function, which corresponds trigonometric equation
+//! a*Cos(x)*Cos(x) + 2*b*Cos(x)*Sin(x) + c*Cos(x) + d*Sin(x) + e = 0
+//! See class math_TrigonometricFunctionRoots
+class math_TrigonometricEquationFunction : public math_FunctionWithDerivative
+{
+  Standard_Real myAA;
+  Standard_Real myBB;
+  Standard_Real myCC;
+  Standard_Real myDD;
+  Standard_Real myEE;
+
+public:
+
+  math_TrigonometricEquationFunction(const Standard_Real A,
+    const Standard_Real B,
+    const Standard_Real C,
+    const Standard_Real D,
+    const Standard_Real E)
+    : myAA(A), myBB(B), myCC(C), myDD(D), myEE(E)
+  {
+  }
+
+  Standard_Boolean Value(const Standard_Real X, Standard_Real& F) {
+    Standard_Real CN = cos(X), SN = sin(X);
+    //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
+    F = CN*(myAA*CN + (myBB + myBB)*SN + myCC) + myDD*SN + myEE;
+    return Standard_True;
+  }
+
+  Standard_Boolean Derivative(const Standard_Real X, Standard_Real& D) {
+    Standard_Real CN = Cos(X), SN = Sin(X);
+    //-- D = -2*AA*CN*SN+2*BB*(CN*CN-SN*SN)-CC*SN+DD*CN;
+    D = -myAA*CN*SN + myBB*(CN*CN - SN*SN);
+    D += D;
+    D += -myCC*SN + myDD*CN;
+    return Standard_True;
+  }
+
+  Standard_Boolean Values(const Standard_Real X, Standard_Real& F, Standard_Real& D) {
+    Standard_Real CN = Cos(X), SN = Sin(X);
+    //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
+    //-- D = -2*AA*CN*SN+2*BB*(CN*CN-SN*SN)-CC*SN+DD*CN;
+    Standard_Real AACN = myAA*CN;
+    Standard_Real BBSN = myBB*SN;
+
+    F = AACN*CN + BBSN*(CN + CN) + myCC*CN + myDD*SN + myEE;
+    D = -AACN*SN + myBB*(CN*CN - SN*SN);
+    D += D;
+    D += -myCC*SN + myDD*CN;
+    return Standard_True;
+  }
+};
+
+#endif // _math_TrigonometricEquationFunction_HeaderFile
index 3d5b991..ff808ec 100644 (file)
 //#endif
 
 #include <math_TrigonometricFunctionRoots.hxx>
+#include <math_TrigonometricEquationFunction.hxx>
 #include <math_DirectPolynomialRoots.hxx>
 #include <Standard_OutOfRange.hxx>
 #include <math_FunctionWithDerivative.hxx>
 #include <math_NewtonFunctionRoot.hxx>
 #include <Precision.hxx>
 
-
-class MyTrigoFunction: public math_FunctionWithDerivative
-{
-  Standard_Real AA;
-  Standard_Real BB;
-  Standard_Real CC;
-  Standard_Real DD;
-  Standard_Real EE;
-
- public:
-  MyTrigoFunction(const Standard_Real A,
-                  const Standard_Real B,
-                  const Standard_Real C,
-                  const Standard_Real D,
-                  const Standard_Real E)
-  : AA(A), BB(B), CC(C), DD(D), EE(E)
-  {
-  }
-
-  Standard_Boolean Value(const Standard_Real X, Standard_Real& F);
-  Standard_Boolean Derivative(const Standard_Real X, Standard_Real& D);
-  Standard_Boolean Values(const Standard_Real X, Standard_Real& F, Standard_Real& D);
-};
-
- Standard_Boolean MyTrigoFunction::Value(const Standard_Real X, Standard_Real& F) {
-   Standard_Real CN= cos(X), SN = sin(X);
-   //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
-   F=CN*(AA*CN + (BB+BB)*SN + CC) + DD*SN + EE;
-   return Standard_True;
- }
-
- Standard_Boolean MyTrigoFunction::Derivative(const Standard_Real X, Standard_Real& D) {
-   Standard_Real CN= Cos(X), SN = Sin(X);
-   //-- D = -2*AA*CN*SN+2*BB*(CN*CN-SN*SN)-CC*SN+DD*CN;
-   D = -AA*CN*SN + BB*(CN*CN-SN*SN);
-   D+=D;
-   D-=CC*SN+DD*CN;
-   return Standard_True;
- }
-
- Standard_Boolean MyTrigoFunction::Values(const Standard_Real X, Standard_Real& F, Standard_Real& D) {
-   Standard_Real CN= Cos(X), SN = Sin(X);
-   //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
-   //-- D = -2*AA*CN*SN+2*BB*(CN*CN-SN*SN)-CC*SN+DD*CN;
-   Standard_Real AACN = AA*CN;
-   Standard_Real BBSN = BB*SN;
-
-   F = AACN*CN + BBSN*(CN+CN) + CC*CN + DD*SN + EE;
-   D = -AACN*SN + BB*(CN*CN+SN*SN);
-   D+=D;
-   D+=-CC*SN+DD*CN;
-   return Standard_True;
- }
-
-
 math_TrigonometricFunctionRoots::math_TrigonometricFunctionRoots
                                 (const Standard_Real theD,
                                  const Standard_Real theE,
@@ -474,7 +420,7 @@ void math_TrigonometricFunctionRoots::Perform(const Standard_Real A,
       // Appel de Newton:
       //OCC541(apo):  Standard_Real TetaNewton=0;  
       Standard_Real TetaNewton = Teta;  
-      MyTrigoFunction MyF(A, B, C, D, E);
+      math_TrigonometricEquationFunction MyF(A, B, C, D, E);
       math_NewtonFunctionRoot Resol(MyF, X, Tol1, Eps, Nit);
       if (Resol.IsDone()) {
        TetaNewton = Resol.Root();
index aa83918..7bb0d99 100644 (file)
@@ -54,7 +54,7 @@ public:
   Standard_EXPORT math_TrigonometricFunctionRoots(const Standard_Real D, const Standard_Real E, const Standard_Real InfBound, const Standard_Real SupBound);
   
   //! Given the three coefficients c, d and e, it performs
-  //! the resolution of 2*b*cos(x)*sin(x) + d*sin(x) + e = 0.
+  //! the resolution of c*Cos(x) + d*sin(x) + e = 0.
   //! The solutions must be contained in [InfBound, SupBound].
   //! InfBound and SupBound can be set by default to 0 and 2*PI.
   Standard_EXPORT math_TrigonometricFunctionRoots(const Standard_Real C, const Standard_Real D, const Standard_Real E, const Standard_Real InfBound, const Standard_Real SupBound);
diff --git a/tests/bugs/modalg_7/bug29289 b/tests/bugs/modalg_7/bug29289
new file mode 100644 (file)
index 0000000..6bb40c8
--- /dev/null
@@ -0,0 +1,11 @@
+puts "========"
+puts "OCC29289"
+puts "========"
+puts ""
+
+################################################################
+# Wrong derivatives in math_TrigonometricFunctionRoots.cxx file.
+################################################################
+pload QAcommands
+OCC29289