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