0023733: PCurve for edge on face creation failure
[occt.git] / src / Extrema / Extrema_ExtPElC.cxx
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
19 #include <Extrema_ExtPElC.ixx>
20 #include <StdFail_NotDone.hxx>
21 #include <math_DirectPolynomialRoots.hxx>
22 #include <math_TrigonometricFunctionRoots.hxx>
23 #include <ElCLib.hxx>
24 #include <Standard_OutOfRange.hxx>
25 #include <Standard_NotImplemented.hxx>
26 #include <Precision.hxx>
27
28
29 //=============================================================================
30
31 Extrema_ExtPElC::Extrema_ExtPElC () { myDone = Standard_False; }
32 //=============================================================================
33
34 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P, 
35                                   const gp_Lin&       L,
36                                   const Standard_Real Tol,
37                                   const Standard_Real Uinf, 
38                                   const Standard_Real Usup)
39 {
40   Perform(P, L, Tol, Uinf, Usup);
41 }
42
43 void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
44                               const gp_Lin&       L,
45                               const Standard_Real Tol,
46                               const Standard_Real Uinf, 
47                               const Standard_Real Usup)
48 {
49   myDone = Standard_False;
50   myNbExt = 0;
51   gp_Vec V1 = gp_Vec(L.Direction());
52   gp_Pnt OR = L.Location();
53   gp_Vec V(OR, P);
54   Standard_Real Mydist = V1.Dot(V);
55   if ((Mydist >= Uinf-Tol) && 
56       (Mydist <= Usup+Tol)){ 
57     
58     gp_Pnt MyP = OR.Translated(Mydist*V1);
59     Extrema_POnCurv MyPOnCurve(Mydist, MyP);
60     mySqDist[0] = P.SquareDistance(MyP);
61     myPoint[0] = MyPOnCurve;
62     myIsMin[0] = Standard_True;
63     myNbExt = 1;
64     myDone = Standard_True;
65   }
66
67 }
68
69
70
71 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
72                                   const gp_Circ&      C,
73                                   const Standard_Real Tol,
74                                   const Standard_Real Uinf,
75                                   const Standard_Real Usup)
76 {
77   Perform(P, C, Tol, Uinf, Usup);
78 }
79
80 void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
81                               const gp_Circ&      C,
82                               const Standard_Real Tol,
83                               const Standard_Real Uinf,
84                               const Standard_Real Usup)
85 /*-----------------------------------------------------------------------------
86 Function:
87   Find values of parameter u such as:
88    - dist(P,C(u)) pass by an extrema,
89    - Uinf <= u <= Usup.
90
91 Method:
92   Pass 3 stages:
93   1- Projection of point P in the plane of the circle,
94   2- Calculation of u solutions in [0.,2.*M_PI]:
95       Let Pp, the projected point and 
96            O, the center of the circle;
97      2 cases:
98      - if Pp is mixed with 0, there is an infinite number of solutions;
99        IsDone() renvoie Standard_False.
100      - otherwise, 2 points are solutions for the complete circle:
101        . Us1 = angle(OPp,OX) corresponds to the minimum,
102        . let Us2 = ( Us1 + M_PI if Us1 < M_PI,
103                     ( Us1 - M_PI otherwise;
104          Us2 corresponds to the maximum.
105   3- Calculate the extrema in [Uinf,Usup].
106 -----------------------------------------------------------------------------*/
107 {
108   myDone = Standard_False;
109   myNbExt = 0;
110
111 // 1- Projection of the point P in the plane of circle -> Pp ...
112
113   gp_Pnt O = C.Location();
114   gp_Vec Axe (C.Axis().Direction());
115   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
116   gp_Pnt Pp = P.Translated(Trsl);
117
118 // 2- Calculate u solutions in [0.,2.*PI] ...
119
120   gp_Vec OPp (O,Pp);
121   if (OPp.Magnitude() < Tol) { return; }
122   Standard_Real Usol[2];
123   Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -M_PI<U1<M_PI
124   Usol[1] = Usol[0] + M_PI;
125
126   Standard_Real myuinf = Uinf;
127   //modified by NIZNHY-PKV Fri Apr 20 15:03:28 2001 f
128   //Standard_Real TolU = Tol*C.Radius();
129   Standard_Real TolU, aR;
130   aR=C.Radius();
131   TolU=Precision::Infinite();
132   if (aR > gp::Resolution()) {
133     TolU= Tol/aR;
134   }
135   //modified by NIZNHY-PKV Fri Apr 20 15:03:32 2001 t
136   ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[0]);
137   ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[1]);
138   if (((Usol[0]-2*M_PI-Uinf) < TolU) && ((Usol[0]-2*M_PI-Uinf) > -TolU)) Usol[0] = Uinf;
139   if (((Usol[1]-2*M_PI-Uinf) < TolU) && ((Usol[1]-2*M_PI-Uinf) > -TolU)) Usol[1] = Uinf;
140
141
142 // 3- Calculate extrema in [Umin,Umax] ...
143
144   gp_Pnt Cu;
145   Standard_Real Us;
146   for (Standard_Integer NoSol = 0; NoSol <= 1; NoSol++) {
147     Us = Usol[NoSol];
148     if (((Uinf-Us) < TolU) && ((Us-Usup) < TolU)) {
149       Cu = ElCLib::Value(Us,C);
150       mySqDist[myNbExt] = Cu.SquareDistance(P);
151       myIsMin[myNbExt] = (NoSol == 0);
152       myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
153       myNbExt++;
154     }
155   }
156   myDone = Standard_True;
157 }
158 //=============================================================================
159
160 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
161                                   const gp_Elips&     C,
162                                   const Standard_Real Tol,
163                                   const Standard_Real Uinf, 
164                                   const Standard_Real Usup)
165 {
166   Perform(P, C, Tol, Uinf, Usup);
167 }
168
169
170
171 void Extrema_ExtPElC::Perform (const gp_Pnt&       P, 
172                                const gp_Elips&     C,
173                                const Standard_Real Tol,
174                                const Standard_Real Uinf, 
175                                const Standard_Real Usup)
176 /*-----------------------------------------------------------------------------
177 Function:
178   Find values of parameter u so that:
179    - dist(P,C(u)) passes by an extremum,
180    - Uinf <= u <= Usup.
181
182 Method:
183   Takes 2 steps:
184   1- Projection of point P in the plane of the ellipse,
185   2- Calculation of the solutions:
186      Let Pp, the projected point; find values u so that:
187        (C(u)-Pp).C'(u) = 0. (1)
188      Let Cos = cos(u) and Sin = sin(u),
189           C(u) = (A*Cos,B*Sin) and Pp = (X,Y);
190      Then, (1) <=> (A*Cos-X,B*Sin-Y).(-A*Sin,B*Cos) = 0.
191                     (B**2-A**2)*Cos*Sin - B*Y*Cos + A*X*Sin = 0.
192      Use algorithm math_TrigonometricFunctionRoots to solve this equation.
193 -----------------------------------------------------------------------------*/
194 {
195   myDone = Standard_False;
196   myNbExt = 0;
197
198 // 1- Projection of point P in the plane of the ellipse -> Pp ...
199
200   gp_Pnt O = C.Location();
201   gp_Vec Axe (C.Axis().Direction());
202   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
203   gp_Pnt Pp = P.Translated(Trsl);
204
205 // 2- Calculation of solutions ...
206
207   Standard_Integer NoSol, NbSol;
208   Standard_Real A = C.MajorRadius();
209   Standard_Real B = C.MinorRadius();
210   gp_Vec OPp (O,Pp);
211   Standard_Real OPpMagn = OPp.Magnitude();
212   if (OPpMagn < Tol) { if (Abs(A-B) < Tol) { return; } }
213   Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
214   Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
215   //  Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
216
217   Standard_Real ko2 = (B*B-A*A)/2., ko3 = -B*Y, ko4 = A*X;
218   if(Abs(ko3) < 1.e-16*Max(Abs(ko2), Abs(ko3))) ko3 = 0.0;
219
220 //  math_TrigonometricFunctionRoots Sol(0.,(B*B-A*A)/2.,-B*Y,A*X,0.,Uinf,Usup);
221   math_TrigonometricFunctionRoots Sol(0.,ko2, ko3, ko4, 0.,Uinf,Usup);
222
223   if (!Sol.IsDone()) { return; }
224   gp_Pnt Cu;
225   Standard_Real Us;
226   NbSol = Sol.NbSolutions();
227   for (NoSol = 1; NoSol <= NbSol; NoSol++) {
228     Us = Sol.Value(NoSol);
229     Cu = ElCLib::Value(Us,C);
230     mySqDist[myNbExt] = Cu.SquareDistance(P);
231     myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
232     Cu = ElCLib::Value(Us + 0.1, C);
233     myIsMin[myNbExt] = mySqDist[myNbExt] < Cu.SquareDistance(P);
234     myNbExt++;
235   }
236   myDone = Standard_True;
237 }
238 //=============================================================================
239
240 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P, 
241                                   const gp_Hypr&      C,
242                                   const Standard_Real Tol,
243                                   const Standard_Real Uinf,
244                                   const Standard_Real Usup)
245 {
246   Perform(P, C, Tol, Uinf, Usup);
247 }
248
249
250 void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
251                               const gp_Hypr&      C,
252                               const Standard_Real Tol,
253                               const Standard_Real Uinf,
254                               const Standard_Real Usup)
255 /*-----------------------------------------------------------------------------
256 Function:
257   Find values of parameter u so that:
258    - dist(P,C(u)) passes by an extremum,
259    - Uinf <= u <= Usup.
260
261 Method:
262   Takes 2 steps:
263   1- Projection of point P in the plane of the hyperbola,
264   2- Calculation of solutions:
265       Let Pp, le point projete; on recherche les valeurs u telles que:
266        (C(u)-Pp).C'(u) = 0. (1)
267       Let R and r be the radiuses of the hyperbola,
268           Chu = Cosh(u) and Shu = Sinh(u),
269           C(u) = (R*Chu,r*Shu) and Pp = (X,Y);
270      Then, (1) <=> (R*Chu-X,r*Shu-Y).(R*Shu,r*Chu) = 0.
271                     (R**2+r**2)*Chu*Shu - X*R*Shu - Y*r*Chu = 0. (2)
272      Let v = e**u;
273      Then, by using Chu = (e**u+e**(-u))/2. and Sh = (e**u-e**(-u)))/2.
274             (2) <=> ((R**2+r**2)/4.) * (v**2-v**(-2)) -
275                     ((X*R+Y*r)/2.) * v +
276                     ((X*R-Y*r)/2.) * v**(-1) = 0.
277       (2)* v**2 <=> ((R**2+r**2)/4.) * v**4 -
278                      ((X*R+Y*r)/2.)  * v**3 +
279                      ((X*R-Y*r)/2.)  * v    -
280                     ((R**2+r**2)/4.)        = 0.
281   Use algorithm math_DirectPolynomialRoots to solve this equation by v.
282 -----------------------------------------------------------------------------*/
283 {
284   myDone = Standard_False;
285   myNbExt = 0;
286
287 // 1- Projection of point P in the plane of hyperbola -> Pp ...
288
289   gp_Pnt O = C.Location();
290   gp_Vec Axe (C.Axis().Direction());
291   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
292   gp_Pnt Pp = P.Translated(Trsl);
293
294 // 2- Calculation of solutions ...
295
296   Standard_Real Tol2 = Tol * Tol;
297   Standard_Real R = C.MajorRadius();
298   Standard_Real r = C.MinorRadius();
299   gp_Vec OPp (O,Pp);
300   Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
301   Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
302
303   Standard_Real C1 = (R*R+r*r)/4.;
304   math_DirectPolynomialRoots Sol(C1,-(X*R+Y*r)/2.,0.,(X*R-Y*r)/2.,-C1);
305   if (!Sol.IsDone()) { return; }
306   gp_Pnt Cu;
307   Standard_Real Us, Vs;
308   Standard_Integer NbSol = Sol.NbSolutions();
309   Standard_Boolean DejaEnr;
310   Standard_Integer NoExt;
311   gp_Pnt TbExt[4];
312   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
313     Vs = Sol.Value(NoSol);
314     if (Vs > 0.) {
315       Us = Log(Vs);
316       if ((Us >= Uinf) && (Us <= Usup)) {
317         Cu = ElCLib::Value(Us,C);
318         DejaEnr = Standard_False;
319         for (NoExt = 0; NoExt < myNbExt; NoExt++) {
320           if (TbExt[NoExt].SquareDistance(Cu) < Tol2) {
321             DejaEnr = Standard_True;
322             break;
323           }
324         }
325         if (!DejaEnr) {
326           TbExt[myNbExt] = Cu;
327           mySqDist[myNbExt] = Cu.SquareDistance(P);
328           myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
329           myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
330           myNbExt++;
331         }
332       } // if ((Us >= Uinf) && (Us <= Usup))
333     } // if (Vs > 0.)
334   } // for (Standard_Integer NoSol = 1; ...
335   myDone = Standard_True;
336 }
337 //=============================================================================
338
339 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
340                                   const gp_Parab&     C,
341                                   const Standard_Real Tol,
342                                   const Standard_Real Uinf,
343                                   const Standard_Real Usup)
344 {
345   Perform(P, C, Tol, Uinf, Usup);
346 }
347
348
349 void Extrema_ExtPElC::Perform(const gp_Pnt&       P, 
350                               const gp_Parab&     C,
351 //                            const Standard_Real Tol,
352                               const Standard_Real ,
353                               const Standard_Real Uinf,
354                               const Standard_Real Usup)
355 /*-----------------------------------------------------------------------------
356 Function:
357   Find values of parameter u so that:
358    - dist(P,C(u)) pass by an extremum,
359    - Uinf <= u <= Usup.
360
361 Method:
362   Takes 2 steps:
363   1- Projection of point P in the plane of the parabola,
364   2- Calculation of solutions:
365      Let Pp, the projected point; find values u so that:
366        (C(u)-Pp).C'(u) = 0. (1)
367      Let F the focus of the parabola,
368           C(u) = ((u*u)/(4.*F),u) and Pp = (X,Y);
369      Alors, (1) <=> ((u*u)/(4.*F)-X,u-Y).(u/(2.*F),1) = 0.
370                     (1./(4.*F)) * U**3 + (2.*F-X) * U - 2*F*Y = 0.
371     Use algorithm math_DirectPolynomialRoots to solve this equation by U.
372 -----------------------------------------------------------------------------*/
373 {
374   myDone = Standard_False;
375   myNbExt = 0;
376
377 // 1- Projection of point P in the plane of the parabola -> Pp ...
378
379   gp_Pnt O = C.Location();
380   gp_Vec Axe (C.Axis().Direction());
381   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
382   gp_Pnt Pp = P.Translated(Trsl);
383
384 // 2- Calculation of solutions ...
385
386   Standard_Real F = C.Focal();
387   gp_Vec OPp (O,Pp);
388   Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
389 //  Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
390   Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
391   math_DirectPolynomialRoots Sol(1./(4.*F),0.,2.*F-X,-2.*F*Y);
392   if (!Sol.IsDone()) { return; }
393   gp_Pnt Cu;
394   Standard_Real Us;
395   Standard_Integer NbSol = Sol.NbSolutions();
396   Standard_Boolean DejaEnr;
397   Standard_Integer NoExt;
398   gp_Pnt TbExt[3];
399   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
400     Us = Sol.Value(NoSol);
401     if ((Us >= Uinf) && (Us <= Usup)) {
402       Cu = ElCLib::Value(Us,C);
403       DejaEnr = Standard_False;
404       for (NoExt = 0; NoExt < myNbExt; NoExt++) {
405     if (TbExt[NoExt].SquareDistance(Cu) < Precision::SquareConfusion()) {
406           DejaEnr = Standard_True;
407           break;
408         }
409       }
410       if (!DejaEnr) {
411         TbExt[myNbExt] = Cu;
412         mySqDist[myNbExt] = Cu.SquareDistance(P);
413         myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
414         myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
415         myNbExt++;
416       }
417     } // if ((Us >= Uinf) && (Us <= Usup))
418   } // for (Standard_Integer NoSol = 1; ...
419   myDone = Standard_True;
420 }
421 //=============================================================================
422
423 Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
424 //=============================================================================
425
426 Standard_Integer Extrema_ExtPElC::NbExt () const
427 {
428   if (!IsDone()) { StdFail_NotDone::Raise(); }
429   return myNbExt;
430 }
431 //=============================================================================
432
433 Standard_Real Extrema_ExtPElC::SquareDistance (const Standard_Integer N) const
434 {
435   if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
436   return mySqDist[N-1];
437 }
438 //=============================================================================
439
440 Standard_Boolean Extrema_ExtPElC::IsMin (const Standard_Integer N) const
441 {
442   if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
443   return myIsMin[N-1];
444 }
445 //=============================================================================
446
447 Extrema_POnCurv Extrema_ExtPElC::Point (const Standard_Integer N) const
448 {
449   if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
450   return myPoint[N-1];
451 }
452 //=============================================================================