0023733: PCurve for edge on face creation failure
[occt.git] / src / Extrema / Extrema_ExtPRevS.cxx
1 // Created on: 1999-09-21
2 // Created by: Edward AGAPOV
3 // Copyright (c) 1999-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22
23 #include <Extrema_ExtPRevS.ixx>
24 #include <Adaptor3d_HCurve.hxx>
25 #include <Extrema_ExtPElC.hxx>
26 #include <Extrema_GenExtPS.hxx>
27 #include <Extrema_POnCurv.hxx>
28 #include <Extrema_POnSurf.hxx>
29 #include <gp_Ax1.hxx>
30 #include <gp_Ax2.hxx>
31 #include <gp_Lin.hxx>
32 #include <gp_Pln.hxx>
33 #include <gp_Pnt.hxx>
34 #include <gp_Trsf.hxx>
35 #include <gp_Vec.hxx>
36 #include <Precision.hxx>
37 #include <ElCLib.hxx>
38
39 static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
40 {
41   Handle(Adaptor3d_HCurve) C = S.BasisCurve();
42   
43   switch (C->GetType()) {
44     
45   case GeomAbs_Line: {
46     gp_Lin L = C->Line();
47     gp_Dir N = S.AxeOfRevolution().Direction();
48     if (N.IsParallel(L.Direction(), Precision::Angular())) {
49       gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
50       if (OO.Magnitude() <= gp::Resolution()) {
51         OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
52         if (N.IsParallel(OO, Precision::Angular()))
53           return gp_Ax2(); // Line and axe of revolution coinside
54       }
55       N ^= OO;
56     }
57     else {
58       N ^= L.Direction();
59     }
60     return gp_Ax2 (L.Location(), N, L.Direction());
61   }
62   case GeomAbs_Circle:
63     return C->Circle().Position();
64   case GeomAbs_Ellipse:
65     return C->Ellipse().Position();
66   case GeomAbs_Hyperbola:
67     return C->Hyperbola().Position();
68   case GeomAbs_Parabola:
69     return C->Parabola().Position();
70   default:
71     return gp_Ax2();
72   }
73 }
74
75 //=======================================================================
76 //function : HasSingularity
77 //purpose  : 
78 //=======================================================================
79
80 static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S) 
81 {
82
83   const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
84   gp_Dir N = S.AxeOfRevolution().Direction();
85   gp_Pnt P = S.AxeOfRevolution().Location();
86
87   gp_Lin L(P, N);
88
89   P = C->Value(C->FirstParameter());
90
91   if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
92
93   P = C->Value(C->LastParameter());
94
95   if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
96   
97   return Standard_False;
98 }
99
100 //=============================================================================
101
102 static void PerformExtPElC (Extrema_ExtPElC& E,
103                             const gp_Pnt& P,
104                             const Handle(Adaptor3d_HCurve)& C,
105                             const Standard_Real Tol)
106 {
107   switch (C->GetType()) {
108   case GeomAbs_Hyperbola:
109     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
110     return;
111   case GeomAbs_Line:
112     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
113     return;
114   case GeomAbs_Circle:
115     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
116     return;
117   case GeomAbs_Ellipse:
118     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
119     return;
120   case GeomAbs_Parabola:
121     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
122     return;
123   default:
124     return ;
125   }
126 }
127
128 //=======================================================================
129 //function : IsCaseAnalyticallyComputable
130 //purpose  : 
131 //=======================================================================
132
133 static Standard_Boolean IsCaseAnalyticallyComputable
134   (const GeomAbs_CurveType& theType,
135    const gp_Ax2& theCurvePos,
136    const gp_Ax1& AxeOfRevolution) 
137 {
138   // check type
139   switch (theType) {
140   case GeomAbs_Line:
141   case GeomAbs_Circle:
142   case GeomAbs_Ellipse:
143   case GeomAbs_Hyperbola:
144   case GeomAbs_Parabola:
145     break;
146   default:
147     return  Standard_False;
148   }
149 //  the axe of revolution must be in the plane of the curve.
150   gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
151   gp_Pnt p1 = AxeOfRevolution.Location();
152   Standard_Real dist = 100., dist2 = dist * dist;
153   Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
154   gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
155
156   if((pl.SquareDistance(p1) < aThreshold) && 
157      (pl.SquareDistance(p2) < aThreshold))
158     return Standard_True;
159   return Standard_False;
160   //   gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
161   //   if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
162   //     return Standard_True;
163   //   else
164   //     return Standard_False;
165 }
166
167 //=======================================================================
168 //function : IsOriginalPnt
169 //purpose  : 
170 //=======================================================================
171
172 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
173                                        const Extrema_POnSurf* Points,
174                                        const Standard_Integer NbPoints)
175 {
176   for (Standard_Integer i=1; i<=NbPoints; i++) {
177     if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
178       return Standard_False;
179     }
180   }
181   return Standard_True;
182 }
183 //=======================================================================
184 //function : IsExtremum
185 //purpose  : 
186 //=======================================================================
187
188 static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
189                                     const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
190                                     gp_Pnt& E,        Standard_Real& Dist2,
191                                     const Standard_Boolean IsVSup,
192                                     const Standard_Boolean IsMin)
193 {
194   E = S->Value(U,V);
195   Dist2 = P.SquareDistance(E);
196   if (IsMin) 
197     return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
198             Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
199             Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
200   else
201     return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
202             Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
203             Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
204 }  
205 //=======================================================================
206 //function : Extrema_ExtPRevS
207 //purpose  : 
208 //=======================================================================
209
210 Extrema_ExtPRevS::Extrema_ExtPRevS() 
211 {
212   myS=NULL;
213   myDone = Standard_False;
214 }
215 //=======================================================================
216 //function : Extrema_ExtPRevS
217 //purpose  : 
218 //=======================================================================
219
220 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
221                                    const Adaptor3d_SurfaceOfRevolution& S,
222                                    const Standard_Real Umin,
223                                    const Standard_Real Usup,
224                                    const Standard_Real Vmin,
225                                    const Standard_Real Vsup,
226                                    const Standard_Real TolU,
227                                    const Standard_Real TolV) 
228 {
229   myS=NULL;
230   Initialize (S,
231               Umin, Usup, Vmin, Vsup,
232               TolU, TolV);
233   Perform(P);
234 }
235 //=======================================================================
236 //function : Extrema_ExtPRevS
237 //purpose  : 
238 //=======================================================================
239
240 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
241                                    const Adaptor3d_SurfaceOfRevolution& S,
242                                    const Standard_Real TolU,
243                                    const Standard_Real TolV) 
244 {
245   myS=NULL;
246   Initialize (S,
247               S.FirstUParameter(), S.LastUParameter(),
248               S.FirstVParameter(), S.LastVParameter(),
249               TolU, TolV);
250   Perform(P);
251 }
252 //=======================================================================
253 //function : Initialize
254 //purpose  : 
255 //=======================================================================
256
257 void Extrema_ExtPRevS::Initialize(const Adaptor3d_SurfaceOfRevolution& S,
258                                   const Standard_Real Umin,
259                                   const Standard_Real Usup,
260                                   const Standard_Real Vmin,
261                                   const Standard_Real Vsup,
262                                   const Standard_Real TolU,
263                                   const Standard_Real TolV)
264 {
265   myvinf=Vmin;
266   myvsup=Vsup;
267   mytolv=TolV;
268
269   Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
270   if (!myS || myS != (Adaptor3d_SurfacePtr)&S) {
271     myS = Adaptor3d_SurfacePtr(&S);
272     myPosition = GetPosition(S);
273     myIsAnalyticallyComputable =
274       IsCaseAnalyticallyComputable (anACurve->GetType(),myPosition,S.AxeOfRevolution());
275   }
276   if (!myIsAnalyticallyComputable) {
277     
278     Standard_Integer nbu = 32, nbv = 32;
279
280     if(HasSingularity(S)) {nbv = 100;}
281
282     myExtPS.Initialize(S, nbu, nbv,
283                        Umin, Usup, Vmin, Vsup,
284                        TolU, TolV);
285   }
286 }
287 //=======================================================================
288 //function : Perform
289 //purpose  : 
290 //=======================================================================
291
292 void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
293 {
294   myDone = Standard_False;
295   myNbExt = 0;
296   
297   if (!myIsAnalyticallyComputable) {
298     
299     myExtPS.Perform(P);
300     myDone = myExtPS.IsDone();
301     myNbExt = myExtPS.NbExt();
302     return;
303   }
304   
305   Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
306
307   gp_Ax1 Ax = myS->AxeOfRevolution();
308   gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
309   gp_Pnt O = Ax.Location();
310
311   Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
312   gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
313   if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
314     return;
315   
316   Standard_Real U,V;
317   gp_Pnt P1, Ppp;
318   Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
319   if (Abs(OPpz) <= gp::Resolution()) {
320     Ppp = Pp;
321     U = 0;
322   }
323   else {
324     Ppp = Pp.Translated(Z.Multiplied(-OPpz));
325     if (O.IsEqual(Ppp,Precision::Confusion())) 
326       U = M_PI/2;
327     else {
328       U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
329     }
330   }
331
332   gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
333   if (U != M_PI/2) {
334     if (Abs(OPq.Magnitude()) <= gp::Resolution()) 
335       OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
336     if (OPpp.AngleWithRef(OPq,Dir) < 0)
337       U += M_PI;
338   }
339   
340   gp_Trsf T;
341   T.SetRotation(Ax, -U);
342   P1 = P.Transformed(T);
343   
344   gp_Pnt E;
345   Standard_Real Dist2;
346   Standard_Integer i;
347   
348   Extrema_ExtPElC anExt;
349   PerformExtPElC(anExt, P1, anACurve, mytolv);
350   
351   if (anExt.IsDone()) {
352     myDone = Standard_True;
353     for (i=1; i<=anExt.NbExt(); i++) {
354       Extrema_POnCurv POC=anExt.Point(i);
355       V = POC.Parameter();
356       if (V > myvsup) {
357         //       if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
358         //                         Standard_True, anExt.IsMin(i))) continue;
359         Standard_Real newV = myvsup;
360
361         if((anACurve->GetType() == GeomAbs_Circle) || 
362            (anACurve->GetType() == GeomAbs_Ellipse)) {
363           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
364
365           if (newV > myvsup) {
366             newV -= 2. * M_PI;
367
368             if (newV + mytolv < myvinf) {
369               newV = myvsup;
370             } else if (newV < myvinf) {
371               newV = myvinf;
372             }
373           }
374         }
375         V = newV;
376
377         if ( !IsExtremum (U, V, P, myS, E, Dist2,
378                            Standard_True, anExt.IsMin(i))) {
379           continue;
380         }
381       } else if (V < myvinf) {
382         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
383         //                        Standard_False, anExt.IsMin(i))) continue;
384
385         Standard_Real newV = myvinf;
386
387         if((anACurve->GetType() == GeomAbs_Circle) || 
388            (anACurve->GetType() == GeomAbs_Ellipse)) {
389           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
390           
391           if(newV < myvinf) {
392             newV += 2. * M_PI;
393  
394             if (newV - mytolv > myvsup) {
395               newV = myvinf;
396             } else if (newV > myvsup) {
397               newV = myvsup;
398             }
399           }
400         }
401         V = newV;
402
403         if ( !IsExtremum (U, V, P, myS, E, Dist2,
404                           Standard_False, anExt.IsMin(i))) continue;
405       } else {
406         E = myS->Value(U,V);
407         Dist2 = P.SquareDistance(E);
408       }
409       if (IsOriginalPnt(E, myPoint, myNbExt)) {
410         myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
411         mySqDist[myNbExt] = Dist2;
412       }
413     }
414   }
415   T.SetRotation(Ax, M_PI);
416   P1.Transform(T);
417   
418   PerformExtPElC(anExt, P1, anACurve, mytolv);
419   if (anExt.IsDone()) {
420     myDone = Standard_True;
421
422     U += M_PI;
423     
424     for (i=1; i<=anExt.NbExt(); i++) {
425       Extrema_POnCurv POC=anExt.Point(i);
426       V = POC.Parameter();
427       if (V > myvsup) {
428         //      if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
429         //                         Standard_True, anExt.IsMin(i))) continue;
430
431         Standard_Real newV = myvsup;
432
433         if((anACurve->GetType() == GeomAbs_Circle) || 
434            (anACurve->GetType() == GeomAbs_Ellipse)) {
435           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
436
437           if (newV > myvsup) {
438             newV -= 2. * M_PI;
439
440             if (newV + mytolv < myvinf) {
441               newV = myvsup;
442             } else if (newV < myvinf) {
443               newV = myvinf;
444             }
445           }
446         }
447         V = newV;
448         
449         if ( !IsExtremum (U, V, P, myS, E, Dist2,
450                            Standard_True, anExt.IsMin(i))) continue;
451       } else if (V < myvinf) {
452         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
453         //                        Standard_False, anExt.IsMin(i))) continue;
454         Standard_Real newV = myvinf;
455
456         if((anACurve->GetType() == GeomAbs_Circle) || 
457            (anACurve->GetType() == GeomAbs_Ellipse)) {
458           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
459           
460           if(newV < myvinf) {
461             newV += 2. * M_PI;
462  
463             if (newV - mytolv > myvsup) {
464               newV = myvinf;
465             } else if (newV > myvsup) {
466               newV = myvsup;
467             }
468           }
469         }
470         V = newV;
471
472         if ( !IsExtremum (U, V, P, myS, E, Dist2,
473                           Standard_False, anExt.IsMin(i))) continue;
474       } else {
475         E = myS->Value(U,V);
476         Dist2 = P.SquareDistance(E);
477       }
478       if (IsOriginalPnt(E, myPoint, myNbExt)) {
479         myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
480         mySqDist[myNbExt] = Dist2;
481       
482       }
483     }
484   }
485 }
486
487
488 //=======================================================================
489 //function : IsDone
490 //purpose  : 
491 //=======================================================================
492
493 Standard_Boolean Extrema_ExtPRevS::IsDone() const
494 {
495   return myDone; 
496 }
497
498
499 //=======================================================================
500 //function : NbExt
501 //purpose  : 
502 //=======================================================================
503
504 Standard_Integer Extrema_ExtPRevS::NbExt() const
505 {
506   if (!IsDone()) { StdFail_NotDone::Raise(); }
507   return myNbExt;
508 }
509
510
511 //=======================================================================
512 //function : Value
513 //purpose  : 
514 //=======================================================================
515
516 Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
517 {
518   if (!IsDone()) { StdFail_NotDone::Raise(); }
519   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
520   if (myIsAnalyticallyComputable)
521     return mySqDist[N];
522   else
523     return myExtPS.SquareDistance(N);
524 }
525 //=======================================================================
526 //function : Point
527 //purpose  : 
528 //=======================================================================
529
530 Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
531 {
532   if (!IsDone()) { StdFail_NotDone::Raise(); }
533   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
534   if (myIsAnalyticallyComputable)
535     return myPoint[N];
536   else
537     return myExtPS.Point(N);
538 }
539
540