OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-08-14 13:30:48 +03:00
Code Activity

0023706: Cannot project point on curve

Browse Source 
1.   Approximation of derivative (by Taylor-series and by three points).
2.   Some methods (Degree(), GetType(), D0(), D3(), DN()) are added.
3.   Getting of subInterval's boundaries.
4.   Algorithm for checking if 1st derivative is equal to zero is amended.
5.   Cases are controlled when extrema or Project point do not exist.
6.   GetNormal() function for gp_Vec2d was added.
7.   Computing of Value, D0, D1, D2 and D3 for offset curves was changed.
8.   Limitation of tolerance for derivative computing was added.
9.   Methods for computing trihedron in singularity point are added.
10. Test tests/bugs/moddata_3/bug23706 is added.
11. Restriction on the LastParameter for visualization of 3-D curves. Calling PlotCurve(...) function for last interval.
12. LProp package is modified for tangent computing in singularity point (LProp_CLProps, LProp_SLProps).
13. Added test cases for issue.
Deleting bad test cases for this fix
This commit is contained in:
nbv
2013-06-13 15:12:06 +04:00
parent 71797c62f1
commit 32ca7a5106
93 changed files with 4498 additions and 1203 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
src/Extrema/Extrema_CurveTool.cdl
View File

@@ -83,6 +83,11 @@ is
U : Real from Standard)
returns Pnt from gp;
---C++: inline
D0 (myclass; C : Curve from Adaptor3d;
U : Real from Standard;
P : out Pnt from gp);
---C++: inline
D1 (myclass; C : Curve from Adaptor3d;
U : Real from Standard;
@@ -90,10 +95,22 @@ is
V : out Vec from gp);
---C++: inline
D2 (myclass; C : Curve from Adaptor3d;
D2 (myclass; C : Curve from Adaptor3d;
U : Real from Standard;
P : out Pnt from gp;
V1, V2 : out Vec from gp);
---C++: inline
D3 (myclass; C : Curve from Adaptor3d;
U : Real from Standard;
P : out Pnt from gp;
V1, V2, V3 : out Vec from gp);
---C++: inline
DN (myclass; C : Curve from Adaptor3d;
U : Real from Standard;
P : out Pnt from gp;
V1, V2 : out Vec from gp);
N : Integer from Standard)
returns Vec from gp;
---C++: inline
Line(myclass; C : Curve from Adaptor3d) returns Lin from gp;

36
src/Extrema/Extrema_CurveTool.lxx
View File

@@ -62,6 +62,17 @@ inline gp_Pnt Extrema_CurveTool::Value(const Adaptor3d_Curve& C,
return C.Value(U);
}
//=======================================================================
//function : D0
//purpose :
//=======================================================================
inline void Extrema_CurveTool::D0(const Adaptor3d_Curve& C,
const Standard_Real U,
gp_Pnt& P)
{
C.D0(U, P);
}
//=======================================================================
//function : D1
@@ -90,6 +101,31 @@ inline void Extrema_CurveTool::D2(const Adaptor3d_Curve& C,
C.D2(U, P, V1, V2);
}
//=======================================================================
//function : D3
//purpose :
//=======================================================================
inline void Extrema_CurveTool::D3(const Adaptor3d_Curve& C,
const Standard_Real U,
gp_Pnt& P,
gp_Vec& V1,
gp_Vec& V2,
gp_Vec& V3)
{
C.D3(U, P, V1, V2, V3);
}
//=======================================================================
//function : DN
//purpose :
//=======================================================================
inline gp_Vec Extrema_CurveTool::DN(const Adaptor3d_Curve& C,
const Standard_Real U,
const Standard_Integer N)
{
return C.DN(U, N);
}
//=======================================================================

19
src/Extrema/Extrema_FuncExtCC.cdl
View File

@@ -50,7 +50,6 @@ is
---Purpose:
SetCurve (me: in out; theRank: Integer; C1: Curve1);
---C++: inline
---Purpose:
SetTolerance (me: in out; theTol: Real);
@@ -103,6 +102,14 @@ is
---Purpose: Returns a tolerance specified in the constructor
-- or in SetTolerance() method.
SubIntervalInitialize(me: in out; theUfirst, theUlast: Vector);
---Purpose: Determines of boundaries of subinterval for find of root.
SearchOfTolerance(me: in out; C: Address from Standard) returns Real from Standard;
---Purpose: Computes a Tol value. If 1st derivative of curve
-- |D1|<Tol, it is considered D1=0.
fields
myC1 : Address from Standard;
myC2 : Address from Standard;
@@ -117,4 +124,14 @@ fields
mySqDist: SequenceOfReal from TColStd;
myPoints: SeqPOnC;
myTolC1,myTolC2: Real from Standard; -- toolerance for derivate
--Supremum of search 1st non-zero derivative
myMaxDerivOrderC1, myMaxDerivOrderC2: Integer from Standard;
--boundaries of subinterval for find of root
myUinfium, myUsupremum: Real from Standard; -- C1 curve
myVinfium, myVsupremum: Real from Standard; -- C2 curve
end FuncExtCC;

