1 // Created on: 1999-09-16
2 // Created by: Edward AGAPOV
3 // Copyright (c) 1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 #include <Standard_NotImplemented.hxx>
18 #include <Standard_OutOfRange.hxx>
19 #include <StdFail_NotDone.hxx>
20 #include <Adaptor3d_HCurve.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>
34 #include <math_FunctionSetRoot.hxx>
35 #include <math_Vector.hxx>
36 #include <GeomAdaptor_HSurfaceOfLinearExtrusion.hxx>
39 IMPLEMENT_STANDARD_RTTIEXT(Extrema_ExtPExtS,Standard_Transient)
41 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C);
43 static void PerformExtPElC (Extrema_ExtPElC& E,
45 const Handle(Adaptor3d_HCurve)& C,
46 const Standard_Real Tol);
48 static Standard_Boolean
49 IsCaseAnalyticallyComputable (const GeomAbs_CurveType& theType,
50 const gp_Ax2& theCurvePos,
51 const gp_Dir& theSurfaceDirection) ;
53 static gp_Pnt GetValue(const Standard_Real U,
54 const Handle(Adaptor3d_HCurve)& C);
55 //=======================================================================
57 //purpose : Returns the projection of a point <Point> on a plane
58 // <ThePlane> along a direction <TheDir>.
59 //=======================================================================
60 static gp_Pnt ProjectPnt(const gp_Ax2& ThePlane,
64 gp_Vec PO(Point,ThePlane.Location());
65 Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
66 Alpha /= TheDir * ThePlane.Direction();
68 P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
72 //=======================================================================
73 //function : IsOriginalPnt
75 //=======================================================================
77 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
78 const Extrema_POnSurf* Points,
79 const Standard_Integer NbPoints)
81 for (Standard_Integer i=1; i<=NbPoints; i++) {
82 if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
83 return Standard_False;
89 //=======================================================================
90 //function : MakePreciser
92 //=======================================================================
94 void Extrema_ExtPExtS::MakePreciser (Standard_Real& U,
96 const Standard_Boolean isMin,
97 const gp_Ax2& OrtogSection) const
101 } else if (U<myuinf) {
105 Standard_Real step = (myusup - myuinf) / 30, D2e, D2next ,D2prev;
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;
115 notFound = (D2e > D2prev || D2e > D2next);
117 notFound = (D2e < D2prev || D2e < D2next);
119 if (notFound && (D2e < D2next && isMin)) {
136 Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
137 D2next = P.SquareDistance(Pnext);
139 notFound = D2e > D2next;
141 notFound = D2e < D2next;
145 //=============================================================================
147 Extrema_ExtPExtS::Extrema_ExtPExtS()
154 myIsAnalyticallyComputable(Standard_False),
155 myDone(Standard_False),
156 myNbExt(Standard_False)
160 //=============================================================================
162 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
163 const Handle(GeomAdaptor_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)
177 myIsAnalyticallyComputable(Standard_False),
178 myDone(Standard_False),
179 myNbExt(Standard_False)
191 //=============================================================================
193 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt& theP,
194 const Handle(GeomAdaptor_HSurfaceOfLinearExtrusion)& theS,
195 const Standard_Real theTolU,
196 const Standard_Real theTolV)
197 : myuinf(theS->FirstUParameter()),
198 myusup(theS->LastUParameter()),
200 myvinf(theS->FirstVParameter()),
201 myvsup(theS->LastVParameter()),
204 myIsAnalyticallyComputable(Standard_False),
205 myDone(Standard_False),
206 myNbExt(Standard_False)
209 theS->FirstUParameter(),
210 theS->LastUParameter(),
211 theS->FirstVParameter(),
212 theS->LastVParameter(),
219 //=======================================================================
220 //function : Initialize
222 //=======================================================================
224 void Extrema_ExtPExtS::Initialize (const Handle(GeomAdaptor_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)
240 Handle(Adaptor3d_HCurve) anACurve = theS->BasisCurve();
242 myF.Initialize (theS->ChangeSurface());
245 myPosition = GetPosition(myC);
246 myDirection = theS->Direction();
247 myIsAnalyticallyComputable = //Standard_False;
248 IsCaseAnalyticallyComputable (myC->GetType(), myPosition, myDirection);
250 if (!myIsAnalyticallyComputable)
252 myExtPS.Initialize (theS->ChangeSurface(),
265 //=======================================================================
268 //=======================================================================
270 void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
272 const Standard_Integer NbExtMax = 4; //dimension of arrays
273 //myPoint[] and mySqDist[]
274 //For "analytical" case
275 myDone = Standard_False;
278 if (!myIsAnalyticallyComputable) {
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___
287 gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
288 Extrema_ExtPElC anExt;
289 PerformExtPElC(anExt, Pp, myC, mytolu);
290 if (!anExt.IsDone()) return;
292 gp_Ax2 anOrtogSection (P, myDirection);
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;
305 for (i=1; i<=aNbExt; i++) {
306 Extrema_POnCurv POC=anExt.Point(i);
308 //// modified by jgv, 23.12.2008 for OCC17194 ////
309 if (myC->IsPeriodic())
311 Standard_Real U2 = U;
312 ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
314 //////////////////////////////////////////////////
315 gp_Pnt E = POC.Value();
316 Pe = ProjectPnt(anOrtogSection, myDirection, E);
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);
326 if(myNbExt == NbExtMax)
330 // modified by NIZHNY-MKK Thu Sep 18 14:46:18 2003.END
334 isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
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, Tol);
342 aFSR.Perform(myF, UV, UVinf, UVsup);
343 // for (Standard_Integer k=1 ; k <= myF.NbExt();
345 for ( k=1 ; k <= myF.NbExt(); k++) {
346 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
347 // modified by NIZHNY-MKK Thu Sep 18 14:46:41 2003.BEGIN
348 // myPoint[++myNbExt] = myF.Point(k);
349 // myValue[myNbExt] = myF.Value(k);
350 myPoint[myNbExt] = myF.Point(k);
351 mySqDist[myNbExt] = myF.SquareDistance(k);
353 if(myNbExt == NbExtMax)
357 // modified by NIZHNY-MKK Thu Sep 18 14:46:43 2003.END
360 if(myNbExt == NbExtMax)
364 // try symmetric point
365 myF.SetPoint(P); //To clear previous solutions
367 MakePreciser(U, P, isMin, anOrtogSection);
368 E = GetValue(U, myC);
369 Pe = ProjectPnt (anOrtogSection, myDirection, E),
370 V = gp_Vec(E,Pe) * gp_Vec(myDirection);
371 UV(1) = U; UV(2) = V;
373 aFSR.Perform (myF,UV,UVinf,UVsup);
375 for (k=1 ; k <= myF.NbExt(); k++) {
376 if(myF.SquareDistance(k) > Precision::Confusion()*Precision::Confusion())
378 //Additional checking solution: FSR sometimes is wrong
379 //when starting point is far from solution.
380 Standard_Real dist = Sqrt(myF.SquareDistance(k));
381 math_Vector Vals(1, 2);
382 const Extrema_POnSurf& PonS=myF.Point(k);
384 PonS.Parameter(u, v);
389 myS->D1(u, v, Pe, du, dv);
390 Standard_Real mdu = du.Magnitude();
391 Standard_Real mdv = dv.Magnitude();
394 if(mdu > Precision::PConfusion())
396 if(u / dist / mdu > Precision::PConfusion())
401 if(mdv > Precision::PConfusion())
403 if(v / dist / mdv > Precision::PConfusion())
410 if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
411 // modified by NIZHNY-MKK Thu Sep 18 14:46:59 2003.BEGIN
412 // myPoint[++myNbExt] = myF.Point(k);
413 // myValue[myNbExt] = myF.Value(k);
414 myPoint[myNbExt] = myF.Point(k);
415 mySqDist[myNbExt] = myF.SquareDistance(k);
417 if(myNbExt == NbExtMax)
421 // modified by NIZHNY-MKK Thu Sep 18 14:47:04 2003.END
424 if(myNbExt == NbExtMax)
430 myDone = Standard_True;
434 //=============================================================================
436 Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
437 //=============================================================================
439 Standard_Integer Extrema_ExtPExtS::NbExt () const
441 if (!IsDone()) { throw StdFail_NotDone(); }
442 if (myIsAnalyticallyComputable)
445 return myExtPS.NbExt();
447 //=============================================================================
449 Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
451 if (!IsDone()) { throw StdFail_NotDone(); }
452 if ((N < 1) || (N > myNbExt)) { throw Standard_OutOfRange(); }
453 if (myIsAnalyticallyComputable)
454 // modified by NIZHNY-MKK Thu Sep 18 14:48:39 2003.BEGIN
455 // return myValue[N];
456 return mySqDist[N-1];
457 // modified by NIZHNY-MKK Thu Sep 18 14:48:42 2003.END
459 return myExtPS.SquareDistance(N);
461 //=============================================================================
463 const Extrema_POnSurf& Extrema_ExtPExtS::Point (const Standard_Integer N) const
465 if (!IsDone()) { throw StdFail_NotDone(); }
466 if ((N < 1) || (N > myNbExt)) { throw Standard_OutOfRange(); }
467 if (myIsAnalyticallyComputable) {
468 // modified by NIZHNY-MKK Thu Sep 18 14:47:40 2003.BEGIN
469 // return myPoint[N];
472 // modified by NIZHNY-MKK Thu Sep 18 14:47:43 2003.END
474 return myExtPS.Point(N);
476 //=============================================================================
479 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C)
481 switch (C->GetType()) {
483 gp_Lin L = C->Line();
484 gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
485 //:abv 30.05.02: OCC - use constructor instead of Set...s() to avoid exception
486 gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
487 // Pos.SetAxis(Pln.XAxis());
488 // Pos.SetXDirection(Pln.YAxis().Direction());
489 // Pos.SetYDirection(Pln.Position().Direction());
493 return C->Circle().Position();
494 case GeomAbs_Ellipse:
495 return C->Ellipse().Position();
496 case GeomAbs_Hyperbola:
497 return C->Hyperbola().Position();
498 case GeomAbs_Parabola:
499 return C->Parabola().Position();
504 //=============================================================================
506 static void PerformExtPElC (Extrema_ExtPElC& E,
508 const Handle(Adaptor3d_HCurve)& C,
509 const Standard_Real Tol)
511 switch (C->GetType()) {
512 case GeomAbs_Hyperbola:
513 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
516 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
519 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
521 case GeomAbs_Ellipse:
522 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
524 case GeomAbs_Parabola:
525 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
532 //=======================================================================
533 //function : IsCaseAnalyticallyComputable
535 //=======================================================================
537 static Standard_Boolean IsCaseAnalyticallyComputable
538 (const GeomAbs_CurveType& theType,
539 const gp_Ax2& theCurvePos,
540 const gp_Dir& theSurfaceDirection)
546 case GeomAbs_Ellipse:
547 case GeomAbs_Hyperbola:
548 case GeomAbs_Parabola:
551 return Standard_False;
553 // check if it is a plane
554 if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
555 return Standard_False;
557 return Standard_True;
559 //=======================================================================
560 //function : GetValue
562 //=======================================================================
564 static gp_Pnt GetValue(const Standard_Real U,
565 const Handle(Adaptor3d_HCurve)& C)
567 switch (C->GetType()) {
569 return ElCLib::Value(U, C->Line());
571 return ElCLib::Value(U, C->Circle());
572 case GeomAbs_Ellipse:
573 return ElCLib::Value(U, C->Ellipse());
574 case GeomAbs_Hyperbola:
575 return ElCLib::Value(U, C->Hyperbola());
576 case GeomAbs_Parabola:
577 return ElCLib::Value(U, C->Parabola());
582 //=======================================================================
585 //=======================================================================
587 //static Standard_Real GetU(const gp_Vec& vec,
589 // const Handle(Adaptor3d_HCurve)& C)
591 // switch (C->GetType()) {
592 // case GeomAbs_Line:
593 // return ElCLib::Parameter(C->Line().Translated(vec), P);
594 // case GeomAbs_Circle:
595 // return ElCLib::Parameter(C->Circle().Translated(vec), P);
596 // case GeomAbs_Ellipse:
597 // return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
598 // case GeomAbs_Hyperbola:
599 // return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
600 // case GeomAbs_Parabola:
601 // return ElCLib::Parameter(C->Parabola().Translated(vec), P);