Commit | Line | Data |
---|---|---|
b311480e | 1 | // Copyright (c) 1995-1999 Matra Datavision |
2 | // Copyright (c) 1999-2012 OPEN CASCADE SAS | |
3 | // | |
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. | |
8 | // | |
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. | |
11 | // | |
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. | |
18 | ||
7fd59977 | 19 | //Modified by MPS : 21-05-97 PRO 7598 |
20 | // Le nombre de solutions rendu est mySqDist.Length() et non | |
21 | // mySqDist.Length()/2 | |
22 | // ajout des valeurs absolues dans le test d'orthogonalite de | |
23 | // GetStateNumber() | |
24 | ||
25 | /*----------------------------------------------------------------------------- | |
26 | Fonctions permettant de rechercher une distance extremale entre 2 courbes | |
27 | C1 et C2 (en partant de points approches C1(u0) et C2(v0)). | |
28 | Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par | |
29 | l'algorithme math_FunctionSetRoot. | |
30 | Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont: | |
31 | { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| } | |
32 | { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| } | |
33 | Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1 | |
34 | et F2 sont egales a: | |
35 | { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2 | |
36 | { Dvf1(u,v) = Dv.Du/||Du|| | |
37 | { Duf2(u,v) = -Du.Dv/||Dv|| | |
38 | { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2 | |
39 | ||
40 | ----------------------------------------------------------------------------*/ | |
7fd59977 | 41 | |
32ca7a51 | 42 | #include <Precision.hxx> |
43 | ||
44 | ||
45 | static const Standard_Real MinTol = 1.e-20; | |
46 | static const Standard_Real TolFactor = 1.e-12; | |
47 | static const Standard_Real MinStep = 1e-7; | |
48 | ||
49 | static const Standard_Integer MaxOrder = 3; | |
50 | ||
51 | ||
52 | ||
53 | //============================================================================= | |
54 | Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C) | |
55 | { | |
56 | const Standard_Integer NPoint = 10; | |
57 | Standard_Real aStartParam, anEndParam; | |
58 | ||
59 | if(C==myC1) | |
60 | { | |
61 | aStartParam = myUinfium; | |
62 | anEndParam = myUsupremum; | |
63 | } | |
64 | else if(C==myC2) | |
65 | { | |
66 | aStartParam = myVinfium; | |
67 | anEndParam = myVsupremum; | |
68 | } | |
69 | else | |
70 | { | |
71 | #ifdef DEB | |
72 | cout << "+++ Function Extrema_FuncExtCC::SearchOfTolerance(...)" << endl; | |
73 | cout << "Warning: No curve for tolerance computing!---"<<endl; | |
74 | #endif | |
75 | return MinTol; | |
76 | } | |
77 | ||
78 | const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint; | |
79 | ||
80 | Standard_Integer aNum = 0; | |
81 | Standard_Real aMax = -Precision::Infinite(); //Maximum value of 1st derivative | |
82 | //(it is computed with using NPoint point) | |
83 | ||
84 | do | |
85 | { | |
86 | Standard_Real u = aStartParam + aNum*aStep; //parameter for every point | |
87 | if(u > anEndParam) | |
88 | u = anEndParam; | |
89 | ||
90 | Pnt Ptemp; //empty point (is not used below) | |
91 | Vec VDer; // 1st derivative vector | |
92 | Tool1::D1(*((Curve1*)C), u, Ptemp, VDer); | |
93 | Standard_Real vm = VDer.Magnitude(); | |
94 | if(vm > aMax) | |
95 | aMax = vm; | |
96 | } | |
97 | while(++aNum < NPoint+1); | |
98 | ||
99 | return Max(aMax*TolFactor,MinTol); | |
100 | } | |
7fd59977 | 101 | |
32ca7a51 | 102 | //============================================================================= |
7fd59977 | 103 | Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol) |
32ca7a51 | 104 | { |
105 | math_Vector V1(1,2), V2(1,2); | |
106 | V1(1) = 0.0; | |
107 | V2(1) = 0.0; | |
108 | V1(2) = 0.0; | |
109 | V2(2) = 0.0; | |
110 | SubIntervalInitialize(V1, V2); | |
111 | myMaxDerivOrderC1 = 0; | |
112 | myTolC1=MinTol; | |
113 | myMaxDerivOrderC2 = 0; | |
114 | myTolC2=MinTol; | |
115 | } | |
7fd59977 | 116 | |
32ca7a51 | 117 | //============================================================================= |
7fd59977 | 118 | Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1, |
119 | const Curve2& C2, | |
120 | const Standard_Real thetol) : | |
121 | myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2), | |
122 | myTol (thetol) | |
32ca7a51 | 123 | { |
124 | math_Vector V1(1,2), V2(1,2); | |
125 | V1(1) = Tool1::FirstParameter(*((Curve1*)myC1)); | |
126 | V2(1) = Tool1::LastParameter(*((Curve1*)myC1)); | |
127 | V1(2) = Tool2::FirstParameter(*((Curve2*)myC2)); | |
128 | V2(2) = Tool2::LastParameter(*((Curve2*)myC2)); | |
129 | SubIntervalInitialize(V1, V2); | |
130 | ||
131 | switch(Tool1::GetType(*((Curve1*)myC1))) | |
132 | { | |
133 | case GeomAbs_BezierCurve: | |
134 | case GeomAbs_BSplineCurve: | |
135 | case GeomAbs_OtherCurve: | |
136 | myMaxDerivOrderC1 = MaxOrder; | |
137 | myTolC1 = SearchOfTolerance((Standard_Address)&C1); | |
138 | break; | |
139 | default: | |
140 | myMaxDerivOrderC1 = 0; | |
141 | myTolC1=MinTol; | |
142 | break; | |
143 | } | |
144 | ||
145 | switch(Tool2::GetType(*((Curve2*)myC2))) | |
146 | { | |
147 | case GeomAbs_BezierCurve: | |
148 | case GeomAbs_BSplineCurve: | |
149 | case GeomAbs_OtherCurve: | |
150 | myMaxDerivOrderC2 = MaxOrder; | |
151 | myTolC2 = SearchOfTolerance((Standard_Address)&C2); | |
152 | break; | |
153 | default: | |
154 | myMaxDerivOrderC2 = 0; | |
155 | myTolC2=MinTol; | |
156 | break; | |
157 | } | |
158 | } | |
7fd59977 | 159 | |
32ca7a51 | 160 | //============================================================================= |
161 | void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C) | |
162 | { | |
163 | Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()") | |
7fd59977 | 164 | |
32ca7a51 | 165 | if (theRank == 1) |
166 | { | |
167 | myC1 = (Standard_Address)&C; | |
168 | switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType()) | |
169 | { | |
170 | case GeomAbs_BezierCurve: | |
171 | case GeomAbs_BSplineCurve: | |
172 | case GeomAbs_OtherCurve: | |
173 | myMaxDerivOrderC1 = MaxOrder; | |
174 | myTolC1 = SearchOfTolerance((Standard_Address)&C); | |
175 | break; | |
176 | default: | |
177 | myMaxDerivOrderC1 = 0; | |
178 | myTolC1=MinTol; | |
179 | break; | |
180 | } | |
181 | } | |
182 | else if (theRank == 2) | |
183 | { | |
184 | myC2 = (Standard_Address)&C; | |
185 | switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType()) | |
186 | { | |
187 | case GeomAbs_BezierCurve: | |
188 | case GeomAbs_BSplineCurve: | |
189 | case GeomAbs_OtherCurve: | |
190 | myMaxDerivOrderC2 = MaxOrder; | |
191 | myTolC2 = SearchOfTolerance((Standard_Address)&C); | |
192 | break; | |
193 | default: | |
194 | myMaxDerivOrderC2 = 0; | |
195 | myTolC2=MinTol; | |
196 | break; | |
197 | } | |
198 | } | |
199 | } | |
200 | ||
201 | //============================================================================= | |
202 | Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F) | |
203 | { | |
7fd59977 | 204 | myU = UV(1); |
205 | myV = UV(2); | |
0d167958 RL |
206 | Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu); |
207 | Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv); | |
32ca7a51 | 208 | |
7fd59977 | 209 | Vec P1P2 (myP1,myP2); |
210 | ||
0d167958 | 211 | Standard_Real Ndu = myDu.Magnitude(); |
32ca7a51 | 212 | |
213 | ||
214 | if(myMaxDerivOrderC1 != 0) | |
215 | { | |
216 | if (Ndu <= myTolC1) | |
217 | { | |
218 | //Derivative is approximated by Taylor-series | |
219 | const Standard_Real DivisionFactor = 1.e-3; | |
220 | Standard_Real du; | |
221 | if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst())) | |
222 | du = 0.0; | |
223 | else | |
224 | du = myUsupremum-myUinfium; | |
225 | ||
226 | const Standard_Real aDelta = Max(du*DivisionFactor,MinStep); | |
227 | ||
228 | Standard_Integer n = 1; //Derivative order | |
229 | Vec V; | |
230 | Standard_Boolean IsDeriveFound; | |
231 | ||
232 | do | |
233 | { | |
234 | V = Tool1::DN(*((Curve1*)myC1),myU,++n); | |
235 | Ndu = V.Magnitude(); | |
236 | IsDeriveFound = (Ndu > myTolC1); | |
237 | } | |
238 | while(!IsDeriveFound && n < myMaxDerivOrderC1); | |
239 | ||
240 | if(IsDeriveFound) | |
241 | { | |
242 | Standard_Real u; | |
243 | ||
244 | if(myU-myUinfium < aDelta) | |
245 | u = myU+aDelta; | |
246 | else | |
247 | u = myU-aDelta; | |
248 | ||
249 | Pnt P1, P2; | |
250 | Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1); | |
251 | Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2); | |
252 | ||
253 | Vec V1(P1,P2); | |
254 | Standard_Real aDirFactor = V.Dot(V1); | |
255 | ||
256 | if(aDirFactor < 0.0) | |
257 | myDu = -V; | |
258 | else | |
259 | myDu = V; | |
260 | }//if(IsDeriveFound) | |
261 | else | |
262 | { | |
263 | //Derivative is approximated by three points | |
264 | ||
265 | Pnt Ptemp; //(0,0,0)-coordinate | |
266 | Pnt P1, P2, P3; | |
267 | Standard_Boolean IsParameterGrown; | |
268 | ||
269 | if(myU-myUinfium < 2*aDelta) | |
270 | { | |
271 | Tool1::D0(*((Curve1*)myC1),myU,P1); | |
272 | Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2); | |
273 | Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3); | |
274 | IsParameterGrown = Standard_True; | |
275 | } | |
276 | else | |
277 | { | |
278 | Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1); | |
279 | Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2); | |
280 | Tool1::D0(*((Curve1*)myC1),myU,P3); | |
281 | IsParameterGrown = Standard_False; | |
282 | } | |
283 | ||
284 | Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); | |
285 | ||
286 | if(IsParameterGrown) | |
287 | myDu=-3*V1+4*V2-V3; | |
288 | else | |
289 | myDu=V1-4*V2+3*V3; | |
290 | }//else of if(IsDeriveFound) | |
291 | Ndu = myDu.Magnitude(); | |
292 | }//if (Ndu <= myTolC1) condition | |
293 | }//if(myMaxDerivOrder != 0) | |
294 | ||
295 | if (Ndu <= MinTol) | |
296 | { | |
297 | #ifdef DEB | |
298 | cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl; | |
299 | cout << "Warning: 1st derivative of C1 is equal to zero!---"<<endl; | |
300 | #endif | |
301 | return Standard_False; | |
7fd59977 | 302 | } |
32ca7a51 | 303 | |
7fd59977 | 304 | |
0d167958 | 305 | Standard_Real Ndv = myDv.Magnitude(); |
32ca7a51 | 306 | |
307 | if(myMaxDerivOrderC2 != 0) | |
308 | { | |
309 | if (Ndv <= myTolC2) | |
310 | { | |
311 | const Standard_Real DivisionFactor = 1.e-3; | |
312 | Standard_Real dv; | |
313 | if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst())) | |
314 | dv = 0.0; | |
315 | else | |
316 | dv = myVsupremum-myVinfium; | |
317 | ||
318 | const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep); | |
319 | ||
320 | //Derivative is approximated by Taylor-series | |
321 | ||
322 | Standard_Integer n = 1; //Derivative order | |
323 | Vec V; | |
324 | Standard_Boolean IsDeriveFound; | |
325 | ||
326 | do | |
327 | { | |
328 | V = Tool2::DN(*((Curve2*)myC2),myV,++n); | |
329 | Ndv = V.Magnitude(); | |
330 | IsDeriveFound = (Ndv > myTolC2); | |
331 | } | |
332 | while(!IsDeriveFound && n < myMaxDerivOrderC2); | |
333 | ||
334 | if(IsDeriveFound) | |
335 | { | |
336 | Standard_Real v; | |
337 | ||
338 | if(myV-myVinfium < aDelta) | |
339 | v = myV+aDelta; | |
340 | else | |
341 | v = myV-aDelta; | |
342 | ||
343 | Pnt P1, P2; | |
344 | Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1); | |
345 | Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2); | |
346 | ||
347 | Vec V1(P1,P2); | |
348 | Standard_Real aDirFactor = V.Dot(V1); | |
349 | ||
350 | if(aDirFactor < 0.0) | |
351 | myDv = -V; | |
352 | else | |
353 | myDv = V; | |
354 | }//if(IsDeriveFound) | |
355 | else | |
356 | { | |
357 | //Derivative is approximated by three points | |
358 | ||
359 | Pnt Ptemp; //(0,0,0)-coordinate | |
360 | Pnt P1, P2, P3; | |
361 | Standard_Boolean IsParameterGrown; | |
362 | ||
363 | if(myV-myVinfium < 2*aDelta) | |
364 | { | |
365 | Tool2::D0(*((Curve2*)myC2),myV,P1); | |
366 | Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2); | |
367 | Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3); | |
368 | IsParameterGrown = Standard_True; | |
369 | } | |
370 | else | |
371 | { | |
372 | Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1); | |
373 | Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2); | |
374 | Tool2::D0(*((Curve2*)myC2),myV,P3); | |
375 | IsParameterGrown = Standard_False; | |
376 | } | |
377 | ||
378 | Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3); | |
379 | ||
380 | if(IsParameterGrown) | |
381 | myDv=-3*V1+4*V2-V3; | |
382 | else | |
383 | myDv=V1-4*V2+3*V3; | |
384 | }//else of if(IsDeriveFound) | |
385 | ||
386 | Ndv = myDv.Magnitude(); | |
387 | }//if (Ndv <= myTolC2) | |
388 | }//if(myMaxDerivOrder != 0) | |
389 | ||
390 | if (Ndv <= MinTol) | |
391 | { | |
392 | #ifdef DEB | |
393 | cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl; | |
394 | cout << "1st derivative of C2 is equal to zero!---"<<endl; | |
395 | #endif | |
396 | return Standard_False; | |
7fd59977 | 397 | } |
7fd59977 | 398 | |
0d167958 RL |
399 | F(1) = P1P2.Dot(myDu)/Ndu; |
400 | F(2) = P1P2.Dot(myDv)/Ndv; | |
7fd59977 | 401 | return Standard_True; |
32ca7a51 | 402 | } |
7fd59977 | 403 | //============================================================================= |
404 | ||
405 | Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV, | |
406 | math_Matrix& Df) | |
32ca7a51 | 407 | { |
7fd59977 | 408 | math_Vector F(1,2); |
409 | return Values(UV,F,Df); | |
32ca7a51 | 410 | } |
7fd59977 | 411 | //============================================================================= |
32ca7a51 | 412 | |
7fd59977 | 413 | Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV, |
32ca7a51 | 414 | math_Vector& F, |
415 | math_Matrix& Df) | |
416 | { | |
7fd59977 | 417 | myU = UV(1); |
418 | myV = UV(2); | |
32ca7a51 | 419 | |
420 | if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv | |
421 | { | |
422 | #ifdef DEB | |
423 | cout << "+++Standard_Boolean Extrema_FuncExtCC::Values(...)." << endl; | |
424 | cout << "Warning: No function value found!---"<<endl; | |
425 | #endif | |
426 | return Standard_False; | |
427 | }//if(Value(UV, F) == Standard_False) | |
428 | ||
429 | Vec Du, Dv, Duu, Dvv; | |
430 | Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu); | |
431 | Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv); | |
432 | ||
433 | //Calling of "Value(...)" function change class member values. | |
434 | //After running it is necessary to return to previous values. | |
435 | const Standard_Real myU_old = myU, myV_old = myV; | |
436 | const Pnt myP1_old = myP1, myP2_old = myP2; | |
437 | const Vec myDu_old = myDu, myDv_old = myDv; | |
438 | ||
439 | ||
440 | //Attention: aDelta value must be greater than same value for "Value(...)" | |
441 | // function to avoid of points' collisions. | |
442 | ||
443 | const Standard_Real DivisionFactor = 0.01; | |
444 | ||
445 | Standard_Real du; | |
446 | if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst())) | |
447 | du = 0.0; | |
448 | else | |
449 | du = myUsupremum-myUinfium; | |
450 | ||
451 | const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep); | |
452 | ||
453 | Standard_Real dv; | |
454 | if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst())) | |
455 | dv = 0.0; | |
456 | else | |
457 | dv = myVsupremum-myVinfium; | |
458 | ||
459 | const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); | |
7fd59977 | 460 | |
461 | Vec P1P2 (myP1,myP2); | |
32ca7a51 | 462 | |
463 | if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) | |
464 | { | |
465 | //Derivative is approximated by three points | |
7fd59977 | 466 | |
32ca7a51 | 467 | math_Vector FF1(1,2), FF2(1,2), FF3(1,2); |
468 | Standard_Real F1, F2, F3; | |
7fd59977 | 469 | |
32ca7a51 | 470 | /////////////////////////// Search of DF1_u derivative (begin) /////////////////// |
471 | if(myU-myUinfium < 2*aDeltaU) | |
472 | { | |
473 | F1=F(1); | |
474 | math_Vector UV2(1,2), UV3(1,2); | |
475 | UV2(1)=myU+aDeltaU; | |
476 | UV2(2)=myV; | |
477 | UV3(1)=myU+2*aDeltaU; | |
478 | UV3(2)=myV; | |
479 | if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) | |
480 | { | |
481 | #ifdef DEB | |
482 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
483 | cout << "There are many points close to singularity points " | |
484 | "and which have zero-derivative." << endl; | |
485 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
486 | #endif | |
487 | return Standard_False; | |
488 | } | |
489 | ||
490 | F2 = FF2(1); | |
491 | F3 = FF3(1); | |
492 | ||
493 | Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU); | |
494 | }//if(myU-myUinfium < 2*aDeltaU) | |
495 | else | |
496 | { | |
497 | F3 = F(1); | |
498 | math_Vector UV2(1,2), UV1(1,2); | |
499 | UV2(1)=myU-aDeltaU; | |
500 | UV2(2)=myV; | |
501 | UV1(1)=myU-2*aDeltaU; | |
502 | UV1(2)=myV; | |
503 | ||
504 | if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) | |
505 | { | |
506 | #ifdef DEB | |
507 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
508 | cout << "There are many points close to singularity points " | |
509 | "and which have zero-derivative." << endl; | |
510 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
511 | #endif | |
512 | return Standard_False; | |
513 | } | |
514 | ||
515 | F1 = FF1(1); | |
516 | F2 = FF2(1); | |
517 | ||
518 | Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU); | |
519 | }//else of if(myU-myUinfium < 2*aDeltaU) condition | |
520 | /////////////////////////// Search of DF1_u derivative (end) /////////////////// | |
521 | ||
522 | //Return to previous values | |
523 | myU = myU_old; | |
524 | myV = myV_old; | |
525 | ||
526 | /////////////////////////// Search of DF1_v derivative (begin) /////////////////// | |
527 | if(myV-myVinfium < 2*aDeltaV) | |
528 | { | |
529 | F1=F(1); | |
530 | math_Vector UV2(1,2), UV3(1,2); | |
531 | UV2(1)=myU; | |
532 | UV2(2)=myV+aDeltaV; | |
533 | UV3(1)=myU; | |
534 | UV3(2)=myV+2*aDeltaV; | |
535 | ||
536 | if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) | |
537 | { | |
538 | #ifdef DEB | |
539 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
540 | cout << "There are many points close to singularity points " | |
541 | "and which have zero-derivative." << endl; | |
542 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
543 | #endif | |
544 | return Standard_False; | |
545 | } | |
546 | F2 = FF2(1); | |
547 | F3 = FF3(1); | |
548 | ||
549 | Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV); | |
550 | }//if(myV-myVinfium < 2*aDeltaV) | |
551 | else | |
552 | { | |
553 | F3 = F(1); | |
554 | math_Vector UV2(1,2), UV1(1,2); | |
555 | UV2(1)=myU; | |
556 | UV2(2)=myV-aDeltaV; | |
557 | UV1(1)=myU; | |
558 | UV1(2)=myV-2*aDeltaV; | |
559 | if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) | |
560 | { | |
561 | #ifdef DEB | |
562 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
563 | cout << "There are many points close to singularity points " | |
564 | "and which have zero-derivative." << endl; | |
565 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
566 | #endif | |
567 | return Standard_False; | |
568 | } | |
569 | ||
570 | F1 = FF1(1); | |
571 | F2 = FF2(1); | |
572 | ||
573 | Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV); | |
574 | }//else of if(myV-myVinfium < 2*aDeltaV) | |
575 | /////////////////////////// Search of DF1_v derivative (end) /////////////////// | |
576 | ||
577 | //Return to previous values | |
578 | myU = myU_old; | |
579 | myV = myV_old; | |
580 | myP1 = myP1_old, myP2 = myP2_old; | |
581 | myDu = myDu_old, myDv = myDv_old; | |
582 | }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) | |
583 | else | |
584 | { | |
585 | const Standard_Real Ndu = myDu.Magnitude(); | |
586 | Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu)); | |
587 | Df(1,2) = myDv.Dot(myDu)/Ndu; | |
588 | }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1)) | |
7fd59977 | 589 | |
32ca7a51 | 590 | if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) |
591 | { | |
592 | //Derivative is approximated by three points | |
593 | ||
594 | math_Vector FF1(1,2), FF2(1,2), FF3(1,2); | |
595 | Standard_Real F1, F2, F3; | |
596 | ||
597 | /////////////////////////// Search of DF2_v derivative (begin) /////////////////// | |
598 | if(myV-myVinfium < 2*aDeltaV) | |
599 | { | |
600 | F1=F(2); | |
601 | math_Vector UV2(1,2), UV3(1,2); | |
602 | UV2(1)=myU; | |
603 | UV2(2)=myV+aDeltaV; | |
604 | UV3(1)=myU; | |
605 | UV3(2)=myV+2*aDeltaV; | |
606 | ||
607 | if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) | |
608 | { | |
609 | #ifdef DEB | |
610 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
611 | cout << "There are many points close to singularity points " | |
612 | "and which have zero-derivative." << endl; | |
613 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
614 | #endif | |
615 | return Standard_False; | |
616 | } | |
617 | ||
618 | F2 = FF2(2); | |
619 | F3 = FF3(2); | |
620 | ||
621 | Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV); | |
622 | ||
623 | }//if(myV-myVinfium < 2*aDeltaV) | |
624 | else | |
625 | { | |
626 | F3 = F(2); | |
627 | math_Vector UV2(1,2), UV1(1,2); | |
628 | UV2(1)=myU; | |
629 | UV2(2)=myV-aDeltaV; | |
630 | UV1(1)=myU; | |
631 | UV1(2)=myV-2*aDeltaV; | |
632 | ||
633 | if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) | |
634 | { | |
635 | #ifdef DEB | |
636 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
637 | cout << "There are many points close to singularity points " | |
638 | "and which have zero-derivative." << endl; | |
639 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
640 | #endif | |
641 | return Standard_False; | |
642 | } | |
643 | ||
644 | F1 = FF1(2); | |
645 | F2 = FF2(2); | |
646 | ||
647 | Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV); | |
648 | }//else of if(myV-myVinfium < 2*aDeltaV) | |
649 | /////////////////////////// Search of DF2_v derivative (end) /////////////////// | |
650 | ||
651 | //Return to previous values | |
652 | myU = myU_old; | |
653 | myV = myV_old; | |
654 | ||
655 | /////////////////////////// Search of DF2_u derivative (begin) /////////////////// | |
656 | if(myU-myUinfium < 2*aDeltaU) | |
657 | { | |
658 | F1=F(2); | |
659 | math_Vector UV2(1,2), UV3(1,2); | |
660 | UV2(1)=myU+aDeltaU; | |
661 | UV2(2)=myV; | |
662 | UV3(1)=myU+2*aDeltaU; | |
663 | UV3(2)=myV; | |
664 | if(!((Value(UV2,FF2)) && (Value(UV3,FF3)))) | |
665 | { | |
666 | #ifdef DEB | |
667 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
668 | cout << "There are many points close to singularity points " | |
669 | "and which have zero-derivative." << endl; | |
670 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
671 | #endif | |
672 | return Standard_False; | |
673 | } | |
674 | ||
675 | F2 = FF2(2); | |
676 | F3 = FF3(2); | |
677 | ||
678 | Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU); | |
679 | }//if(myU-myUinfium < 2*aDelta) | |
680 | else | |
681 | { | |
682 | F3 = F(2); | |
683 | math_Vector UV2(1,2), UV1(1,2); | |
684 | UV2(1)=myU-aDeltaU; | |
685 | UV2(2)=myV; | |
686 | UV1(1)=myU-2*aDeltaU; | |
687 | UV1(2)=myV; | |
688 | ||
689 | if(!((Value(UV2,FF2)) && (Value(UV1,FF1)))) | |
690 | { | |
691 | #ifdef DEB | |
692 | cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl; | |
693 | cout << "There are many points close to singularity points " | |
694 | "and which have zero-derivative." << endl; | |
695 | cout << "Try to decrease aDelta variable's value. ---" << endl; | |
696 | #endif | |
697 | return Standard_False; | |
698 | } | |
699 | ||
700 | F1 = FF1(2); | |
701 | F2 = FF2(2); | |
702 | ||
703 | Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU); | |
704 | }//else of if(myU-myUinfium < 2*aDeltaU) | |
705 | /////////////////////////// Search of DF2_u derivative (end) /////////////////// | |
706 | ||
707 | //Return to previous values | |
708 | myU = myU_old; | |
709 | myV = myV_old; | |
710 | myP1 = myP1_old; | |
711 | myP2 = myP2_old; | |
712 | myDu = myDu_old; | |
713 | myDv = myDv_old; | |
714 | }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) | |
715 | else | |
716 | { | |
717 | Standard_Real Ndv = myDv.Magnitude(); | |
718 | Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv)); | |
719 | Df(2,1) = -myDu.Dot(myDv)/Ndv; | |
720 | }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2)) | |
7fd59977 | 721 | |
7fd59977 | 722 | return Standard_True; |
723 | ||
32ca7a51 | 724 | }//end of function |
7fd59977 | 725 | //============================================================================= |
726 | ||
727 | Standard_Integer Extrema_FuncExtCC::GetStateNumber () | |
728 | { | |
0d167958 RL |
729 | Vec Du (myDu), Dv (myDv); |
730 | Vec P1P2 (myP1, myP2); | |
731 | ||
7fd59977 | 732 | Standard_Real mod = Du.Magnitude(); |
32ca7a51 | 733 | if(mod > myTolC1) { |
7fd59977 | 734 | Du /= mod; |
735 | } | |
736 | mod = Dv.Magnitude(); | |
32ca7a51 | 737 | if(mod > myTolC2) { |
7fd59977 | 738 | Dv /= mod; |
739 | } | |
740 | ||
741 | if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol) { | |
742 | mySqDist.Append(myP1.SquareDistance(myP2)); | |
743 | myPoints.Append(POnC(myU, myP1)); | |
744 | myPoints.Append(POnC(myV, myP2)); | |
745 | } | |
746 | return 0; | |
747 | } | |
748 | //============================================================================= | |
749 | ||
750 | void Extrema_FuncExtCC::Points (const Standard_Integer N, | |
751 | POnC& P1, | |
752 | POnC& P2) const | |
753 | { | |
754 | P1 = myPoints.Value(2*N-1); | |
755 | P2 = myPoints.Value(2*N); | |
756 | } | |
757 | //============================================================================= | |
758 | ||
32ca7a51 | 759 | void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound, |
760 | const math_Vector& theSupBound) | |
761 | { | |
762 | myUinfium = theInfBound(1); | |
763 | myUsupremum = theSupBound(1); | |
764 | myVinfium = theInfBound(2); | |
765 | myVsupremum = theSupBound(2); | |
766 | } | |
767 | //============================================================================= |