739
src/Extrema/Extrema_FuncExtCC.gxx
View File

@@ -38,126 +38,690 @@ et F2 sont egales a:
{ Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
----------------------------------------------------------------------------*/
#include <Precision.hxx>
static const Standard_Real MinTol = 1.e-20;
static const Standard_Real TolFactor = 1.e-12;
static const Standard_Real MinStep = 1e-7;
static const Standard_Integer MaxOrder = 3;
//=============================================================================
Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
{
const Standard_Integer NPoint = 10;
Standard_Real aStartParam, anEndParam;
if(C==myC1)
{
aStartParam = myUinfium;
anEndParam = myUsupremum;
}
else if(C==myC2)
{
aStartParam = myVinfium;
anEndParam = myVsupremum;
}
else
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::SearchOfTolerance(...)" << endl;
cout << "Warning: No curve for tolerance computing!---"<<endl;
#endif
return MinTol;
}
const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint;
Standard_Integer aNum = 0;
Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
//(it is computed with using NPoint point)
do
{
Standard_Real u = aStartParam + aNum*aStep; //parameter for every point
if(u > anEndParam)
u = anEndParam;
Pnt Ptemp; //empty point (is not used below)
Vec VDer; // 1st derivative vector
Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
Standard_Real vm = VDer.Magnitude();
if(vm > aMax)
aMax = vm;
}
while(++aNum < NPoint+1);
return Max(aMax*TolFactor,MinTol);
}
#define Tol 1.e-20
#define delta 1.e-9
//=============================================================================
Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
{
}
{
math_Vector V1(1,2), V2(1,2);
V1(1) = 0.0;
V2(1) = 0.0;
V1(2) = 0.0;
V2(2) = 0.0;
SubIntervalInitialize(V1, V2);
myMaxDerivOrderC1 = 0;
myTolC1=MinTol;
myMaxDerivOrderC2 = 0;
myTolC2=MinTol;
}
//=============================================================================
Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
const Curve2& C2,
const Standard_Real thetol) :
myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
myTol (thetol)
{
}
{
math_Vector V1(1,2), V2(1,2);
V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
SubIntervalInitialize(V1, V2);
switch(Tool1::GetType(*((Curve1*)myC1)))
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC1 = MaxOrder;
myTolC1 = SearchOfTolerance((Standard_Address)&C1);
break;
default:
myMaxDerivOrderC1 = 0;
myTolC1=MinTol;
break;
}
switch(Tool2::GetType(*((Curve2*)myC2)))
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC2 = MaxOrder;
myTolC2 = SearchOfTolerance((Standard_Address)&C2);
break;
default:
myMaxDerivOrderC2 = 0;
myTolC2=MinTol;
break;
}
}
Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV,
math_Vector& F)
{
//=============================================================================
void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
{
Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
if (theRank == 1)
{
myC1 = (Standard_Address)&C;
switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC1 = MaxOrder;
myTolC1 = SearchOfTolerance((Standard_Address)&C);
break;
default:
myMaxDerivOrderC1 = 0;
myTolC1=MinTol;
break;
}
}
else if (theRank == 2)
{
myC2 = (Standard_Address)&C;
switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType())
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrderC2 = MaxOrder;
myTolC2 = SearchOfTolerance((Standard_Address)&C);
break;
default:
myMaxDerivOrderC2 = 0;
myTolC2=MinTol;
break;
}
}
}
//=============================================================================
Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F)
{
myU = UV(1);
myV = UV(2);
Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);
Vec P1P2 (myP1,myP2);
Standard_Real Ndu = myDu.Magnitude();
if (Ndu <= Tol) {
Pnt P1, P2;
P1 = Tool1::Value(*((Curve1*)myC1), myU-delta);
P2 = Tool1::Value(*((Curve1*)myC1), myU+delta);
Vec V(P1,P2);
myDu = V;
Ndu = myDu.Magnitude();
if (Ndu <= Tol) {
return Standard_False;
if(myMaxDerivOrderC1 != 0)
{
if (Ndu <= myTolC1)
{
//Derivative is approximated by Taylor-series
const Standard_Real DivisionFactor = 1.e-3;
Standard_Real du;
if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
du = 0.0;
else
du = myUsupremum-myUinfium;
const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
Standard_Integer n = 1; //Derivative order
Vec V;
Standard_Boolean IsDeriveFound;
do
{
V = Tool1::DN(*((Curve1*)myC1),myU,++n);
Ndu = V.Magnitude();
IsDeriveFound = (Ndu > myTolC1);
}
while(!IsDeriveFound && n < myMaxDerivOrderC1);
if(IsDeriveFound)
{
Standard_Real u;
if(myU-myUinfium < aDelta)
u = myU+aDelta;
else
u = myU-aDelta;
Pnt P1, P2;
Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1);
Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2);
Vec V1(P1,P2);
Standard_Real aDirFactor = V.Dot(V1);
if(aDirFactor < 0.0)
myDu = -V;
else
myDu = V;
}//if(IsDeriveFound)
else
{
//Derivative is approximated by three points
Pnt Ptemp; //(0,0,0)-coordinate
Pnt P1, P2, P3;
Standard_Boolean IsParameterGrown;
if(myU-myUinfium < 2*aDelta)
{
Tool1::D0(*((Curve1*)myC1),myU,P1);
Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2);
Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3);
IsParameterGrown = Standard_True;
}
else
{
Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1);
Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2);
Tool1::D0(*((Curve1*)myC1),myU,P3);
IsParameterGrown = Standard_False;
}
Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
if(IsParameterGrown)
myDu=-3*V1+4*V2-V3;
else
myDu=V1-4*V2+3*V3;
}//else of if(IsDeriveFound)
Ndu = myDu.Magnitude();
}//if (Ndu <= myTolC1) condition
}//if(myMaxDerivOrder != 0)
if (Ndu <= MinTol)
{
#ifdef DEB
cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
cout << "Warning: 1st derivative of C1 is equal to zero!---"<<endl;
#endif
return Standard_False;
}
}
Standard_Real Ndv = myDv.Magnitude();
if (Ndv <= Tol) {
// Traitement des singularite, on approche la Tangente
// par une corde
Pnt P1, P2;
P1 = Tool2::Value(*((Curve2*)myC2), myV-delta);
P2 = Tool2::Value(*((Curve2*)myC2), myV+delta);
Vec V(P1,P2);
myDv = V;
Ndv = myDv.Magnitude();
if (Ndv <= Tol) {
return Standard_False;
if(myMaxDerivOrderC2 != 0)
{
if (Ndv <= myTolC2)
{
const Standard_Real DivisionFactor = 1.e-3;
Standard_Real dv;
if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
dv = 0.0;
else
dv = myVsupremum-myVinfium;
const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep);
//Derivative is approximated by Taylor-series
Standard_Integer n = 1; //Derivative order
Vec V;
Standard_Boolean IsDeriveFound;
do
{
V = Tool2::DN(*((Curve2*)myC2),myV,++n);
Ndv = V.Magnitude();
IsDeriveFound = (Ndv > myTolC2);
}
while(!IsDeriveFound && n < myMaxDerivOrderC2);
if(IsDeriveFound)
{
Standard_Real v;
if(myV-myVinfium < aDelta)
v = myV+aDelta;
else
v = myV-aDelta;
Pnt P1, P2;
Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1);
Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2);
Vec V1(P1,P2);
Standard_Real aDirFactor = V.Dot(V1);
if(aDirFactor < 0.0)
myDv = -V;
else
myDv = V;
}//if(IsDeriveFound)
else
{
//Derivative is approximated by three points
Pnt Ptemp; //(0,0,0)-coordinate
Pnt P1, P2, P3;
Standard_Boolean IsParameterGrown;
if(myV-myVinfium < 2*aDelta)
{
Tool2::D0(*((Curve2*)myC2),myV,P1);
Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2);
Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3);
IsParameterGrown = Standard_True;
}
else
{
Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1);
Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2);
Tool2::D0(*((Curve2*)myC2),myV,P3);
IsParameterGrown = Standard_False;
}
Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
if(IsParameterGrown)
myDv=-3*V1+4*V2-V3;
else
myDv=V1-4*V2+3*V3;
}//else of if(IsDeriveFound)
Ndv = myDv.Magnitude();
}//if (Ndv <= myTolC2)
}//if(myMaxDerivOrder != 0)
if (Ndv <= MinTol)
{
#ifdef DEB
cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
cout << "1st derivative of C2 is equal to zero!---"<<endl;
#endif
return Standard_False;
}
}
F(1) = P1P2.Dot(myDu)/Ndu;
F(2) = P1P2.Dot(myDv)/Ndv;
return Standard_True;
}
}
//=============================================================================
Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV,
math_Matrix& Df)
{
{
math_Vector F(1,2);
return Values(UV,F,Df);
}
}
//=============================================================================
Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV,
math_Vector& F,
math_Matrix& Df)
{
math_Vector& F,
math_Matrix& Df)
{
myU = UV(1);
myV = UV(2);
Vec Duu, Dvv;
Tool1::D2(*((Curve1*)myC1), myU,myP1,myDu,Duu);
Tool2::D2(*((Curve2*)myC2), myV,myP2,myDv,Dvv);
if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv
{
#ifdef DEB
cout << "+++Standard_Boolean Extrema_FuncExtCC::Values(...)." << endl;
cout << "Warning: No function value found!---"<<endl;
#endif
return Standard_False;
}//if(Value(UV, F) == Standard_False)
Vec Du, Dv, Duu, Dvv;
Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu);
Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv);
//Calling of "Value(...)" function change class member values.
//After running it is necessary to return to previous values.
const Standard_Real myU_old = myU, myV_old = myV;
const Pnt myP1_old = myP1, myP2_old = myP2;
const Vec myDu_old = myDu, myDv_old = myDv;
//Attention: aDelta value must be greater than same value for "Value(...)"
// function to avoid of points' collisions.
const Standard_Real DivisionFactor = 0.01;
Standard_Real du;
if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
du = 0.0;
else
du = myUsupremum-myUinfium;
const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
Standard_Real dv;
if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
dv = 0.0;
else
dv = myVsupremum-myVinfium;
const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
Vec P1P2 (myP1,myP2);
Standard_Real Ndu = myDu.Magnitude();
if (Ndu <= Tol) {
Pnt P1, P2;
Vec V1;
Tool1::D1(*((Curve1*)myC1),myU+delta, P2, Duu);
Tool1::D1(*((Curve1*)myC1),myU-delta, P1, V1);
Vec V(P1,P2);
myDu = V;
Duu -= V1;
Ndu = myDu.Magnitude();
if (Ndu <= Tol) {
return Standard_False;
}
}
Standard_Real Ndv = myDv.Magnitude();
if (Ndv <= Tol) {
Pnt P1, P2;
Vec V1;
Tool2::D1(*((Curve2*)myC2),myV+delta, P2, Dvv);
Tool2::D1(*((Curve2*)myC2),myV-delta, P1, V1);
Vec V(P1,P2);
myDv = V;
Dvv -= V1;
Ndv = myDv.Magnitude();
if (Ndv <= Tol) {
return Standard_False;
}
}
F(1) = P1P2.Dot(myDu)/Ndu;
F(2) = P1P2.Dot(myDv)/Ndv;
if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
{
//Derivative is approximated by three points
math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
Standard_Real F1, F2, F3;
/////////////////////////// Search of DF1_u derivative (begin) ///////////////////
if(myU-myUinfium < 2*aDeltaU)
{
F1=F(1);
math_Vector UV2(1,2), UV3(1,2);
UV2(1)=myU+aDeltaU;
UV2(2)=myV;
UV3(1)=myU+2*aDeltaU;
UV3(2)=myV;
if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F2 = FF2(1);
F3 = FF3(1);
Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
}//if(myU-myUinfium < 2*aDeltaU)
else
{
F3 = F(1);
math_Vector UV2(1,2), UV1(1,2);
UV2(1)=myU-aDeltaU;
UV2(2)=myV;
UV1(1)=myU-2*aDeltaU;
UV1(2)=myV;
if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F1 = FF1(1);
F2 = FF2(1);
Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
}//else of if(myU-myUinfium < 2*aDeltaU) condition
/////////////////////////// Search of DF1_u derivative (end) ///////////////////
//Return to previous values
myU = myU_old;
myV = myV_old;
/////////////////////////// Search of DF1_v derivative (begin) ///////////////////
if(myV-myVinfium < 2*aDeltaV)
{
F1=F(1);
math_Vector UV2(1,2), UV3(1,2);
UV2(1)=myU;
UV2(2)=myV+aDeltaV;
UV3(1)=myU;
UV3(2)=myV+2*aDeltaV;
if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F2 = FF2(1);
F3 = FF3(1);
Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
}//if(myV-myVinfium < 2*aDeltaV)
else
{
F3 = F(1);
math_Vector UV2(1,2), UV1(1,2);
UV2(1)=myU;
UV2(2)=myV-aDeltaV;
UV1(1)=myU;
UV1(2)=myV-2*aDeltaV;
if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F1 = FF1(1);
F2 = FF2(1);
Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
}//else of if(myV-myVinfium < 2*aDeltaV)
/////////////////////////// Search of DF1_v derivative (end) ///////////////////
//Return to previous values
myU = myU_old;
myV = myV_old;
myP1 = myP1_old, myP2 = myP2_old;
myDu = myDu_old, myDv = myDv_old;
}//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
else
{
const Standard_Real Ndu = myDu.Magnitude();
Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
Df(1,2) = myDv.Dot(myDu)/Ndu;
}//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
{
//Derivative is approximated by three points
math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
Standard_Real F1, F2, F3;
/////////////////////////// Search of DF2_v derivative (begin) ///////////////////
if(myV-myVinfium < 2*aDeltaV)
{
F1=F(2);
math_Vector UV2(1,2), UV3(1,2);
UV2(1)=myU;
UV2(2)=myV+aDeltaV;
UV3(1)=myU;
UV3(2)=myV+2*aDeltaV;
if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F2 = FF2(2);
F3 = FF3(2);
Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
}//if(myV-myVinfium < 2*aDeltaV)
else
{
F3 = F(2);
math_Vector UV2(1,2), UV1(1,2);
UV2(1)=myU;
UV2(2)=myV-aDeltaV;
UV1(1)=myU;
UV1(2)=myV-2*aDeltaV;
if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F1 = FF1(2);
F2 = FF2(2);
Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
}//else of if(myV-myVinfium < 2*aDeltaV)
/////////////////////////// Search of DF2_v derivative (end) ///////////////////
//Return to previous values
myU = myU_old;
myV = myV_old;
/////////////////////////// Search of DF2_u derivative (begin) ///////////////////
if(myU-myUinfium < 2*aDeltaU)
{
F1=F(2);
math_Vector UV2(1,2), UV3(1,2);
UV2(1)=myU+aDeltaU;
UV2(2)=myV;
UV3(1)=myU+2*aDeltaU;
UV3(2)=myV;
if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F2 = FF2(2);
F3 = FF3(2);
Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
}//if(myU-myUinfium < 2*aDelta)
else
{
F3 = F(2);
math_Vector UV2(1,2), UV1(1,2);
UV2(1)=myU-aDeltaU;
UV2(2)=myV;
UV1(1)=myU-2*aDeltaU;
UV1(2)=myV;
if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
return Standard_False;
}
F1 = FF1(2);
F2 = FF2(2);
Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
}//else of if(myU-myUinfium < 2*aDeltaU)
/////////////////////////// Search of DF2_u derivative (end) ///////////////////
//Return to previous values
myU = myU_old;
myV = myV_old;
myP1 = myP1_old;
myP2 = myP2_old;
myDu = myDu_old;
myDv = myDv_old;
}//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
else
{
Standard_Real Ndv = myDv.Magnitude();
Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
Df(2,1) = -myDu.Dot(myDv)/Ndv;
}//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
Df(1,2) = myDv.Dot(myDu)/Ndu;
Df(2,1) = -myDu.Dot(myDv)/Ndv;
Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
return Standard_True;
}
}//end of function
//=============================================================================
Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
@@ -166,11 +730,11 @@ Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
Vec P1P2 (myP1, myP2);
Standard_Real mod = Du.Magnitude();
if(mod > Tol) {
if(mod > myTolC1) {
Du /= mod;
}
mod = Dv.Magnitude();
if(mod > Tol) {
if(mod > myTolC2) {
Dv /= mod;
}
@@ -192,9 +756,12 @@ void Extrema_FuncExtCC::Points (const Standard_Integer N,
}
//=============================================================================
void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
const math_Vector& theSupBound)
{
myUinfium = theInfBound(1);
myUsupremum = theSupBound(1);
myVinfium = theInfBound(2);
myVsupremum = theSupBound(2);
}
//=============================================================================

10
src/Extrema/Extrema_FuncExtCC.lxx
View File

@@ -16,16 +16,6 @@
// and conditions governing the rights and limitations under the License.
//=============================================================================
inline void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
{
Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
if (theRank == 1) {myC1 = (Standard_Address)&C;}
else {myC2 = (Standard_Address)&C;}
}
//=============================================================================
inline void Extrema_FuncExtCC::SetTolerance (const Standard_Real theTol)
{
myTol = theTol;

17
src/Extrema/Extrema_FuncExtPC.cdl
View File

@@ -94,6 +94,14 @@ is
TypeMismatch from Standard;
-- if N < 1 or N > NbExt(me).
SubIntervalInitialize(me: in out; theUfirst, theUlast: Real from Standard);
---Purpose: Determines boundaries of subinterval for find of root.
SearchOfTolerance(me: in out) returns Real from Standard;
---Purpose: Computes a Tol value. If 1st derivative of curve
-- |D1|<Tol, it is considered D1=0.
fields
myP : Pnt;
myC : Address from Standard;
@@ -108,5 +116,14 @@ fields
myPinit: Boolean;
myCinit: Boolean;
myD1Init: Boolean;
myTol: Real from Standard; -- toolerance for derivate
--Supremum of search 1st non-zero derivative
myMaxDerivOrder: Integer from Standard;
--boundaries of subinterval for find of root
myUinfium, myUsupremum: Real from Standard;
end FuncExtPC;

311
src/Extrema/Extrema_FuncExtPC.gxx
View File

@@ -18,8 +18,13 @@
#include <Standard_TypeMismatch.hxx>
#include <Precision.hxx>
#define delta 1.e-9
#define Tol 1.e-20
static const Standard_Real TolFactor = 1.e-12;
static const Standard_Real MinTol = 1.e-20;
static const Standard_Real MinStep = 1e-7;
static const Standard_Integer MaxOrder = 3;
/*-----------------------------------------------------------------------------
Fonction permettant de rechercher une distance extremale entre un point P et
@@ -36,6 +41,35 @@ les algorithmes math_FunctionRoot et math_FunctionRoots.
= ||D1c(u)|| ** 2 + (C(u)-P).D2c(u) }
----------------------------------------------------------------------------*/
Standard_Real Extrema_FuncExtPC::SearchOfTolerance()
{
const Standard_Integer NPoint = 10;
const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint;
Standard_Integer aNum = 0;
Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
//(it is computed with using NPoint point)
do
{
Standard_Real u = myUinfium + aNum*aStep; //parameter for every point
if(u > myUsupremum)
u = myUsupremum;
Pnt Ptemp; //empty point (is not used below)
Vec VDer; // 1st derivative vector
Tool::D1(*((Curve*)myC), u, Ptemp, VDer);
Standard_Real vm = VDer.Magnitude();
if(vm > aMax)
aMax = vm;
}
while(++aNum < NPoint+1);
return Max(aMax*TolFactor,MinTol);
}
//=============================================================================
Extrema_FuncExtPC::Extrema_FuncExtPC():
myU(0.),
@@ -44,31 +78,67 @@ myD1f(0.)
myPinit = Standard_False;
myCinit = Standard_False;
myD1Init = Standard_False;
SubIntervalInitialize(0.0,0.0);
myMaxDerivOrder = 0;
myTol=MinTol;
}
//=============================================================================
Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P,
const Curve& C):
myU(0.),
myD1f(0.)
{
const Curve& C): myU(0.), myD1f(0.)
{
myP = P;
myC = (Standard_Address)&C;
myPinit = Standard_True;
myCinit = Standard_True;
myD1Init = Standard_False;
}
SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
Tool::LastParameter(*((Curve*)myC)));
switch(Tool::GetType(*((Curve*)myC)))
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrder = MaxOrder;
myTol = SearchOfTolerance();
break;
default:
myMaxDerivOrder = 0;
myTol=MinTol;
break;
}
}
//=============================================================================
void Extrema_FuncExtPC::Initialize(const Curve& C)
{
{
myC = (Standard_Address)&C;
myCinit = Standard_True;
myPoint.Clear();
mySqDist.Clear();
myIsMin.Clear();
}
SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
Tool::LastParameter(*((Curve*)myC)));
switch(Tool::GetType(*((Curve*)myC)))
{
case GeomAbs_BezierCurve:
case GeomAbs_BSplineCurve:
case GeomAbs_OtherCurve:
myMaxDerivOrder = MaxOrder;
myTol = SearchOfTolerance();
break;
default:
myMaxDerivOrder = 0;
myTol=MinTol;
break;
}
}
//=============================================================================
@@ -83,30 +153,106 @@ void Extrema_FuncExtPC::SetPoint(const Pnt& P)
//=============================================================================
Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
{
if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
if (!myPinit || !myCinit)
Standard_TypeMismatch::Raise("No init");
myU = U;
Vec D1c;
Tool::D1(*((Curve*)myC),myU,myPc,D1c);
Standard_Real Ndu = D1c.Magnitude();
if (Ndu <= Tol) { // Cas Singulier (PMN 22/04/1998)
Pnt P1, P2;
P2 = Tool::Value(*((Curve*)myC),myU + delta);
P1 = Tool::Value(*((Curve*)myC),myU - delta);
Vec V(P1,P2);
D1c = V;
Ndu = D1c.Magnitude();
if (Ndu <= Tol) {
Vec aD2;
Tool::D2(*((Curve*)myC),myU,myPc,D1c,aD2);
Ndu = aD2.Magnitude();
if(Ndu <= Tol)
return Standard_False;
D1c = aD2;
if(myMaxDerivOrder != 0)
{
if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
{
const Standard_Real DivisionFactor = 1.e-3;
Standard_Real du;
if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
du = 0.0;
else
du = myUsupremum-myUinfium;
const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
//Derivative is approximated by Taylor-series
Standard_Integer n = 1; //Derivative order
Vec V;
Standard_Boolean IsDeriveFound;
do
{
V = Tool::DN(*((Curve*)myC),myU,++n);
Ndu = V.Magnitude();
IsDeriveFound = (Ndu > myTol);
}
while(!IsDeriveFound && n < myMaxDerivOrder);
if(IsDeriveFound)
{
Standard_Real u;
if(myU-myUinfium < aDelta)
u = myU+aDelta;
else
u = myU-aDelta;
Pnt P1, P2;
Tool::D0(*((Curve*)myC),Min(myU, u),P1);
Tool::D0(*((Curve*)myC),Max(myU, u),P2);
Vec V1(P1,P2);
Standard_Real aDirFactor = V.Dot(V1);
if(aDirFactor < 0.0)
D1c = -V;
else
D1c = V;
}//if(IsDeriveFound)
else
{
//Derivative is approximated by three points
Pnt Ptemp; //(0,0,0)-coordinate
Pnt P1, P2, P3;
Standard_Boolean IsParameterGrown;
if(myU-myUinfium < 2*aDelta)
{
Tool::D0(*((Curve*)myC),myU,P1);
Tool::D0(*((Curve*)myC),myU+aDelta,P2);
Tool::D0(*((Curve*)myC),myU+2*aDelta,P3);
IsParameterGrown = Standard_True;
}
else
{
Tool::D0(*((Curve*)myC),myU-2*aDelta,P1);
Tool::D0(*((Curve*)myC),myU-aDelta,P2);
Tool::D0(*((Curve*)myC),myU,P3);
IsParameterGrown = Standard_False;
}
Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
if(IsParameterGrown)
D1c=-3*V1+4*V2-V3;
else
D1c=V1-4*V2+3*V3;
}
Ndu = D1c.Magnitude();
}//(if (Ndu <= myTol)) condition
}//if(myMaxDerivOrder != 0)
if (Ndu <= MinTol)
{
#ifdef DEB
cout << "+++Function Extrema_FuncExtPC::Value(...)." << endl;
cout << "Warning: 1st derivative is equal to zero!---"<<endl;
#endif
return Standard_False;
}
}
Vec PPc (myP,myPc);
F = PPc.Dot(D1c)/Ndu;
@@ -125,36 +271,104 @@ Standard_Boolean Extrema_FuncExtPC::Derivative (const Standard_Real U, Standard_
//=============================================================================
Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U, Standard_Real& F, Standard_Real& D1f)
{
if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
{
if (!myPinit || !myCinit)
Standard_TypeMismatch::Raise("No init");
Pnt myPc_old = myPc, myP_old = myP;
if(Value(U,F) == Standard_False)
{
#ifdef DEB
cout << "+++Function Extrema_FuncExtPC::Values(...)." << endl;
cout << "Warning: No function value found!---"<<endl;
#endif
myD1Init = Standard_False;
return Standard_False;
}
myU = U;
myPc = myPc_old;
myP = myP_old;
Vec D1c,D2c;
Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
Standard_Real Ndu = D1c.Magnitude();
if (Ndu <= Tol) {// Cas Singulier (PMN 22/04/1998)
Pnt P1, P2;
Vec V1;
Tool::D1(*((Curve*)myC),myU+delta, P2, V1);
Tool::D1(*((Curve*)myC),myU-delta, P1, D2c);
Vec V(P1,P2);
D1c = V;
D2c -= V1;
Ndu = D1c.Magnitude();
if (Ndu <= Tol) {
myD1Init = Standard_False;
return Standard_False;
if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
{
//Derivative is approximated by three points
//Attention: aDelta value must be greater than same value for "Value(...)"
// function to avoid of points' collisions.
const Standard_Real DivisionFactor = 0.01;
Standard_Real du;
if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
du = 0.0;
else
du = myUsupremum-myUinfium;
const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
Standard_Real F1, F2, F3;
if(myU-myUinfium < 2*aDelta)
{
F1=F;
const Standard_Real U1 = myU, U2 = myU+aDelta, U3 = myU+2*aDelta;
if(!((Value(U2,F2)) && (Value(U3,F3))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
myD1Init = Standard_False;
return Standard_False;
}
//After calling of Value(...) function variable myU will be redeterminated.
//So we must return it previous value.
D1f=(-3*F1+4*F2-F3)/(2.0*aDelta);
}
else
{
F3 = F;
const Standard_Real U1 = myU-2*aDelta, U2 = myU-aDelta, U3 = myU;
if(!((Value(U2,F2)) && (Value(U1,F1))))
{
#ifdef DEB
cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl;
cout << "There are many points close to singularity points "
"and which have zero-derivative." << endl;
cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
myD1Init = Standard_False;
return Standard_False;
}
//After calling of Value(...) function variable myU will be redeterminated.
//So we must return it previous value.
D1f=(F1-4*F2+3*F3)/(2.0*aDelta);
}
myU = U;
myPc = myPc_old;
myP = myP_old;
}
else
{
Vec PPc (myP,myPc);
D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
}
}
Vec PPc (myP,myPc);
F = PPc.Dot(D1c)/Ndu;
D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
myD1f = D1f;
myD1Init = Standard_True;
return Standard_True;
}
}
//=============================================================================
Standard_Integer Extrema_FuncExtPC::GetStateNumber ()
@@ -200,3 +414,8 @@ POnC Extrema_FuncExtPC::Point (const Standard_Integer N) const
}
//=============================================================================
void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast)
{
myUinfium = theUfirst;
myUsupremum = theUlast;
}

19
src/Extrema/Extrema_GExtCC2d.gxx
View File

@@ -40,6 +40,20 @@
#include <Precision.hxx>
//=======================================================================
//function : IsParallelDot
//purpose : This function returns True if angle between <theV1> and
//<theV2> vectors or between <theV1> and <-theV2> vectors is less than
//AngTol.
//=======================================================================
static Standard_Boolean IsParallelDot( gp_Vec2d theV1,
gp_Vec2d theV2,
Standard_Real AngTol)
{
return Abs(theV1.Dot(theV2)) >=
theV1.Magnitude()*theV2.Magnitude()*cos(AngTol);
}
Extrema_GExtCC2d::Extrema_GExtCC2d() {}
@@ -535,7 +549,8 @@ void Extrema_GExtCC2d::Results(const Extrema_ECC2d& AlgExt,
gp_Vec2d v1, v2;
Tool1::D1(C1,U,p, v1);
Tool2::D1(*((Curve2*)myC),U2,p, v2);
if (v1.IsParallel(v2, Precision::Angular())) {
if (IsParallelDot(v1, v2, Precision::Angular()))
{
mynbext++;
Val = AlgExt.SquareDistance(i);
P1.SetValues(U, P1.Value());
@@ -543,7 +558,7 @@ void Extrema_GExtCC2d::Results(const Extrema_ECC2d& AlgExt,
mySqDist.Append(Val);
mypoints.Append(P1);
mypoints.Append(P2);
}
}
// modified by NIZHNY-EAP Thu Jan 27 16:41:00 2000 ___END___
}
}

6
src/Extrema/Extrema_GExtPC.gxx
View File

@@ -131,8 +131,9 @@ void Extrema_GExtPC::Perform(const ThePoint& P)
if (myuinf >= anInfToCheck) anInfToCheck = myuinf;
if (myusup <= aSupToCheck) aSupToCheck = myusup;
if((aSupToCheck - anInfToCheck) <= mytolu) continue;
if (i != 1) {
if (i != 1)
{
TheCurveTool::D1(*((TheCurve*)myC), myintuinf, PP, V1);
s1 = (TheVector(P, PP))*V1;
if (s1*s2 < 0.0) {
@@ -145,6 +146,7 @@ void Extrema_GExtPC::Perform(const ThePoint& P)
TheCurveTool::D1(*((TheCurve*)myC), myintusup, PP, V1);
s2 = (TheVector(P, PP))*V1;
}
IntervalPerform(P);
IntExtIsDone = IntExtIsDone || mydone;
}

2
src/Extrema/Extrema_GenExtCC.gxx
View File

@@ -184,6 +184,8 @@ b- Calcul des minima:
UVinf(2) = aCache2->TrimFirstParameter();
UVsup(1) = aCache1->TrimLastParameter();
UVsup(2) = aCache2->TrimLastParameter();
myF.SubIntervalInitialize(UVinf,UVsup);
// - des 'bords' du tableau TbDist2
for (NoV = 0; NoV <= aNbV+1; NoV++) {

1
src/Extrema/Extrema_GenExtPC.gxx
View File

@@ -171,6 +171,7 @@ Methode:
-----------------------------------------------------------------------------*/
{
myF.SetPoint(P);
myF.SubIntervalInitialize(myumin,myusup);
myDone = Standard_False;
math_FunctionRoots S (myF, myumin, myusup, mynbsample, mytolu, mytolF, mytolF);