projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0025621: CAST analysis - Avoid constructors not supplying an initial value for all...
[occt.git] / src / Extrema / Extrema_ExtPExtS.cxx
CommitLineData
b311480e 1// Created on: 1999-09-16
2// Created by: Edward AGAPOV
3// Copyright (c) 1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17#include <Standard_NotImplemented.hxx>
18#include <Standard_OutOfRange.hxx>
19#include <StdFail_NotDone.hxx>
20#include <Adaptor3d_HCurve.hxx>
21#include <ElCLib.hxx>
22#include <Extrema_ExtPElC.hxx>
23#include <Extrema_ExtPExtS.hxx>
24#include <Extrema_POnCurv.hxx>
25#include <Extrema_POnSurf.hxx>
26#include <Precision.hxx>
27#include <gp.hxx>
28#include <gp_Ax2.hxx>
29#include <gp_Dir.hxx>
30#include <gp_Lin.hxx>
31#include <gp_Pln.hxx>
32#include <gp_Pnt.hxx>
33#include <gp_Vec.hxx>
34#include <math_FunctionSetRoot.hxx>
35#include <math_Vector.hxx>
0734c53d 36#include <Adaptor3d_HSurfaceOfLinearExtrusion.hxx>
37
38IMPLEMENT_STANDARD_HANDLE (Extrema_ExtPExtS, Standard_Transient)
39IMPLEMENT_STANDARD_RTTIEXT(Extrema_ExtPExtS, Standard_Transient)
7fd59977 40
41static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C);
42
43static void PerformExtPElC (Extrema_ExtPElC& E,
44 const gp_Pnt& P,
45 const Handle(Adaptor3d_HCurve)& C,
46 const Standard_Real Tol);
47
48static Standard_Boolean
49 IsCaseAnalyticallyComputable (const GeomAbs_CurveType& theType,
50 const gp_Ax2& theCurvePos,
51 const gp_Dir& theSurfaceDirection) ;
52
53static gp_Pnt GetValue(const Standard_Real U,
54 const Handle(Adaptor3d_HCurve)& C);
55//=======================================================================
56//function : Project
57//purpose : Returns the projection of a point <Point> on a plane
58// <ThePlane> along a direction <TheDir>.
59//=======================================================================
60static gp_Pnt ProjectPnt(const gp_Ax2& ThePlane,
61 const gp_Dir& TheDir,
62 const gp_Pnt& Point)
63{
64 gp_Vec PO(Point,ThePlane.Location());
65 Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
66 Alpha /= TheDir * ThePlane.Direction();
67 gp_Pnt P;
68 P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
69 return P;
70}
71
72//=======================================================================
73//function : IsOriginalPnt
74//purpose :
75//=======================================================================
76
77static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
78 const Extrema_POnSurf* Points,
79 const Standard_Integer NbPoints)
80{
81 for (Standard_Integer i=1; i<=NbPoints; i++) {
82 if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
83 return Standard_False;
84 }
85 }
86 return Standard_True;
87}
88
89//=======================================================================
90//function : MakePreciser
91//purpose :
92//=======================================================================
93
94void Extrema_ExtPExtS::MakePreciser (Standard_Real& U,
95 const gp_Pnt& P,
96 const Standard_Boolean isMin,
97 const gp_Ax2& OrtogSection) const
98{
99 if (U>myusup) {
100 U = myusup;
101 } else if (U<myuinf) {
102 U = myuinf;
103 } else {
104
105 Standard_Real step = (myusup - myuinf) / 30, D2e, D2next ,D2prev;
106 gp_Pnt
107 Pe = ProjectPnt (OrtogSection, myDirection, GetValue(U,myC)),
108 Pprev = ProjectPnt (OrtogSection, myDirection, GetValue(U-step, myC)),
109 Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
110 D2e = P.SquareDistance(Pe),
111 D2next = P.SquareDistance(Pnext),
112 D2prev = P.SquareDistance(Pprev);
113 Standard_Boolean notFound;
114 if (isMin)
115 notFound = (D2e > D2prev || D2e > D2next);
116 else
117 notFound = (D2e < D2prev || D2e < D2next);
118
119 if (notFound && (D2e < D2next && isMin)) {
120 step = -step;
121 D2next = D2prev;
122 Pnext = Pprev;
123 }
124 while (notFound) {
125 U = U + step;
126 if (U > myusup) {
aabe3a17 127 U = myusup;
128 break;
7fd59977 129 }
130 if (U < myuinf) {
aabe3a17 131 U = myuinf;
132 break;
7fd59977 133 }
134 D2e = D2next;
135 Pe = Pnext;
136 Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
137 D2next = P.SquareDistance(Pnext);
138 if (isMin)
aabe3a17 139 notFound = D2e > D2next;
7fd59977 140 else
aabe3a17 141 notFound = D2e < D2next;
7fd59977 142 }
143 }
144}
145//=============================================================================
146
cbff1e55 147Extrema_ExtPExtS::Extrema_ExtPExtS()
148: myuinf(0.0),
149 myusup(0.0),
150 mytolu(0.0),
151 myvinf(0.0),
152 myvsup(0.0),
153 mytolv(0.0),
154 myIsAnalyticallyComputable(Standard_False),
155 myDone(Standard_False),
156 myNbExt(Standard_False)
7fd59977 157{
7fd59977 158}
159
160//=============================================================================
161
0734c53d 162Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
163 const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)& theS,
164 const Standard_Real theUmin,
165 const Standard_Real theUsup,
166 const Standard_Real theVmin,
167 const Standard_Real theVsup,
168 const Standard_Real theTolU,
169 const Standard_Real theTolV)
cbff1e55 170: myuinf(theUmin),
171 myusup(theUsup),
172 mytolu(theTolU),
173 myvinf(theVmin),
174 myvsup(theVsup),
175 mytolv(theTolV),
176 myS (theS),
177 myIsAnalyticallyComputable(Standard_False),
178 myDone(Standard_False),
179 myNbExt(Standard_False)
7fd59977 180{
0734c53d 181 Initialize (theS,
182 theUmin,
183 theUsup,
184 theVmin,
185 theVsup,
186 theTolU,
187 theTolV);
188
189 Perform (theP);
7fd59977 190}
191//=============================================================================
192
0734c53d 193Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
194 const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)& theS,
195 const Standard_Real theTolU,
196 const Standard_Real theTolV)
cbff1e55 197: myuinf(theS->FirstUParameter()),
198 myusup(theS->LastUParameter()),
199 mytolu(theTolU),
200 myvinf(theS->FirstVParameter()),
201 myvsup(theS->LastVParameter()),
202 mytolv(theTolV),
203 myS (theS),
204 myIsAnalyticallyComputable(Standard_False),
205 myDone(Standard_False),
206 myNbExt(Standard_False)
7fd59977 207{
0734c53d 208 Initialize (theS,
209 theS->FirstUParameter(),
210 theS->LastUParameter(),
211 theS->FirstVParameter(),
212 theS->LastVParameter(),
213 theTolU,
214 theTolV);
215
216 Perform (theP);
7fd59977 217}
cbff1e55 218
7fd59977 219//=======================================================================
220//function : Initialize
221//purpose :
222//=======================================================================
223
0734c53d 224void Extrema_ExtPExtS::Initialize (const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)& theS,
225 const Standard_Real theUinf,
226 const Standard_Real theUsup,
227 const Standard_Real theVinf,
228 const Standard_Real theVsup,
229 const Standard_Real theTolU,
230 const Standard_Real theTolV)
7fd59977 231{
0734c53d 232 myuinf = theUinf;
233 myusup = theUsup;
234 mytolu = theTolU;
7fd59977 235
0734c53d 236 myvinf = theVinf;
237 myvsup = theVsup;
238 mytolv = theTolV;
7fd59977 239
0734c53d 240 Handle(Adaptor3d_HCurve) anACurve = theS->BasisCurve();
7fd59977 241
0734c53d 242 myF.Initialize (theS->ChangeSurface());
7fd59977 243 myC = anACurve;
0734c53d 244 myS = theS;
7fd59977 245 myPosition = GetPosition(myC);
0734c53d 246 myDirection = theS->Direction();
7fd59977 247 myIsAnalyticallyComputable = //Standard_False;
0734c53d 248 IsCaseAnalyticallyComputable (myC->GetType(), myPosition, myDirection);
7fd59977 249
250 if (!myIsAnalyticallyComputable)
0734c53d 251 {
252 myExtPS.Initialize (theS->ChangeSurface(),
253 32,
254 32,
255 theUinf,
256 theUsup,
257 theVinf,
258 theVsup,
259 theTolU,
260 theTolV);
261 }
7fd59977 262}
263
264
265//=======================================================================
266//function : Perform
267//purpose :
268//=======================================================================
269
270void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
271{
aabe3a17 272 const Standard_Integer NbExtMax = 4; //dimension of arrays
273 //myPoint[] and mySqDist[]
274 //For "analytical" case
7fd59977 275 myDone = Standard_False;
276 myNbExt = 0;
277
278 if (!myIsAnalyticallyComputable) {
279 myExtPS.Perform(P);
280 myDone = myExtPS.IsDone();
281// modified by NIZHNY-EAP Wed Nov 17 12:59:08 1999 ___BEGIN___
282 myNbExt = myExtPS.NbExt();
283// modified by NIZHNY-EAP Wed Nov 17 12:59:09 1999 ___END___
284 return;
285 }
286
287 gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
288 Extrema_ExtPElC anExt;
289 PerformExtPElC(anExt, Pp, myC, mytolu);
290 if (!anExt.IsDone()) return;
291
292 gp_Ax2 anOrtogSection (P, myDirection);
293 Standard_Real U,V;
294 Standard_Boolean
295 isMin,
296 isSimpleCase =
297 myDirection.IsParallel(myPosition.Direction(),Precision::Angular());
298 Standard_Integer i, aNbExt = anExt.NbExt();
299 math_Vector UV(1,2), Tol(1,2), UVinf(1,2), UVsup(1,2);
300 Tol(1) = mytolu; Tol(2) = mytolv;
301 UVinf(1) = myuinf; UVinf(2) = myvinf;
302 UVsup(1) = myusup; UVsup(2) = myvsup;
303
304
305 for (i=1; i<=aNbExt; i++) {
306 Extrema_POnCurv POC=anExt.Point(i);
307 U = POC.Parameter();
308 //// modified by jgv, 23.12.2008 for OCC17194 ////
309 if (myC->IsPeriodic())
aabe3a17 310 {
311 Standard_Real U2 = U;
312 ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
313 }
7fd59977 314 //////////////////////////////////////////////////
315 gp_Pnt E = POC.Value();
316 Pe = ProjectPnt(anOrtogSection, myDirection, E);
317
318 if (isSimpleCase) {
319 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
320 // modified by NIZHNY-MKK Thu Sep 18 14:46:14 2003.BEGIN
321 // myPoint[++myNbExt] = Extrema_POnSurf(U, V, Pe);
322 // myValue[myNbExt] = anExt.Value(i);
323 myPoint[myNbExt] = Extrema_POnSurf(U, V, Pe);
324 mySqDist[myNbExt] = anExt.SquareDistance(i);
325 myNbExt++;
aabe3a17 326 if(myNbExt == NbExtMax)
327 {
328 break;
329 }
7fd59977 330 // modified by NIZHNY-MKK Thu Sep 18 14:46:18 2003.END
331 }
332 else {
333 myF.SetPoint(P);
334 isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
335
336 MakePreciser(U, P, isMin, anOrtogSection);
337 E = GetValue(U, myC);
338 Pe = ProjectPnt (anOrtogSection, myDirection, E),
339 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
340 UV(1) = U; UV(2) = V;
341 math_FunctionSetRoot aFSR (myF,UV,Tol,UVinf,UVsup);
342// for (Standard_Integer k=1 ; k <= myF.NbExt();
343 Standard_Integer k;
344 for ( k=1 ; k <= myF.NbExt(); k++) {
aabe3a17 345 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
346 // modified by NIZHNY-MKK Thu Sep 18 14:46:41 2003.BEGIN
347 // myPoint[++myNbExt] = myF.Point(k);
348 // myValue[myNbExt] = myF.Value(k);
349 myPoint[myNbExt] = myF.Point(k);
350 mySqDist[myNbExt] = myF.SquareDistance(k);
351 myNbExt++;
352 if(myNbExt == NbExtMax)
353 {
354 break;
355 }
356 // modified by NIZHNY-MKK Thu Sep 18 14:46:43 2003.END
357 }
358 }
359 if(myNbExt == NbExtMax)
360 {
361 break;
7fd59977 362 }
363 // try symmetric point
aabe3a17 364 myF.SetPoint(P); //To clear previous solutions
7fd59977 365 U *= -1;
366 MakePreciser(U, P, isMin, anOrtogSection);
367 E = GetValue(U, myC);
368 Pe = ProjectPnt (anOrtogSection, myDirection, E),
369 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
370 UV(1) = U; UV(2) = V;
371
372 aFSR.Perform (myF,UV,UVinf,UVsup);
373
374 for (k=1 ; k <= myF.NbExt(); k++) {
aabe3a17 375 if(myF.SquareDistance(k) > Precision::Confusion()*Precision::Confusion())
376 {
377 //Additional checking solution: FSR sometimes is wrong
378 //when starting point is far from solution.
379 Standard_Real dist = Sqrt(myF.SquareDistance(k));
380 math_Vector Vals(1, 2);
381 const Extrema_POnSurf& PonS=myF.Point(k);
382 Standard_Real u, v;
383 PonS.Parameter(u, v);
384 UV(1) = u;
385 UV(2) = v;
386 myF.Value(UV, Vals);
387 gp_Vec du, dv;
388 myS->D1(u, v, Pe, du, dv);
389 Standard_Real mdu = du.Magnitude();
390 Standard_Real mdv = dv.Magnitude();
391 u = Abs(Vals(1));
392 v = Abs(Vals(2));
393 if(mdu > Precision::PConfusion())
394 {
395 if(u / dist / mdu > Precision::PConfusion())
396 {
397 continue;
398 }
399 }
400 if(mdv > Precision::PConfusion())
401 {
402 if(v / dist / mdv > Precision::PConfusion())
403 {
404 continue;
405 }
406 }
407
408 }
409 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
410 // modified by NIZHNY-MKK Thu Sep 18 14:46:59 2003.BEGIN
411 // myPoint[++myNbExt] = myF.Point(k);
412 // myValue[myNbExt] = myF.Value(k);
413 myPoint[myNbExt] = myF.Point(k);
414 mySqDist[myNbExt] = myF.SquareDistance(k);
415 myNbExt++;
416 if(myNbExt == NbExtMax)
417 {
418 break;
419 }
420 // modified by NIZHNY-MKK Thu Sep 18 14:47:04 2003.END
421 }
422 }
423 if(myNbExt == NbExtMax)
424 {
425 break;
7fd59977 426 }
427 }
428 }
429 myDone = Standard_True;
430 return;
431}
432
433//=============================================================================
434
435Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
436//=============================================================================
437
438Standard_Integer Extrema_ExtPExtS::NbExt () const
439{
440 if (!IsDone()) { StdFail_NotDone::Raise(); }
441 if (myIsAnalyticallyComputable)
442 return myNbExt;
443 else
444 return myExtPS.NbExt();
445}
446//=============================================================================
447
448Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
449{
450 if (!IsDone()) { StdFail_NotDone::Raise(); }
451 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
452 if (myIsAnalyticallyComputable)
453 // modified by NIZHNY-MKK Thu Sep 18 14:48:39 2003.BEGIN
454 // return myValue[N];
455 return mySqDist[N-1];
456 // modified by NIZHNY-MKK Thu Sep 18 14:48:42 2003.END
457 else
458 return myExtPS.SquareDistance(N);
459}
460//=============================================================================
461
5d99f2c8 462const Extrema_POnSurf& Extrema_ExtPExtS::Point (const Standard_Integer N) const
7fd59977 463{
464 if (!IsDone()) { StdFail_NotDone::Raise(); }
465 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
466 if (myIsAnalyticallyComputable) {
467 // modified by NIZHNY-MKK Thu Sep 18 14:47:40 2003.BEGIN
468 // return myPoint[N];
469 return myPoint[N-1];
470 }
471 // modified by NIZHNY-MKK Thu Sep 18 14:47:43 2003.END
472 else
473 return myExtPS.Point(N);
474}
475//=============================================================================
476
477
478static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C)
479{
480 switch (C->GetType()) {
481 case GeomAbs_Line: {
482 gp_Lin L = C->Line();
483 gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
484 //:abv 30.05.02: OCC - use constructor instead of Set...s() to avoid exception
485 gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
486// Pos.SetAxis(Pln.XAxis());
487// Pos.SetXDirection(Pln.YAxis().Direction());
488// Pos.SetYDirection(Pln.Position().Direction());
489 return Pos;
490 }
491 case GeomAbs_Circle:
492 return C->Circle().Position();
493 case GeomAbs_Ellipse:
494 return C->Ellipse().Position();
495 case GeomAbs_Hyperbola:
496 return C->Hyperbola().Position();
497 case GeomAbs_Parabola:
498 return C->Parabola().Position();
499 default:
500 return gp_Ax2 ();
501 }
502}
503//=============================================================================
504
505static void PerformExtPElC (Extrema_ExtPElC& E,
506 const gp_Pnt& P,
507 const Handle(Adaptor3d_HCurve)& C,
508 const Standard_Real Tol)
509{
510 switch (C->GetType()) {
511 case GeomAbs_Hyperbola:
512 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
513 return;
514 case GeomAbs_Line:
515 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
516 return;
517 case GeomAbs_Circle:
c6541a0c 518 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
7fd59977 519 return;
520 case GeomAbs_Ellipse:
c6541a0c 521 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
7fd59977 522 return;
523 case GeomAbs_Parabola:
524 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
525 return;
7fd59977 526 default:
527 return;
7fd59977 528 }
529}
530
531//=======================================================================
532//function : IsCaseAnalyticallyComputable
533//purpose :
534//=======================================================================
535
536static Standard_Boolean IsCaseAnalyticallyComputable
537 (const GeomAbs_CurveType& theType,
538 const gp_Ax2& theCurvePos,
539 const gp_Dir& theSurfaceDirection)
540{
541 // check type
542 switch (theType) {
543 case GeomAbs_Line:
544 case GeomAbs_Circle:
545 case GeomAbs_Ellipse:
546 case GeomAbs_Hyperbola:
547 case GeomAbs_Parabola:
548 break;
549 default:
550 return Standard_False;
551 }
552 // check if it is a plane
553 if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
554 return Standard_False;
555 else
556 return Standard_True;
557}
558//=======================================================================
559//function : GetValue
560//purpose :
561//=======================================================================
562
563static gp_Pnt GetValue(const Standard_Real U,
564 const Handle(Adaptor3d_HCurve)& C)
565{
566 switch (C->GetType()) {
567 case GeomAbs_Line:
568 return ElCLib::Value(U, C->Line());
569 case GeomAbs_Circle:
570 return ElCLib::Value(U, C->Circle());
571 case GeomAbs_Ellipse:
572 return ElCLib::Value(U, C->Ellipse());
573 case GeomAbs_Hyperbola:
574 return ElCLib::Value(U, C->Hyperbola());
575 case GeomAbs_Parabola:
576 return ElCLib::Value(U, C->Parabola());
577 default:
578 return gp_Pnt();
579 }
580}
581//=======================================================================
582//function : GetU
583//purpose :
584//=======================================================================
0797d9d3 585//#ifdef OCCT_DEBUG
7fd59977 586//static Standard_Real GetU(const gp_Vec& vec,
587// const gp_Pnt& P,
588// const Handle(Adaptor3d_HCurve)& C)
589//{
590// switch (C->GetType()) {
591// case GeomAbs_Line:
592// return ElCLib::Parameter(C->Line().Translated(vec), P);
593// case GeomAbs_Circle:
594// return ElCLib::Parameter(C->Circle().Translated(vec), P);
595// case GeomAbs_Ellipse:
596// return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
597// case GeomAbs_Hyperbola:
598// return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
599// case GeomAbs_Parabola:
600// return ElCLib::Parameter(C->Parabola().Translated(vec), P);
601// default:
602// return 0;
603// }
604//}
605//#endif
Open CASCADE Technology repository