1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
19 #include <Standard_TypeMismatch.hxx>
20 #include <Precision.hxx>
21 static const Standard_Real TolFactor = 1.e-12;
22 static const Standard_Real MinTol = 1.e-20;
23 static const Standard_Real MinStep = 1e-7;
25 static const Standard_Integer MaxOrder = 3;
29 /*-----------------------------------------------------------------------------
30 Fonction permettant de rechercher une distance extremale entre un point P et
31 une courbe C (en partant d'un point approche C(u0)).
32 Cette classe herite de math_FunctionWithDerivative et est utilisee par
33 les algorithmes math_FunctionRoot et math_FunctionRoots.
34 Si on note D1c et D2c les derivees premiere et seconde:
35 { F(u) = (C(u)-P).D1c(u)/ ||Du||}
36 { DF(u) = ||Du|| + (C(u)-P).D2c(u)/||Du|| - F(u)*Duu*Du/||Du||**2
39 { F(u) = (C(u)-P).D1c(u) }
40 { DF(u) = D1c(u).D1c(u) + (C(u)-P).D2c(u)
41 = ||D1c(u)|| ** 2 + (C(u)-P).D2c(u) }
42 ----------------------------------------------------------------------------*/
44 Standard_Real Extrema_FuncExtPC::SearchOfTolerance()
46 const Standard_Integer NPoint = 10;
47 const Standard_Real aStep = (myUsupremum - myUinfium)/(Standard_Real)NPoint;
49 Standard_Integer aNum = 0;
50 Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative
51 //(it is computed with using NPoint point)
55 Standard_Real u = myUinfium + aNum*aStep; //parameter for every point
59 Pnt Ptemp; //empty point (is not used below)
60 Vec VDer; // 1st derivative vector
61 Tool::D1(*((Curve*)myC), u, Ptemp, VDer);
62 Standard_Real vm = VDer.Magnitude();
66 while(++aNum < NPoint+1);
68 return Max(aMax*TolFactor,MinTol);
72 //=============================================================================
74 Extrema_FuncExtPC::Extrema_FuncExtPC():
78 myPinit = Standard_False;
79 myCinit = Standard_False;
80 myD1Init = Standard_False;
82 SubIntervalInitialize(0.0,0.0);
87 //=============================================================================
89 Extrema_FuncExtPC::Extrema_FuncExtPC (const Pnt& P,
90 const Curve& C): myU(0.), myD1f(0.)
93 myC = (Standard_Address)&C;
94 myPinit = Standard_True;
95 myCinit = Standard_True;
96 myD1Init = Standard_False;
98 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
99 Tool::LastParameter(*((Curve*)myC)));
101 switch(Tool::GetType(*((Curve*)myC)))
103 case GeomAbs_BezierCurve:
104 case GeomAbs_BSplineCurve:
105 case GeomAbs_OtherCurve:
106 myMaxDerivOrder = MaxOrder;
107 myTol = SearchOfTolerance();
115 //=============================================================================
117 void Extrema_FuncExtPC::Initialize(const Curve& C)
119 myC = (Standard_Address)&C;
120 myCinit = Standard_True;
125 SubIntervalInitialize(Tool::FirstParameter(*((Curve*)myC)),
126 Tool::LastParameter(*((Curve*)myC)));
128 switch(Tool::GetType(*((Curve*)myC)))
130 case GeomAbs_BezierCurve:
131 case GeomAbs_BSplineCurve:
132 case GeomAbs_OtherCurve:
133 myMaxDerivOrder = MaxOrder;
134 myTol = SearchOfTolerance();
143 //=============================================================================
145 void Extrema_FuncExtPC::SetPoint(const Pnt& P)
148 myPinit = Standard_True;
154 //=============================================================================
157 Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
159 if (!myPinit || !myCinit)
160 Standard_TypeMismatch::Raise("No init");
164 Tool::D1(*((Curve*)myC),myU,myPc,D1c);
165 Standard_Real Ndu = D1c.Magnitude();
167 if(myMaxDerivOrder != 0)
169 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
171 const Standard_Real DivisionFactor = 1.e-3;
173 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
176 du = myUsupremum-myUinfium;
178 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
179 //Derivative is approximated by Taylor-series
181 Standard_Integer n = 1; //Derivative order
183 Standard_Boolean IsDeriveFound;
187 V = Tool::DN(*((Curve*)myC),myU,++n);
189 IsDeriveFound = (Ndu > myTol);
191 while(!IsDeriveFound && n < myMaxDerivOrder);
197 if(myU-myUinfium < aDelta)
203 Tool::D0(*((Curve*)myC),Min(myU, u),P1);
204 Tool::D0(*((Curve*)myC),Max(myU, u),P2);
207 Standard_Real aDirFactor = V.Dot(V1);
216 //Derivative is approximated by three points
218 Pnt Ptemp; //(0,0,0)-coordinate
220 Standard_Boolean IsParameterGrown;
222 if(myU-myUinfium < 2*aDelta)
224 Tool::D0(*((Curve*)myC),myU,P1);
225 Tool::D0(*((Curve*)myC),myU+aDelta,P2);
226 Tool::D0(*((Curve*)myC),myU+2*aDelta,P3);
227 IsParameterGrown = Standard_True;
231 Tool::D0(*((Curve*)myC),myU-2*aDelta,P1);
232 Tool::D0(*((Curve*)myC),myU-aDelta,P2);
233 Tool::D0(*((Curve*)myC),myU,P3);
234 IsParameterGrown = Standard_False;
237 Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
244 Ndu = D1c.Magnitude();
245 }//(if (Ndu <= myTol)) condition
246 }//if(myMaxDerivOrder != 0)
251 cout << "+++Function Extrema_FuncExtPC::Value(...)." << endl;
252 cout << "Warning: 1st derivative is equal to zero!---"<<endl;
254 return Standard_False;
258 F = PPc.Dot(D1c)/Ndu;
259 return Standard_True;
262 //=============================================================================
264 Standard_Boolean Extrema_FuncExtPC::Derivative (const Standard_Real U, Standard_Real& D1f)
266 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
268 return Values(U,F,D1f); /* on fait appel a Values pour simplifier la
269 sauvegarde de l'etat. */
271 //=============================================================================
273 Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U, Standard_Real& F, Standard_Real& D1f)
275 if (!myPinit || !myCinit)
276 Standard_TypeMismatch::Raise("No init");
278 Pnt myPc_old = myPc, myP_old = myP;
280 if(Value(U,F) == Standard_False)
283 cout << "+++Function Extrema_FuncExtPC::Values(...)." << endl;
284 cout << "Warning: No function value found!---"<<endl;
287 myD1Init = Standard_False;
288 return Standard_False;
296 Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
298 Standard_Real Ndu = D1c.Magnitude();
299 if (Ndu <= myTol) // Cas Singulier (PMN 22/04/1998)
301 //Derivative is approximated by three points
303 //Attention: aDelta value must be greater than same value for "Value(...)"
304 // function to avoid of points' collisions.
305 const Standard_Real DivisionFactor = 0.01;
307 if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
310 du = myUsupremum-myUinfium;
312 const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
314 Standard_Real F1, F2, F3;
316 if(myU-myUinfium < 2*aDelta)
319 //const Standard_Real U1 = myU;
320 const Standard_Real U2 = myU + aDelta;
321 const Standard_Real U3 = myU + aDelta * 2.0;
323 if(!((Value(U2,F2)) && (Value(U3,F3))))
326 cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl;
327 cout << "There are many points close to singularity points "
328 "and which have zero-derivative." << endl;
329 cout << "Try to decrease aDelta variable's value. ---" << endl;
331 myD1Init = Standard_False;
332 return Standard_False;
335 //After calling of Value(...) function variable myU will be redeterminated.
336 //So we must return it previous value.
337 D1f=(-3*F1+4*F2-F3)/(2.0*aDelta);
342 const Standard_Real U1 = myU - aDelta * 2.0;
343 const Standard_Real U2 = myU - aDelta;
344 //const Standard_Real U3 = myU;
346 if(!((Value(U2,F2)) && (Value(U1,F1))))
349 cout << "+++ Function Extrema_FuncExtPC::Values(...)" << endl;
350 cout << "There are many points close to singularity points "
351 "and which have zero-derivative." << endl;
352 cout << "Try to decrease aDelta variable's value. ---" << endl;
354 myD1Init = Standard_False;
355 return Standard_False;
357 //After calling of Value(...) function variable myU will be redeterminated.
358 //So we must return it previous value.
359 D1f=(F1-4*F2+3*F3)/(2.0*aDelta);
368 D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
373 myD1Init = Standard_True;
374 return Standard_True;
376 //=============================================================================
378 Standard_Integer Extrema_FuncExtPC::GetStateNumber ()
380 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
381 mySqDist.Append(myPc.SquareDistance(myP));
382 Standard_Integer IntVal;
384 myD1Init = Standard_True;
385 Standard_Real FF, DD;
388 if (!myD1Init) IntVal = 0;
390 if (myD1f > 0.) { IntVal = 1; }
393 myIsMin.Append(IntVal);
394 myPoint.Append(POnC(myU,myPc));
397 //=============================================================================
399 Standard_Integer Extrema_FuncExtPC::NbExt () const { return mySqDist.Length(); }
400 //=============================================================================
402 Standard_Real Extrema_FuncExtPC::SquareDistance (const Standard_Integer N) const
404 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
405 return mySqDist.Value(N);
407 //=============================================================================
408 Standard_Boolean Extrema_FuncExtPC::IsMin (const Standard_Integer N) const
410 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
411 return (myIsMin.Value(N) == 1);
413 //=============================================================================
414 const POnC & Extrema_FuncExtPC::Point (const Standard_Integer N) const
416 if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
417 return myPoint.Value(N);
419 //=============================================================================
421 void Extrema_FuncExtPC::SubIntervalInitialize(const Standard_Real theUfirst, const Standard_Real theUlast)
423 myUinfium = theUfirst;
424 myUsupremum = theUlast;