0023733: PCurve for edge on face creation failure
[occt.git] / src / Extrema / Extrema_ExtPElC.cxx
CommitLineData
b311480e 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
7fd59977 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
31Extrema_ExtPElC::Extrema_ExtPElC () { myDone = Standard_False; }
32//=============================================================================
33
34Extrema_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
43void 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
71Extrema_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
80void 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/*-----------------------------------------------------------------------------
0d969553
Y
86Function:
87 Find values of parameter u such as:
88 - dist(P,C(u)) pass by an extrema,
7fd59977 89 - Uinf <= u <= Usup.
90
0d969553
Y
91Method:
92 Pass 3 stages:
93 1- Projection of point P in the plane of the circle,
c6541a0c 94 2- Calculation of u solutions in [0.,2.*M_PI]:
0d969553
Y
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;
7fd59977 99 IsDone() renvoie Standard_False.
0d969553
Y
100 - otherwise, 2 points are solutions for the complete circle:
101 . Us1 = angle(OPp,OX) corresponds to the minimum,
c6541a0c
D
102 . let Us2 = ( Us1 + M_PI if Us1 < M_PI,
103 ( Us1 - M_PI otherwise;
0d969553
Y
104 Us2 corresponds to the maximum.
105 3- Calculate the extrema in [Uinf,Usup].
7fd59977 106-----------------------------------------------------------------------------*/
107{
108 myDone = Standard_False;
109 myNbExt = 0;
110
0d969553 111// 1- Projection of the point P in the plane of circle -> Pp ...
7fd59977 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
0d969553 118// 2- Calculate u solutions in [0.,2.*PI] ...
7fd59977 119
120 gp_Vec OPp (O,Pp);
121 if (OPp.Magnitude() < Tol) { return; }
122 Standard_Real Usol[2];
c6541a0c
D
123 Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -M_PI<U1<M_PI
124 Usol[1] = Usol[0] + M_PI;
7fd59977 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
c6541a0c
D
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;
7fd59977 140
141
0d969553 142// 3- Calculate extrema in [Umin,Umax] ...
7fd59977 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
160Extrema_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
171void 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/*-----------------------------------------------------------------------------
0d969553
Y
177Function:
178 Find values of parameter u so that:
179 - dist(P,C(u)) passes by an extremum,
7fd59977 180 - Uinf <= u <= Usup.
181
0d969553
Y
182Method:
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:
7fd59977 187 (C(u)-Pp).C'(u) = 0. (1)
0d969553
Y
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.
7fd59977 191 (B**2-A**2)*Cos*Sin - B*Y*Cos + A*X*Sin = 0.
0d969553 192 Use algorithm math_TrigonometricFunctionRoots to solve this equation.
7fd59977 193-----------------------------------------------------------------------------*/
194{
195 myDone = Standard_False;
196 myNbExt = 0;
197
0d969553 198// 1- Projection of point P in the plane of the ellipse -> Pp ...
7fd59977 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
0d969553 205// 2- Calculation of solutions ...
7fd59977 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);
7fd59977 231 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
a20400b3 232 Cu = ElCLib::Value(Us + 0.1, C);
233 myIsMin[myNbExt] = mySqDist[myNbExt] < Cu.SquareDistance(P);
7fd59977 234 myNbExt++;
235 }
236 myDone = Standard_True;
237}
238//=============================================================================
239
240Extrema_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
250void 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/*-----------------------------------------------------------------------------
0d969553
Y
256Function:
257 Find values of parameter u so that:
258 - dist(P,C(u)) passes by an extremum,
7fd59977 259 - Uinf <= u <= Usup.
260
0d969553
Y
261Method:
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:
7fd59977 266 (C(u)-Pp).C'(u) = 0. (1)
0d969553
Y
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.
7fd59977 271 (R**2+r**2)*Chu*Shu - X*R*Shu - Y*r*Chu = 0. (2)
0d969553
Y
272 Let v = e**u;
273 Then, by using Chu = (e**u+e**(-u))/2. and Sh = (e**u-e**(-u)))/2.
7fd59977 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.
0d969553 281 Use algorithm math_DirectPolynomialRoots to solve this equation by v.
7fd59977 282-----------------------------------------------------------------------------*/
283{
284 myDone = Standard_False;
285 myNbExt = 0;
286
0d969553 287// 1- Projection of point P in the plane of hyperbola -> Pp ...
7fd59977 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
0d969553 294// 2- Calculation of solutions ...
7fd59977 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);
7fd59977 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
339Extrema_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
349void 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/*-----------------------------------------------------------------------------
0d969553
Y
356Function:
357 Find values of parameter u so that:
358 - dist(P,C(u)) pass by an extremum,
7fd59977 359 - Uinf <= u <= Usup.
360
0d969553
Y
361Method:
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:
7fd59977 366 (C(u)-Pp).C'(u) = 0. (1)
0d969553
Y
367 Let F the focus of the parabola,
368 C(u) = ((u*u)/(4.*F),u) and Pp = (X,Y);
7fd59977 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.
0d969553 371 Use algorithm math_DirectPolynomialRoots to solve this equation by U.
7fd59977 372-----------------------------------------------------------------------------*/
373{
374 myDone = Standard_False;
375 myNbExt = 0;
376
0d969553 377// 1- Projection of point P in the plane of the parabola -> Pp ...
7fd59977 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
0d969553 384// 2- Calculation of solutions ...
7fd59977 385
386 Standard_Real F = C.Focal();
387 gp_Vec OPp (O,Pp);
7fd59977 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++) {
08cd2f6b 405 if (TbExt[NoExt].SquareDistance(Cu) < Precision::SquareConfusion()) {
7fd59977 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
423Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
424//=============================================================================
425
426Standard_Integer Extrema_ExtPElC::NbExt () const
427{
428 if (!IsDone()) { StdFail_NotDone::Raise(); }
429 return myNbExt;
430}
431//=============================================================================
432
433Standard_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
440Standard_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
447Extrema_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//=============================================================================