projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023952: Improving thread-safety of intersections, approximations and other modeling...
[occt.git] / src / IntAna / IntAna_Curve.cxx
CommitLineData
b311480e 1// Created on: 1992-06-30
2// Created by: Laurent BUCHARD
3// Copyright (c) 1992-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
7fd59977 21
22#ifndef DEB
23#define No_Standard_RangeError
24#define No_Standard_OutOfRange
25#endif
26
27//----------------------------------------------------------------------
28//-- Differents constructeurs sont proposes qui correspondent aux
29//-- polynomes en Z :
30//-- A(Sin(Theta),Cos(Theta)) Z**2
31//-- + B(Sin(Theta),Cos(Theta)) Z
32//-- + C(Sin(Theta),Cos(Theta))
33//--
34//-- Une Courbe est definie sur un domaine
35//--
36//-- Value retourne le point de parametre U(Theta),V(Theta)
37//-- ou V est la solution du polynome A V**2 + B V + C
38//-- (Selon les cas, on prend V+ ou V-)
39//--
40//-- D1u calcule le vecteur tangent a la courbe
41//-- et retourne le booleen Standard_False si ce calcul ne peut
42//-- pas etre mene a bien.
43//----------------------------------------------------------------------
44
45#include <Precision.hxx>
46
47#include <IntAna_Curve.ixx>
48
49#include <Standard_DomainError.hxx>
50#include <math_DirectPolynomialRoots.hxx>
51#include <ElSLib.hxx>
52#include <gp_XYZ.hxx>
53
54//=======================================================================
55//function : IntAna_Curve
56//purpose :
57//=======================================================================
58 IntAna_Curve::IntAna_Curve()
59{
60 typequadric=GeomAbs_OtherSurface;
61 firstbounded=Standard_False;
62 lastbounded=Standard_False;
63}
64//=======================================================================
65//function : SetConeQuadValues
66//purpose : Description de l intersection Cone Quadrique
67//=======================================================================
68 void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
69 const Standard_Real Qxx,
70 const Standard_Real Qyy,
71 const Standard_Real Qzz,
72 const Standard_Real Qxy,
73 const Standard_Real Qxz,
74 const Standard_Real Qyz,
75 const Standard_Real Qx,
76 const Standard_Real Qy,
77 const Standard_Real Qz,
78 const Standard_Real Q1,
79 const Standard_Real TOL,
80 const Standard_Real DomInf,
81 const Standard_Real DomSup,
82 const Standard_Boolean twocurves,
83 const Standard_Boolean takezpositive)
84{
85
86 Ax3 = Cone.Position();
87 RCyl = Cone.RefRadius();
88
89 Angle = Cone.SemiAngle();
90 Standard_Real UnSurTgAngle = 1.0/(Tan(Cone.SemiAngle()));
91
92 typequadric= GeomAbs_Cone;
93
94 TwoCurves = twocurves; //-- deux Z pour un meme parametre
95 TakeZPositive = takezpositive; //-- Prendre sur la courbe le Z Positif
96 //-- ( -B + Sqrt()) et non (-B - Sqrt())
97
98
99 Z0Cte = Q1; //-- Attention On a Z?Cos Cos(t)
100 Z0Sin = 0.0; //-- et Non 2 Z?Cos Cos(t) !!!
101 Z0Cos = 0.0; //-- Ce pour tous les Parametres
102 Z0CosCos = 0.0; //-- ie pas de Coefficient 2
103 Z0SinSin = 0.0; //-- devant les termes CS C S
104 Z0CosSin = 0.0;
105
106 Z1Cte = 2.0*(UnSurTgAngle)*Qz;
107 Z1Sin = Qy+Qy;
108 Z1Cos = Qx+Qx;
109 Z1CosCos = 0.0;
110 Z1SinSin = 0.0;
111 Z1CosSin = 0.0;
112
113 Z2Cte = Qzz * UnSurTgAngle*UnSurTgAngle;
114 Z2Sin = (UnSurTgAngle+UnSurTgAngle)*Qyz;
115 Z2Cos = (UnSurTgAngle+UnSurTgAngle)*Qxz;
116 Z2CosCos = Qxx;
117 Z2SinSin = Qyy;
118 Z2CosSin = Qxy+Qxy;
119
120 Tolerance = TOL;
121 DomainInf = DomInf;
122 DomainSup = DomSup;
123
124 RestrictedInf = RestrictedSup = Standard_True; //-- Le Domaine est Borne
125 firstbounded = lastbounded = Standard_False;
126}
127
128//=======================================================================
129//function : SetCylinderQuadValues
130//purpose : Description de l intersection Cylindre Quadrique
131//=======================================================================
132 void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
133 const Standard_Real Qxx,
134 const Standard_Real Qyy,
135 const Standard_Real Qzz,
136 const Standard_Real Qxy,
137 const Standard_Real Qxz,
138 const Standard_Real Qyz,
139 const Standard_Real Qx,
140 const Standard_Real Qy,
141 const Standard_Real Qz,
142 const Standard_Real Q1,
143 const Standard_Real TOL,
144 const Standard_Real DomInf,
145 const Standard_Real DomSup,
146 const Standard_Boolean twocurves,
147 const Standard_Boolean takezpositive)
148{
149
150 Ax3 = Cyl.Position();
151 RCyl = Cyl.Radius();
152 typequadric= GeomAbs_Cylinder;
153
154 TwoCurves = twocurves; //-- deux Z pour un meme parametre
155 TakeZPositive = takezpositive; //-- Prendre sur la courbe le Z Positif
156 Standard_Real RCylmul2 = RCyl+RCyl; //-- ( -B + Sqrt())
157
158 Z0Cte = Q1;
159 Z0Sin = RCylmul2*Qy;
160 Z0Cos = RCylmul2*Qx;
161 Z0CosCos = Qxx*RCyl*RCyl;
162 Z0SinSin = Qyy*RCyl*RCyl;
163 Z0CosSin = RCylmul2*RCyl*Qxy;
164
165 Z1Cte = Qz+Qz;
166 Z1Sin = RCylmul2*Qyz;
167 Z1Cos = RCylmul2*Qxz;
168 Z1CosCos = 0.0;
169 Z1SinSin = 0.0;
170 Z1CosSin = 0.0;
171
172 Z2Cte = Qzz;
173 Z2Sin = 0.0;
174 Z2Cos = 0.0;
175 Z2CosCos = 0.0;
176 Z2SinSin = 0.0;
177 Z2CosSin = 0.0;
178
179 Tolerance = TOL;
180 DomainInf = DomInf;
181 DomainSup = DomSup;
182
183 RestrictedInf = RestrictedSup = Standard_True;
184 firstbounded = lastbounded = Standard_False;
185}
186
187//=======================================================================
188//function : IsOpen
189//purpose :
190//=======================================================================
191 Standard_Boolean IntAna_Curve::IsOpen() const
192{
193 return(RestrictedInf && RestrictedSup);
194}
195
196//=======================================================================
197//function : Domain
198//purpose :
199//=======================================================================
200 void IntAna_Curve::Domain(Standard_Real& DInf,
201 Standard_Real& DSup) const
202{
203 if(RestrictedInf && RestrictedSup) {
204 DInf=DomainInf;
205 DSup=DomainSup;
206 if(TwoCurves) {
207 DSup+=DSup-DInf;
208 }
209 }
210 else {
211 Standard_DomainError::Raise("IntAna_Curve::Domain");
212 }
213}
214//=======================================================================
215//function : IsConstant
216//purpose :
217//=======================================================================
218 Standard_Boolean IntAna_Curve::IsConstant() const
219{
220 //-- ??? Pas facile de decider a la seule vue des Param.
221 return(Standard_False);
222}
223
224//=======================================================================
225//function : IsFirstOpen
226//purpose :
227//=======================================================================
228 Standard_Boolean IntAna_Curve::IsFirstOpen() const
229{
230 return(firstbounded);
231}
232
233//=======================================================================
234//function : IsLastOpen
235//purpose :
236//=======================================================================
237 Standard_Boolean IntAna_Curve::IsLastOpen() const
238{
239 return(lastbounded);
240}
241//=======================================================================
242//function : SetIsFirstOpen
243//purpose :
244//=======================================================================
245 void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
246{
247 firstbounded = Flag;
248}
249
250//=======================================================================
251//function : SetIsLastOpen
252//purpose :
253//=======================================================================
254 void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
255{
256 lastbounded = Flag;
257}
258
259//=======================================================================
260//function : InternalUVValue
261//purpose :
262//=======================================================================
263 void IntAna_Curve::InternalUVValue(const Standard_Real theta,
264 Standard_Real& Param1,
265 Standard_Real& Param2,
266 Standard_Real& A,
267 Standard_Real& B,
268 Standard_Real& C,
269 Standard_Real& cost,
270 Standard_Real& sint,
271 Standard_Real& SigneSqrtDis) const
272{
273 Standard_Real Theta=theta;
274 Standard_Boolean SecondSolution=Standard_False;
275
276 if((Theta<DomainInf) ||
277 ((Theta>DomainSup) && (!TwoCurves)) ||
278 (Theta>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
279 Standard_DomainError::Raise("IntAna_Curve::Domain");
280 }
281
282 if(Theta>DomainSup) {
283 Theta=DomainSup+DomainSup-Theta;
284 SecondSolution=Standard_True;
285 }
286
287 Param1=Theta;
288
289 if(!TwoCurves) {
290 SecondSolution=TakeZPositive;
291 }
292 //
293 cost = Cos(Theta);
294 sint = Sin(Theta);
295 Standard_Real costsint = cost*sint;
296
297 A=Z2Cte+sint*(Z2Sin+sint*Z2SinSin)+cost*(Z2Cos+cost*Z2CosCos)
298 +Z2CosSin*costsint;
299
300 B=Z1Cte+sint*(Z1Sin+sint*Z1SinSin)+cost*(Z1Cos+cost*Z1CosCos)
301 +Z1CosSin*costsint;
302
303 C=Z0Cte+sint*(Z0Sin+sint*Z0SinSin)+cost*(Z0Cos+cost*Z0CosCos)
304 +Z0CosSin*costsint;
305
306
307 Standard_Real Discriminant = B*B-4.0*A*C;
308
309 if(Abs(A)<=0.000000001) {
310 //-- cout<<" IntAna_Curve:: Internal UV Value : A="<<A<<" -> Abs(A)="<<Abs(A)<<endl;
311 if(Abs(B)<=0.0000000001) {
312 //-- cout<<" Probleme : Pas de solutions "<<endl;
313 Param2=0.0;
314 }
315 else {
316 //modified by NIZNHY-PKV Fri Dec 2 16:02:46 2005f
317 Param2 = -C/B;
318 /*
319 if(!SecondSolution) {
320 //-- Cas Param2 = (-B+Sqrt(Discriminant))/(A+A);
321 //-- = (-B+Sqrt(B**2 - Eps)) / 2A
322 //-- = -C / B
323 Param2 = -C/B;
324 }
325 else {
326 //-- Cas Param2 = (-B-Sqrt(Discriminant))/(A+A);
327 //-- = (-B-Sqrt(B**2 - Eps)) / 2A
328 if(A) {
329 Param2 = -B/A;
330 }
331 else {
332 Param2 = -B*10000000.0;
333 }
334 }
335 */
336 //modified by NIZNHY-PKV Fri Dec 2 16:02:54 2005t
337 }
338 }
339 else {
340 if(Discriminant<=0.0000000001 ||
341 Abs(Discriminant/(4*A))<=0.0000000001) Discriminant=0.0;
342 SigneSqrtDis = (SecondSolution)? Sqrt(Discriminant)
343 : -Sqrt(Discriminant);
344 Param2=(-B+SigneSqrtDis)/(A+A);
345 }
346}
347
348//=======================================================================
349//function : Value
350//purpose :
351//=======================================================================
352 gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
353{
354 Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
355 //
356 A=0.0; B=0.0; C=0.0;
357 U=0.0; V=0.0;
358 sint=0.0; cost=0.0;
359 SigneSqrtDis=0.0;
360 InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
361 //-- checked the parameter U and Raises Domain Error if Error
362 return(InternalValue(U,V));
363}
364//=======================================================================
365//function : D1u
366//purpose :
367//=======================================================================
368 Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
369 gp_Pnt& Pt,
370 gp_Vec& Vec)
371{
372 //-- Pour detecter le cas ou le calcul est impossible
373 Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
374 A=0.0; B=0.0; C=0.0;
375 U=0.0; V=0.0;
376 sint=0.0; cost=0.0;
377 //
378 InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
379 //
380 Pt = Value(theta);
1ef32e96 381 if(Abs(A)<0.0000001 || Abs(SigneSqrtDis)<0.0000000001) return(Standard_False);
7fd59977 382
383
384 //-- Approximation de la derivee (mieux que le calcul mathematique!)
385 Standard_Real dtheta = (DomainSup-DomainInf)*0.000001;
386 Standard_Real theta2 = theta+dtheta;
387 if((theta2<DomainInf) || ((theta2>DomainSup) && (!TwoCurves))
388 || (theta2>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
389 dtheta = -dtheta;
390 theta2 = theta+dtheta;
391 }
392 gp_Pnt P2 = Value(theta2);
393 dtheta = 1.0/dtheta;
394 Vec.SetCoord((P2.X()-Pt.X())*dtheta,
395 (P2.Y()-Pt.Y())*dtheta,
396 (P2.Z()-Pt.Z())*dtheta);
397
398 return(Standard_True);
399}
400//=======================================================================
401//function : FindParameter
402//purpose : Para est en sortie le parametre sur la courbe
403//=======================================================================
404 Standard_Boolean IntAna_Curve::FindParameter (const gp_Pnt& P,
405 Standard_Real& Para) const
406{
407 Standard_Real theta,z, aTolPrecision=0.0001;
c6541a0c 408 Standard_Real PIpPI = M_PI + M_PI;
7fd59977 409 //
410 switch (typequadric) {
411
412 case GeomAbs_Cylinder:
413 {
414 ElSLib::CylinderParameters(Ax3,RCyl,P,theta,z);
415 }
416 break;
417
418 case GeomAbs_Cone :
419 {
420 ElSLib::ConeParameters(Ax3,RCyl,Angle,P,theta,z);
421 }
422 break;
423
424 default:
841d6672 425 return Standard_False;
7fd59977 426 break;
427 }
428 //
429 Standard_Real epsAng = 1.e-8;
430 Standard_Real tmin = DomainInf;
431 Standard_Real tmax = DomainSup;
432 Standard_Real U,V,A,B,C,sint,cost,SigneSqrtDis;
433 Standard_Real z1,z2;
434
435 A=0.0; B=0.0; C=0.0;
436 U=0.0; V=0.0;
437 sint=0.0; cost=0.0;
438 SigneSqrtDis=0.0;
439 //U=V=A=B=C=sint=cost=SigneSqrtDis=0.0;
440 //
441 if (!firstbounded && tmin > theta && (tmin-theta) <= epsAng) {
442 theta = tmin;
443 }
444 else if (!lastbounded && theta > tmax && (theta-tmax) <= epsAng) {
445 theta = tmax;
446 }
447 //
448 if (theta < tmin ) {
449 theta = theta + PIpPI;
450 }
451 else if (theta > tmax) {
452 theta = theta - PIpPI;
453 }
454 if (theta < tmin || theta > tmax) {
455 if(theta>tmax) {
456 InternalUVValue(tmax,U,V,A,B,C,cost,sint,SigneSqrtDis);
457 gp_Pnt PMax(InternalValue(U,V));
458 if(PMax.Distance(P) < aTolPrecision) {
459 Para = tmax;
460 return(Standard_True);
461 }
462 }
463 if(theta<tmin) {
464 InternalUVValue(tmin,U,V,A,B,C,cost,sint,SigneSqrtDis);
465 gp_Pnt PMin(InternalValue(U,V));
466 if(PMin.Distance(P) < aTolPrecision) {
467 Para = tmin;
468 return(Standard_True);
469 }
470 }
471 //-- lbr le 14 Fev 96 : On teste malgre tout si le point n est pas le
472 //-- point de debut ou de fin
473 //-- cout<<"False 1 "<<endl;
474 // theta = tmin; le 25 Nov 96
475 }
476
477 if (TwoCurves) {
478 if(theta > tmax)
479 theta = tmax;
480 if(theta < tmin)
481 theta = tmin;
482 InternalUVValue(theta,U,z1,A,B,C,cost,sint,SigneSqrtDis);
483 A = B = C = sint = cost = SigneSqrtDis = 0.0;
484 InternalUVValue(tmax+tmax - theta,U,z2,A,B,C,cost,sint,SigneSqrtDis);
485
486 if (Abs(z-z1) <= Abs(z-z2)) {
487 Para = theta;
488 }
489 else {
490 Para = tmax+tmax - theta;
491 }
492 }
493 else {
494 Para = theta;
495 }
496
497 if((Para<DomainInf) || ((Para>DomainSup) && (!TwoCurves))
498 || (Para>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
499 return(Standard_False);
500 }
501
502 InternalUVValue(Para,U,V,A,B,C,cost,sint,SigneSqrtDis);
503 gp_Pnt PPara = InternalValue(U,V);
504 Standard_Real Dist = PPara.Distance(P);
505 if(Dist > aTolPrecision) {
506 //-- Il y a eu un probleme
507 //-- On teste si le point est un point double
508 InternalUVValue(tmin,U,V,A,B,C,cost,sint,SigneSqrtDis);
509 PPara = InternalValue(U,V);
510 Dist = PPara.Distance(P);
511 if(Dist <= aTolPrecision) {
512 Para = tmin;
513 return(Standard_True);
514 }
515
516 InternalUVValue(tmax,U,V,A,B,C,cost,sint,SigneSqrtDis);
517 PPara = InternalValue(U,V);
518 Dist = PPara.Distance(P);
519 if(Dist <= aTolPrecision) {
520 Para = tmax;
521 return(Standard_True);
522 }
523 if (TwoCurves) {
524 Standard_Real Theta = DomainSup+DomainSup-DomainInf;
525 InternalUVValue(Theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
526 PPara = InternalValue(U,V);
527 Dist = PPara.Distance(P);
528 if(Dist <= aTolPrecision) {
529 Para = Theta;
530 return(Standard_True);
531 }
532 }
533 return(Standard_False);
534 }
535 return(Standard_True);
536}
537//=======================================================================
538//function : InternalValue
539//purpose :
540//=======================================================================
541 gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
542 const Standard_Real _V) const
543{
544 //-- cout<<" ["<<U<<","<<V<<"]";
545 Standard_Real V = _V;
546 if(V > 100000.0 ) { V= 100000.0; }
547 if(V < -100000.0 ) { V=-100000.0; }
548
549 switch(typequadric) {
550 case GeomAbs_Cone:
551 {
552 //------------------------------------------------
553 //-- Parametrage : X = V * Cos(U) ---
554 //-- Y = V * Sin(U) ---
555 //-- Z = (V-RCyl) / Tan(SemiAngle)--
556 //------------------------------------------------
557 //-- Angle Vaut Cone.SemiAngle()
558 return(ElSLib::ConeValue(U,(V-RCyl)/Sin(Angle),Ax3,RCyl,Angle));
559 }
560 break;
561
562 case GeomAbs_Cylinder:
563 return(ElSLib::CylinderValue(U,V,Ax3,RCyl));
564 case GeomAbs_Sphere:
565 return(ElSLib::SphereValue(U,V,Ax3,RCyl));
566 default:
567 return(gp_Pnt(0.0,0.0,0.0));
568 }
569 return(gp_Pnt(0.0,0.0,0.0));
570}
571
572//
573//=======================================================================
574//function : SetDomain
575//purpose :
576//=======================================================================
577 void IntAna_Curve::SetDomain(const Standard_Real DInf,
578 const Standard_Real DSup)
579{
580 if(DInf>=DSup) {
581 Standard_DomainError::Raise("IntAna_Curve::Domain");
582 }
583 //
584 DomainInf=DInf;
585 DomainSup=DSup;
586}
Open CASCADE Technology repository