1
0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-03 17:56:21 +03:00

0029289: Wrong derivatives in math_TrigonometricFunctionRoots.cxx file

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
This commit is contained in:
ifv 2017-11-09 17:20:10 +03:00 committed by bugmaster
parent b8bf959578
commit 18d8e3e794
8 changed files with 295 additions and 67 deletions

View 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.

View File

@ -41,7 +41,10 @@
#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);
}

View 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;
}

View File

@ -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

View File

@ -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

View File

@ -25,67 +25,13 @@
//#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();

View 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);

View File

@ -0,0 +1,11 @@
puts "========"
puts "OCC29289"
puts "========"
puts ""
################################################################
# Wrong derivatives in math_TrigonometricFunctionRoots.cxx file.
################################################################
pload QAcommands
OCC29289